generate N(0,1) using uniform(0,1) in R

Question

I am trying to generate N(0,1) using uniform(0,1) for a simulation but can't get the code to run.

Firstly, my x is found by making X the subject for the CDF of normal followed by getting out the histogram. This is followed by imposing a normal curve to see if it fits. Below is my code.

sigma=1; mu=0
u<-runif(n)
x<-mu + sqrt(2*sigma(log(u*sigma*sqrt(2*pi))))
hist(x, Freq=F)
xpt<-seq(-5,5,0.1)
ypt<-dnorm(xpt,0,1)
lines(xpt,ypt,col=2)

Answer 1

You seem to invert PDF (probability density function) not CDF (cumulative density function).

Actually, Normal random variable is not generated with inverse CDF, as its CDF is not in closed form.

Check out Box-Muller transform. You simulate two independent sets of uniform random variables:

u <- runif(1000)
v <- runif(1000)
x <- sqrt(-2 * log(u)) * cos(2 * pi * v)
# y <- sqrt(-2 * log(u)) * sin(2 * pi * v)

Then x is from N(0, 1), and (x, y) are bivariate normal, with zero mean and identity covariance.

Answer 2

Source for all my comments:

Averill M. Law, W. David Kelton, Simulation Modeling and Analysis, third edition, McGraw-Hill, 2000. ISBN: 0-07-058290-4

The Box-Muller transform could be used by generating two U(0, 1) random variates.

`x1 = sqrt(-2 * log(u1)) * cos(2 * pi * u2)`
`x2 = sqrt(-2 * log(u1)) * sin(2 * pi * u2)`

To generalize:

box_muller <- function(n = 1, mean = 0, sd = 1)
{
  x <- vector("numeric", n)

  i <- 1
  while(i <= n)
  {
    u1 <- runif(1, 0, 1)
    u2 <- runif(1, 0, 1)

    x[i] <- sqrt(-2 * log(u1)) * cos(2 * pi * u2)

    if ((i + 1) <= n)
    {
      x[i + 1] <- sqrt(-2 * log(u1)) * sin(2 * pi * u2)
      i <- i + 1
    }

    i <- i + 1
  }

  x * sd + mean
}

hist(box_muller(1000))

Note, however, that

While this method is valid in principle, i.e. if U1 and U2 are truly IID U(0, 1) random variables, there is a serious difficulty if U1 and U2 are actually adjacent random numbers produced by a linear congruential generator. (p. 465)

Marsaglia Bray Polar Method

This method is a two-step process where

1. Generate U1 and U2 as IID U(0, 1); let V[i] = 2 * U[i] - 1 for i = 1, 2; and let W = V[1]^2 + V[2]^2.
2. If W > 1, go back to step 1. Otherwise, let Y = sqrt(-2 log(W)/W), X[1] = V[1] * Y, and X[2] = V2 * Y. Then X[1] and X[2] are IID N(0, 1) random variates. (p. 466)

Implementation:

marsaglia_bray <- function(n = 1, mean = 0, sd = 1)
{
  x <- vector("numeric", n)

  i <- 1

  while(i <= n)
  {
    u <- runif(2, 0, 1)
    v <- 2 * u - 1
    w <- sum(v^2)

    if (w < 1)
    {
      y <- sqrt(-2 * log(w) / w)

      z <- v * y

      x[i] <- z[1] 

      if ((i + 1) <= n)
      {
        x[i + 1] <- z[2]
        i <- i + 1
      }

      i <- i + 1
    }
  }

  x * sd + mean
}

hist(marsaglia_bray(1000))

You may also consider the Ziggurat Alogrithm