n<-50000
w<-c()
x1<-c()
x2<-c()
z<-c()
for (i in
1:n) {
u <-
runif(1,-1,1)
v <-
runif(1,-1,1)
w[i] <-
u^2 + v^2
if(w[i]<=1){
z[i]<-(-2*log(w[i])/w[i])^0.5
x1[i]<-u*z[i]
x2[i]<-v*z[i]
}
else
if(w[i]>1){
x1[i]<-NA
x2[i]<-NA
}
}
y1<-x1[!is.na(x1)]
y2<-x2[!is.na(x2)]
par(mfrow=c(1,2))
hist(y1,freq=F,
breaks=100)
hist(y2,freq=F,
breaks=100)
library(MASS)
bivn <-
mvrnorm(50000, mu = c(0, 0), Sigma = matrix(c(1, 0, 0, 1), 2))
bivn.kde
<- kde2d(bivn[,1], bivn[,2], n = 50)
Marsaglia
<- kde2d(y1, y2, n = 50)
par(mfrow=c(3,2))
contour(Marsaglia,main="Marsaglia")
contour(bivn.kde,main="mvrnorm")
image(Marsaglia,main="Marsaglia")
image(bivn.kde,main="mvrnorm")
persp(Marsaglia,
phi = 0, theta = 0,main="Marsaglia")
persp(bivn.kde,
phi = 0, theta = 0, main="mvrnorm")
Compare the Marsaglia and mvrnorm:
