شنبه, ۱۷ شهریور ۱۳۹۷

آموزش توزیع نرمال دو متغیره در R

همانطور که می دانید نرم افزار R یکی از قدرتمندترین نرم افزارهای آماری است که پیدا کردن کدهای آن برای کارهای مختلف کار سختی نیست. یکی از کاربردهای R انجام اعمال مختلف بر روی توزیع نرمال دو متغیره (Bivariate Normal Distribution) است که در این آموزش به بحث شبیه سازی توزیع نرمال دو متغیره با R پرداخته ایم.

مدرس: افشین متولی
سایت: آر استودیو
سطح آموزش: متوسط

مقدمه

توزیع نرمال دو متغیره یکی از حوزه هایی است که اکثر دانشجویان در رسم کردن نمودارهای مربوط به آن در نرم افزار R مشکل دارند زیرا کمی گنگ است. در این آموزش نحوه رسم این نمودار را به صورت گام به گام به شما آموزش می دهیم. قبل از معرفی روش ها، مقادیر اولیه را به صورت فرضی برای توزیع تعریف می کنیم:

شبیه سازی نرمال دو متغیره

فرض کنید N تعداد نمونه های تصادفی مورد نیاز، mu1 میانگین متغیر تصادفی اول و mu2 میانگین متغیر تصادفی دوم و sigma ماتریس کوواریانس باشد.

library(mixtools)  #for ellipse
# r-studio.ir
N <- 200     # Number of random samples
set.seed(123)
# Target parameters for univariate normal distributions
rho <- -0.6
mu1 <- 1; s1 <- 2
mu2 <- 1; s2 <- 8

# Parameters for bivariate normal distribution
mu <- c(mu1,mu2) # Mean

sigma <- matrix(c(s1^2, s1*s2*rho, s1*s2*rho, s2^2),2) # Covariance matrix

library(mixtools) #for ellipse

# r-studio.ir

N <- 200 # Number of random samples

set.seed(123)

# Target parameters for univariate normal distributions

rho <- -0.6

mu1 <- 1; s1 <- 2

mu2 <- 1; s2 <- 8

# Parameters for bivariate normal distribution

mu <- c(mu1,mu2) # Mean

sigma <- matrix(c(s1^2, s1*s2*rho, s1*s2*rho, s2^2),2) # Covariance matrix

حال برای رسم کردن بیضی های مربوط به داده های نرمال دو متغیره از تابع زیر استفاده کنید:

ellipse_bvn <- function(bvn, alpha){
Xbar <- apply(bvn,2,mean)
S <- cov(bvn)
ellipse(Xbar, S, alpha = alpha, col="red")
}

ellipse_bvn <- function(bvn, alpha){

Xbar <- apply(bvn,2,mean)

S <- cov(bvn)

ellipse(Xbar, S, alpha = alpha, col="red")

}

برای شبیه سازی توزیع نرمال دو متغیره با R مراحل زیر را طی کنید:

#Install.pachages(MASS)
library(MASS)
bvn1 <- mvrnorm(N, mu = mu, Sigma = sigma ) # from MASS package
colnames(bvn1) <- c("bvn1_X1","bvn1_X2")
rbvn<-function (n, m1, s1, m2, s2, rho)
{
X1 <- rnorm(n, mu1, s1)
X2 <- rnorm(n, mu2 + (s2/s1) * rho * (X1 – mu1), sqrt((1 – rho^2)*s2^2))
cbind(X1, X2)
}
bvn3 <- rbvn(N,mu1,s1,mu2,s2,rho)
colnames(bvn3) <- c("bvn3_X1","bvn3_X2")
gibbs <- function (n, mu1, s1, mu2, s2, rho)
{
mat <- matrix(ncol = 2, nrow = n)
x <- 0
y <- 0
mat[1, ] <- c(x, y)
for (i in 2:n) {
x <- rnorm(1, mu1 + (s1/s2) * rho * (y – mu2), sqrt((1 – rho^2)*s1^2))
y <- rnorm(1, mu2 + (s2/s1) * rho * (x – mu1), sqrt((1 – rho^2)*s2^2))
mat[i, ] <- c(x, y)
}
mat
}
bvn4 <- gibbs(N,mu1,s1,mu2,s2,rho)
colnames(bvn4) <-  c("bvn4_X1","bvn4_X2")

library (mvtnorm)
bvn5 <- mvtnorm::rmvnorm(N,mu,sigma, method="svd")
colnames(bvn5) <- c("bvn5_X1","bvn5_X2")

bvn <- list(bvn1,bvn2,bvn3,bvn4,bvn5)

par(mfrow=c(3,2))
plot(bvn1, xlab="X1",ylab="X2",main= "All Samples")
for(i in 2:5){
points(bvn[[i]],col=i)
}
for(i in 1:5){
item <- paste("bvn",i,sep="")
plot(bvn[[i]],xlab="X1",ylab="X2",main=item, col=i)
ellipse_bvn(bvn[[i]],.5)
ellipse_bvn(bvn[[i]],.05)
}
par(mfrow=c(1,1))

#Install.pachages(MASS)

library(MASS)

bvn1 <- mvrnorm(N, mu = mu, Sigma = sigma ) # from MASS package

colnames(bvn1) <- c("bvn1_X1","bvn1_X2")

rbvn<-function (n, m1, s1, m2, s2, rho)

{

X1 <- rnorm(n, mu1, s1)

X2 <- rnorm(n, mu2 + (s2/s1) * rho * (X1 – mu1), sqrt((1 – rho^2)*s2^2))

cbind(X1, X2)

}

bvn3 <- rbvn(N,mu1,s1,mu2,s2,rho)

colnames(bvn3) <- c("bvn3_X1","bvn3_X2")

gibbs <- function (n, mu1, s1, mu2, s2, rho)

{

mat <- matrix(ncol = 2, nrow = n)

x <- 0

y <- 0

mat[1, ] <- c(x, y)

for (i in 2:n) {

x <- rnorm(1, mu1 + (s1/s2) * rho * (y – mu2), sqrt((1 – rho^2)*s1^2))

y <- rnorm(1, mu2 + (s2/s1) * rho * (x – mu1), sqrt((1 – rho^2)*s2^2))

mat[i, ] <- c(x, y)

}

mat

}

bvn4 <- gibbs(N,mu1,s1,mu2,s2,rho)

colnames(bvn4) <- c("bvn4_X1","bvn4_X2")

library (mvtnorm)

bvn5 <- mvtnorm::rmvnorm(N,mu,sigma, method="svd")

colnames(bvn5) <- c("bvn5_X1","bvn5_X2")

bvn <- list(bvn1,bvn2,bvn3,bvn4,bvn5)

par(mfrow=c(3,2))

plot(bvn1, xlab="X1",ylab="X2",main= "All Samples")

for(i in 2:5){

points(bvn[[i]],col=i)

}

for(i in 1:5){

item <- paste("bvn",i,sep="")

plot(bvn[[i]],xlab="X1",ylab="X2",main=item, col=i)

ellipse_bvn(bvn[[i]],.5)

ellipse_bvn(bvn[[i]],.05)

}

par(mfrow=c(1,1))

نتایج شبیه سازی توزیع نرمال دو متغیره با R:

            bvn2_X1       bvn2_X2
bvn2_X1     3.9               -9.5
bvn2_X2    -9.5               64.5

             bvn3_X1     bvn3_X2
bvn3_X1     4.1           -9.8
bvn3_X2    -9.8           63.7

               bvn4_X1       bvn4_X2
bvn4_X1      4.0              -9.7
bvn4_X2     -9.7               64.6

                bvn5_X1     bvn5_X2
bvn5_X1     4.0             -9.6
bvn5_X2    -9.6             65.3