Robert & Casella. Chpt 2, Ex. 2.3

An antiquated generator of the normal distribution is:

Generate U_{1}, U_{2}, U_{3}, .... , U_{12} \sim \mathcal{U}[-1/2, 1/2]
Set Z=\sum_{i=1}^{12}U_{i}

The argument being that the CLT normality is sufficiently accurate with 12 terms

(a) Show that \mathbb{E}(Z)=0 and \sigma^{2}(Z)=1

For any random variable U_{i}\sim \mathcal{U}[-1/2, 1/2] it follows from the definition of expected value and variance of the uniform distribution that \mathbb{E}(U_{i})=0 and \sigma^{2}(U_{i})=1/12 . Since U_{1}, U_{2}, U_{3}, .... , U_{12} are i.i.d. then we can do:

\mathbb{E}(Z)=\sum_{i=1}^{12}\mathbb{E}(U_{i})=(12)(0)=0

\sigma^{2}(Z)= \sum_{i=1}^{12}\sigma^{2}(U_{i})=(12)(1/12)=1

(b) Using histograms, compare this CLT-normal generator with the Box-Muller algorithm.

(c)Compare both generators with rnorm

Here, I generated 10^{5} samples from three different normal generators and compared them by estimated PDFs, histograms and Q-Q plots. The samples generated from CLT-normal generator lacked fitness to the tail of N(0,1) while the Box-Muller transformation worked quite well compared with the samples generated from the R function rnorm

fig01

fig02

fig03

The pictures are kinda crappy as you can see, so here’s the R code for them in case you’d like to reproduce them:


#ex 2.3
#Box-Muller's algorithm
Box.Muller<-function(n){


U1<-runif(n)
U2<-runif(n)
sqrt(-2*log(U1))*cos(2*pi*U2)


}


#CLT normal generator
CLT.Norm<-function(n){

matrix(runif(n*12, -.5, .5), ncol=12) %*% rep(1,12)

}


set.seed(100)
n<-1e5
rv.bm<-Box.Muller(n)
rv.clt<-CLT.Norm(n)
rv.rnorm<-rnorm(n)


#Histograms and Q-Q plots
par(mfrow=c(2,3))
hist(rv.rnorm, 50, main="Histogram of rnorm")
hist(rv.bm, 50, main="Histogram of Box-Muller generator")
hist(rv.clt, 50, main="Histogram of CLT generator")
qqplot(rv.bm, rv.clt, pch=".", main="Q-Q Plot: Box-Muller vs CLT",
xlab="Box-Muller", ylab="CLT")
abline(0,1,col="red")
qqplot(rv.rnorm, rv.bm, pch=".", main="Q-Q Plot: rnorm vs Box-Muller",
xlab="rnorm", ylab="Box-Muller")
abline(0,1,col="red")
qqplot(rv.rnorm, rv.clt, pch=".", main="Q-Q Plot: rnorm vs CLT",
xlab="rnorm", ylab="Box-Muller")
abline(0,1,col="red")


#overlay estimated densities
par(mfrow=c(1,1),lwd=2)
plot(density(rv.rnorm), main="Estimated PDFs", col = "red")
lines(density(rv.bm), col="green")
lines(density(rv.clt), col="blue")
legend(-4,.4,legend=c("rnorm", "Box-Muller", "CLT"), lty=1, lwd=2, col=c("red","green","blue"))