Statistics 6401, Fall 2017

Project:

This R program transforms independent BVN \((X1,X2)\) to correlated BVN \((W1,W2)\).

Simulation:

Set up the parameters for the simulations. The goal is to simulate the the \(BVN(10, 25, 2, 3, -0.4)\)

n <- 2000  # number of values simulated
muw1 <- 10
muw2 <- 25
sigmaw1 <- 2
sigmaw2 <- 3
sigmasqw1 <- sigmaw1^2
sigmasqw2 <- sigmaw2^2
rhow1w2 <- -0.4
covw1w2 <- rhow1w2*sigmaw1*sigmaw2

Figure 1:

Plot the BVN

mu1 <- muw1 # setting the expected value of x1
mu2 <- muw2 # setting the expected value of x2
s11 <- sigmasqw1 # setting the variance of x1
s12 <- covw1w2 # setting the covariance between x1 and x2
s22 <- sigmasqw2 # setting the variance of x2
rho <- rhow1w2 # setting the correlation coefficient between x1 and x2
x1 <- seq(mu1-10,mu1+10,length=41) # generating the vector series x1 
x2 <- seq(mu2-10,mu2+10,length=41) 
f <- function(x1,x2){
    term1 <- 1/(2*pi*sqrt(s11*s22*(1-rho^2))) 
    term2 <- -1/(2*(1-rho^2))
    term3 <- (x1-mu1)^2/s11
    term4 <- (x2-mu2)^2/s22
    term5 <- -2*rho*((x1-mu1)*(x2-mu2))/(sqrt(s11)*sqrt(s22)) 
    term1*exp(term2*(term3+term4-term5)) 
} # setting up the function of the multivariate normal density
z <- outer(x1,x2,f) # calculating the density values 
persp(x1, x2, z, 
      main="Two dimensional Normal Distribution",
      sub=expression(italic(f)~(bold(x))==frac(1,2~pi~sqrt(sigma[11]~ 
                     sigma[22]~(1-rho^2)))~phantom(0)~exp~bgroup("{", 
                 list(-frac(1,2(1-rho^2)), 
                 bgroup("[", frac((x[1]~-~mu[1])^2, sigma[11])~-~2~rho~frac(x[1]~-~mu[1],
                 sqrt(sigma[11]))~ frac(x[2]~-~mu[2],sqrt(sigma[22]))~+~ 
                 frac((x[2]~-~mu[2])^2, sigma[22]),"]")),"}")),
      col="lightgreen", 
      theta=30, phi=20,
      r=50,
      d=0.1,
      expand=0.5,
      ltheta=90, lphi=180,
      shade=0.75,
      ticktype="detailed",
      nticks=5) # produces the 3-D plot
# adding a text line to the graph
mtext(expression(list(mu[1]==10,mu[2]==25,sigma[11]==2,sigma[22]==3,sigma[12]==-2.4,rho==-0.4)), side=3)

Figure 2

Simulate independent random uniforms

u1 <- runif(n)
u2 <- runif(n)
par(mfrow=c(2,2))
hist(u1)
hist(u2)

plot(u1,u2)

Check the means, standard deviation, and the correlation between the uniforms

cat('means')

means

umean = c(mean(u1),mean(u2)); umean

[1] 0.4993098 0.4953063

cat('standard deviations')

standard deviations

usd <-c(sd(u1),sd(u2)); usd

[1] 0.2899142 0.2876048

cat('correlation')

correlation

ucor <- cor(u1,u2); ucor

[1] 0.001807912

Figure 3

Simulate independent standard normals \(BVN(0,0,1,1,0)\) - Box-Muller Method

x1 <- sqrt(-2*log(u1))*cos(2*pi*u2)
x2 <- sqrt(-2*log(u1))*sin(2*pi*u2)
par(mfrow=c(2,2))
hist(x1)
hist(x2)

plot(x1,x2)

Check the means, standard deviation, and the correlation between the independent standard normals.

cat('means')

means

xmean <- c(mean(x1),mean(x2));xmean

[1] -0.0004314568  0.0223948236

cat('standard deviationss')

standard deviationss

xsd <- c(sd(x1),sd(x2));xsd

[1] 1.000200 1.009022

cat('correlations')

correlations

xcor <- cor(x1,x2);xcor

[1] 0.01658797

Figure 4

Transform to correlated normals \(BVN(0,0,\sigma_{W_1}^2,\sigma_{W_2}^2,\rho)\)

c11 <- sigmaw1
c21 <- rhow1w2*sigmaw2
c22 <- sigmaw2*sqrt(1-rhow1w2^2)
w1 <- c11*x1
w2 <- c21*x1 + c22*x2
par(mfrow=c(2,2))
hist(w1)
hist(w2)

plot(w1,w2)

Check the means, standard deviation, and the correlation between the correlations normals with mean zero.

cat('means')

means

wmean <- c(mean(w1),mean(w2))
wmean

[1] -0.0008629136  0.0620933326

cat('standard deviations')

standard deviations

wsd <- c(sd(w1),sd(w2))
wsd

[1] 2.00040 3.00452

cat('correlation')

correlation

wcor <- cor(w1,w2)
wcor

[1] -0.3841609

Figure 5

Transform to add means \(BVN(\mu_{W_1},\mu_{W_2},\sigma_{W_1}^2,\sigma_{W_2}^2,\rho)\)

y1 <- muw1 + w1
y2 <- muw2 + w2
par(mfrow=c(2,2))
hist(y1)
hist(y2)

plot(y1,y2)

Check the means, standard deviation, and the correlation between the correlations normals

cat('means')

means

ymean <- c(mean(y1),mean(y2))
ymean

[1]  9.999137 25.062093

cat('standard deviation')

standard deviation

ysd <- c(sd(y1),sd(y2))
ysd

[1] 2.00040 3.00452

cat('correlation')

correlation

ycor <- cor(y1,y2);ycor

[1] -0.3841609

Figure 6

Now rotate by \(\theta\) to make independent \(BVN(\mu_{W_1},\mu_{W_2},\sigma_{W_1}^2,\sigma_{W_2}^2,0)\) again

theta <- 0.5*atan((2*rhow1w2*sigmaw1*sigmaw2)/(sigmasqw1-sigmasqw2))
cat('theta')

theta

theta

[1] 0.3824964

v1 <- y1*cos(theta) + y2*sin(theta)
v2 <- -y1*sin(theta) + y2*cos(theta)
par(mfrow=c(2,2))
hist(v1)
hist(v2)

plot(v1,v2)

Check the means, standard deviation, and the correlation between the bivariate normals with zero correlation

cat('means')

means

vmean = c(mean(v1),mean(v2))
vmean

[1] 18.63067 19.51895

cat('standard deviation')

standard deviation

vsd = c(sd(v1),sd(v2))
vsd

[1] 1.761449 3.150561

cat('correlation')

correlation

vcor = cor(v1,v2)
vcor

[1] 0.01343562

Figure 7

Using the mvrnorm( ) function from the library MASS

library(MASS)
cat('means')

means

mu <- c(muw1,muw2);mu

[1] 10 25

cat('Variance-Covariance Matrix')

Variance-Covariance Matrix

Sigma <- matrix(c(sigmasqw1,covw1w2,covw1w2,sigmasqw2),2,2);Sigma

     [,1] [,2]
[1,]  4.0 -2.4
[2,] -2.4  9.0

aa <- mvrnorm(n, mu, Sigma)
par(mfrow=c(2,2))
hist(aa[,1])
hist(aa[,2])

plot(aa)

Check the means, standard deviation, and the correlation between the correlations normals

cat('means')

means

aamean <- c(mean(aa[,1]),mean(aa[,2]))
aamean

[1] 10.06245 24.96034

cat('standard deviations')

standard deviations

aasd <- c(sd(aa[,1]),sd(aa[,2]))
aasd

[1] 1.968776 2.921684

cat('correlation')

correlation

aacor <- cor(aa[,1],aa[,2]);aacor

[1] -0.3822073

Scatterplots for large data

If you do not have the following R packages installed, you should uncomment the following two lines in the next cell.

#install.packages("IDPmisc")
#install.packages("gplots")
library(IDPmisc)
library(gplots)

Figure 8

See Figure 2. Plot of independent uniforms

ipairs(matrix(c(u1,u2),n,2))

Figure 9

See Figure 3. Plot of independent standard normals

ipairs(matrix(c(x1,x2),n,2))

Figure 10

See Figure 4. Plot of correlated normals with mean zero

ipairs(matrix(c(w1,w2),n,2))

Figure 11

See Figure 5. Plot of correlated normals

ipairs(matrix(c(y1,y2),n,2))

Figure 12

See Figure 6. Plot of uncorrelated normals

ipairs(matrix(c(v1,v2),n,2))

Figure 13

See Figure 7. Plot of correlated normals

ipairs(aa)

Figure 14

The same plots but using the hist2d() function

par(mfrow=c(2,2))
hist2d(u1,u2, nbins=50, col = c("white",heat.colors(16))) 

----------------------------
2-D Histogram Object
----------------------------

Call: hist2d(x = u1, y = u2, nbins = 50, col = c("white", heat.colors(16)))

Number of data points:  2000 
Number of grid bins:  50 x 50 
X range: ( 0.0007641783 , 0.9997216 )
Y range: ( 0.0009787537 , 0.9985697 )

box()

hist2d(x1,x2, nbins=50, col = c("white",heat.colors(16)))

----------------------------
2-D Histogram Object
----------------------------

Call: hist2d(x = x1, y = x2, nbins = 50, col = c("white", heat.colors(16)))

Number of data points:  2000 
Number of grid bins:  50 x 50 
X range: ( -3.714062 , 3.299808 )
Y range: ( -3.75093 , 3.598531 )

box()

hist2d(w1,w2, nbins=50, col = c("white",heat.colors(16)))

----------------------------
2-D Histogram Object
----------------------------

Call: hist2d(x = w1, y = w2, nbins = 50, col = c("white", heat.colors(16)))

Number of data points:  2000 
Number of grid bins:  50 x 50 
X range: ( -7.428124 , 6.599615 )
Y range: ( -10.6697 , 10.338 )

box()

hist2d(y1,y2, nbins=50, col = c("white",heat.colors(16)))

----------------------------
2-D Histogram Object
----------------------------

Call: hist2d(x = y1, y = y2, nbins = 50, col = c("white", heat.colors(16)))

Number of data points:  2000 
Number of grid bins:  50 x 50 
X range: ( 2.571876 , 16.59962 )
Y range: ( 14.3303 , 35.338 )

box()