[PYTHON] Random number generator with normal distribution N (0,1)

In this article, I will write the code to generate a variable that follows $ N (0,1) $ by the Box-Muller method.

Box-Muller Algorithm Overview

Let the two independent uniforms be $ U_1 and U_2 $. At this time, it is known that the following two variables defined by $ X_1 and X_2 $ follow a normal distribution (only the result or w).

\begin{eqnarray}
X_1 &=& \sqrt{-2\log(U_1)}\cos(2\pi U_2) \\
X_2 &=& \sqrt{-2\log(U_2)}\sin(2\pi U_1)
\end{eqnarray}

→ Pay attention to the "mixing condition" of $ U_1 and U_2 $ on the right side. (If you make a mistake here, the distribution will look strange.)

Preparation

This time we will use numpy, seaborn, so let's get ready.

## preparation
import numpy as np
import seaborn as sns

Step1: Generate uniform distribution

Use the function np.random.uniform (0,1) which takes a value between $ 0 \ sim 1 $. Since we need to create independent uniform distributions $ U_1 and U_2 $, we will put them in an array.

Let's set the number of samples to 10000.

## Sample size
N = 10000
## Uniform distributions
random_list1 = [np.random.uniform(0,1) for i in range(N)]
random_list2 = [np.random.uniform(0,1) for i in range(N)]

Step2: Confirm that it is uniform

For the time being, (just a visual check, but ...)

sns.tsplot(random_list1)

→ Somehow it looks okay!

Step3: Two variable bumps-zip function

As described in "Box-Muller Algorithm Overview", the "mixing condition" of $ U_1 and U_2 $ is troublesome, but let's use the zip function here.

zip(random_list1,random_list2))

If so, it will fetch one element at the same position in the two arrays random_list1 and random_list2 and create a new list. In terms of image

\begin{eqnarray}
list1 &=& [a_0, a_1, a_2, \ldots] \\
list2 &=& [b_0, b_1, b_2, \ldots]
\end{eqnarray}

Against

[[a_0,b_0], [a_1,b_1], [a_2,b_2] \ldots] ]

It is an image that returns.

Step4: Generate -map function

This time we will use the map function and the lambda expression.

norm_list = map(lambda x: np.sqrt(-2*np.log(x[0]))*np.sin(2*np.pi*x[1]), zip(random_list1,random_list2))

→ Try $ X_2 $, which is made by turning over $ U_1 and U_2 $!

Step5 ,, and what is the result?

sns.distplot(norm_list)

Somehow it looks like a normal distribution w

Step6 An easier way?

Well, do you say that everyone thinks together, isn't it prepared? !!

rnd = np.random.RandomState(0)
random_element = [rnd.randn() for i in range(N)]
sns.distplot(random_element)

Reference -Box-Muller method