# Generating the standard normals
par(mfrow = c(2, 3)) # landscape
z1=rnorm(150,0,1)
z2=rnorm(150,0,1)
# The parameters
m1=10
m2=15
s1=3
s2=2
rho=0.67

# Calculating the bivariate normal random numbers
x1<-s1*z1+m1
x2<-s2*(rho*z1+sqrt(1-rho^2)*z2)+m2

# Ploting the generated data
plot(x1, x2, xlim=c(-2,22), ylim=c(5,24), col="grey40", pch=20, xlab="x1", ylab="x2")
title(main=expression(mu[1]*"="*10*", "*mu[2]*"="*15*", "*sigma[1]^2*"="*9*", "*sigma[2]^2*"="*4*" and "*rho*"="*0.67*"." ))

#Overlaying the contours at levels=0.1, 0.5, 2, 5
x1=seq(0.5,19.5,0.1)
x2=seq(8,22,0.1)
z=outer(((x1-m1)/s1)^2,((x2-m2)/s2)^2,FUN="+")-2*rho*outer((x1-m1)/s1,(x2-m2)/s2,FUN="*")
contour(x1,x2,z,add=TRUE,level=c(0.1,2),lwd=2,col="blue")
contour(x1,x2,z,add=TRUE,level=c(0.5,5),lwd=2,col="red")

#----- Figure: x1, x2 have increased variance
z1=rnorm(150,0,1)
z2=rnorm(150,0,1)
# The parameters
m1=10
m2=15
s1=5
s2=4
rho=0.67

# Calculating the bivariate normal random numbers
x1<-s1*z1+m1
x2<-s2*(rho*z1+sqrt(1-rho^2)*z2)+m2

# Ploting the generated data
plot(x1, x2, xlim=c(-2,22), ylim=c(5,24), col="grey40", pch=20, xlab="x1", ylab="x2")
title(main=expression(mu[1]*"="*10*", "*mu[2]*"="*15*", "*sigma[1]^2*"="*25*", "*sigma[2]^2*"="*16*" and "*rho*"="*0.67*"." ))

#Overlaying the contours at levels=0.1, 0.5, 2
x1=seq(0,20,0.1)
x2=seq(6,24,0.1)
z=outer(((x1-m1)/s1)^2,((x2-m2)/s2)^2,FUN="+")-2*rho*outer((x1-m1)/s1,(x2-m2)/s2,FUN="*")
contour(x1,x2,z,add=TRUE,level=c(0.1,2),lwd=2,col="blue")
contour(x1,x2,z,add=TRUE,level=c(0.5),lwd=2,col="red")

#----- Figure: x1, x2 have decreased variance
z1=rnorm(150,0,1)
z2=rnorm(150,0,1)
# The parameters
m1=10
m2=15
s1=2
s2=1.5
rho=0.67

# Calculating the bivariate normal random numbers
x1<-s1*z1+m1
x2<-s2*(rho*z1+sqrt(1-rho^2)*z2)+m2

# Ploting the generated data
plot(x1, x2,xlim=c(-2,22), ylim=c(5,24), col="grey40", pch=20, xlab="x1", ylab="x2")
title(main=expression(mu[1]*"="*10*", "*mu[2]*"="*15*", "*sigma[1]^2*"="*4*", "*sigma[2]^2*"="*2.25*" and "*rho*"="*0.67*"." ))

#Overlaying the contours at levels=0.1, 0.5, 2
x1=seq(4,16,0.1)
x2=seq(11,20,0.1)
z=outer(((x1-m1)/s1)^2,((x2-m2)/s2)^2,FUN="+")-2*rho*outer((x1-m1)/s1,(x2-m2)/s2,FUN="*")
contour(x1,x2,z,add=TRUE,level=c(0.1,2),lwd=2,col="blue")
contour(x1,x2,z,add=TRUE,level=c(0.5),lwd=2,col="red")

# --- Figure: negative rho
z1=rnorm(150,0,1)
z2=rnorm(150,0,1)
# The parameters
m1=10
m2=15
s1=3
s2=2
rho=-0.67

# Calculating the bivariate normal random numbers
x1<-s1*z1+m1
x2<-s2*(rho*z1+sqrt(1-rho^2)*z2)+m2

# Ploting the generated data
plot(x1, x2, xlim=c(-2,22), ylim=c(5,24), col="grey40", pch=20, xlab="x1", ylab="x2")
title(main=expression(mu[1]*"="*10*", "*mu[2]*"="*15*", "*sigma[1]^2*"="*9*", "*sigma[2]^2*"="*4*" and "*rho*"="*-0.67*"." ))

#Overlaying the contours at levels=0.1, 0.5, 2, 5
x1=seq(0.5,19.5,0.1)
x2=seq(8,22,0.1)
z=outer(((x1-m1)/s1)^2,((x2-m2)/s2)^2,FUN="+")-2*rho*outer((x1-m1)/s1,(x2-m2)/s2,FUN="*")
contour(x1,x2,z,add=TRUE,level=c(0.1,2),lwd=2,col="blue")
contour(x1,x2,z,add=TRUE,level=c(0.5,5),lwd=2,col="red")

# ------- Figure, increase rho
z1=rnorm(150,0,1)
z2=rnorm(150,0,1)
# The parameters
m1=10
m2=15
s1=3
s2=2
rho=0.80

# Calculating the bivariate normal random numbers
x1<-s1*z1+m1
x2<-s2*(rho*z1+sqrt(1-rho^2)*z2)+m2

# Ploting the generated data
plot(x1, x2, xlim=c(-2,22), ylim=c(5,24), col="grey40", pch=20, xlab="x1", ylab="x2")
title(main=expression(mu[1]*"="*10*", "*mu[2]*"="*15*", "*sigma[1]^2*"="*9*", "*sigma[2]^2*"="*4*" and "*rho*"="*0.80*"." ))

#Overlaying the contours at levels=0.1, 0.5, 2
x1=seq(2,18,0.1)
x2=seq(9,20,0.1)
z=outer(((x1-m1)/s1)^2,((x2-m2)/s2)^2,FUN="+")-2*rho*outer((x1-m1)/s1,(x2-m2)/s2,FUN="*")
contour(x1,x2,z,add=TRUE,level=c(0.1,2),lwd=2,col="blue")
contour(x1,x2,z,add=TRUE,level=c(0.5),lwd=2,col="red")

#------- Figure: decrease rho
#par(mfrow = c(3, 2))
z1=rnorm(150,0,1)
z2=rnorm(150,0,1)
# The parameters
m1=10
m2=15
s1=3
s2=2
rho=0.20

# Calculating the bivariate normal random numbers
x1<-s1*z1+m1
x2<-s2*(rho*z1+sqrt(1-rho^2)*z2)+m2

# Ploting the generated data
plot(x1, x2, xlim=c(-2,22), ylim=c(5,24), col="grey40", pch=20, xlab="x1", ylab="x2")
title(main=expression(mu[1]*"="*10*", "*mu[2]*"="*15*", "*sigma[1]^2*"="*9*", "*sigma[2]^2*"="*4*" and "*rho*"="*0.20*"." ))

#Overlaying the contours at levels=0.1, 0.5, 2, 5
x1=seq(2,17,0.1)
x2=seq(10,21,0.1)
z=outer(((x1-m1)/s1)^2,((x2-m2)/s2)^2,FUN="+")-2*rho*outer((x1-m1)/s1,(x2-m2)/s2,FUN="*")
contour(x1,x2,z,add=TRUE,level=c(0.1,2),lwd=2,col="blue")
contour(x1,x2,z,add=TRUE,level=c(0.5,5),lwd=2,col="red")