How can you sample from a gaussian (normal) distribution? Here’s an algorithm for generating a standard normal random variable Z:

Step 1: Generate a candidate output Y with an exponential distribution (with rate 1); that is, generate \(U_1\) (0,1) and set \(Y = -\ln{U}\) ; that means Y is in the range \((0,\infty)\)

Step 2: Generate $ U_2 $ (0,1). If \(U_2 {\leq} e^{-\frac {(Y-1)^2}{2}}\), set \(\lvert Z \rvert=Y\), otherwise go back to Step 1 and try again.

Step 3:Choose a sign: generate \(U_3\) (0,1) . Output Z if \(U_3\leq{0.5}\), else output \(-{Z}\) if \(U_3>0.5\) \(\)

\(\)

We would like to check the behaviour of \(P(Z{\leq}z)\), if it is indeed a normal distribution. Since z is generated only when the algorithm terminates (accepted), we can assert that: \(P(Z{\leq}z)= P(Z{\leq}z| Accept)\) Now, ‘accept’ means that: \(U_2{\leq}e^{-\frac {(Y-1)^2}{2}}\) and according the Bayes rule:

\[P(Z{\leq}z)= P(Z{\leq}z| Accept) = \frac {P(Z{\leq}z,U_2 \leq e^{-\frac {(Y-1)^2}{2}})}{P(U_2{\leq}e^{-\frac {(Y-1)^2}{2}})}\]

Let’s evaluate the denominator:

\[P(U_2{\leq}e^{-\frac {(Y-1)^2}{2}}) = \int_{u_1=0}^{1} P(U_2{\leq}e^{-\frac {(Y-1)^2}{2}} | u_1) du_1\]

Change variables, instead of \(du_1\) use y, so because \(y=-\ln{u_1}\), \(du_1 = -e^{-y}dy\), hence:

\[= \int_{y=\infty}^{0} P(U_2{\leq}e^{-\frac {(Y-1)^2}{2}} | Y=y) -e^{-y}dy\] \[= \int_{y=0}^{\infty} P(U_2{\leq}e^{-\frac {(Y-1)^2}{2}} | Y=y) e^{-y}dy = \int_{y=0}^{\infty} e^{-\frac {(y-1)^2}{2}} e^{-y}dy = \int_{y=0}^{\infty} e^{\frac {-y^2-1}{2}} dy = \sqrt{\frac{\pi}{2e}}\]

Now let’s deal with the nominator and first check what happens when z is positive, \(z \geq 0\) :

\[P(Z{\leq}z,U_2 \leq e^{-\frac {(Y-1)^2}{2}}) = ?\]

Let’s split the probablity to when Z is positive and when it is negative. When Z is negative it means we’ve output \(-Y\):

\[=P(Y{\leq}z,U_2 \leq e^{-\frac {(Y-1)^2}{2}} | Z>0)P(Z>0) +P(-Y{\leq}z,U_2 \leq e^{-\frac {(Y-1)^2}{2}} | Z\leq0)P(Z\leq0) =\]

Since \(P(Z)\leq0\) is 0.5, and since \(P(-Y\leq z) = 1\) :

\[=0.5 P(Y{\leq}z,U_2 \leq e^{-\frac {(Y-1)^2}{2}}) +0.5P(U_2 \leq e^{-\frac {(Y-1)^2}{2}})\] \[=0.5\int_{y=0}^{z} e^{-\frac {(y-1)^2}{2}} e^{-y}dy + 0.5\int_{y=0}^{\infty} e^{-\frac {(y-1)^2}{2}} e^{-y}dy\] \[=0.5 ( \int_{y=0}^{z} e^{\frac {-y^2-1}{2}} dy + \int_{y=0}^{\infty} e^{\frac {-y^2-1}{2}} dy ) = 0.5e^{-0.5} ( \int_{y=0}^{z} e^{\frac {-y^2}{2}} dy + \int_{y=0}^{\infty} e^{\frac {-y^2}{2}} dy )\]

Replacing the last term \(y \rightarrow -y\)

\[= 0.5e^{-0.5} ( \int_{y=0}^{z} e^{\frac {-y^2}{2}} dy + \int_{y=-\infty}^{0} e^{\frac {-y^2}{2}} dy ) = 0.5e^{-0.5} ( \int_{y=-\infty}^{z} e^{\frac {-y^2}{2}} dy )\]

Now, taking both the nominator and the denominator:

\[P(Z{\leq}z) = \frac{0.5e^{-0.5} \int_{y=-\infty}^{z} e^{\frac {-y^2}{2}} dy}{\sqrt{\frac{\pi}{2e}} = 0.5e^{-0.5} \sqrt{2\pi}} = \frac{1}{\sqrt{2\pi}} \int_{y=-\infty}^{z} e^{\frac {-y^2}{2}} dy\]

And this is exactly the definition of a standard normal distribution, with mean 0 and variance of 1. Solving for the case \(z\leq0\) is similar.