Implementing Box-Muller Transform

After writing [the previous post, I think I understand Box-Muller transform.

I’m implementing it in this post.

Original form of Box-Muller transform is as follows:

$U_1, U_2$ : random variables from uniform distribution (0,1)
$y_1 \, = \, \sqrt{-2 ln U_1} cos(2\pi U_2)$
$y_2 \, = \, \sqrt{-2 ln U_1} sin(2\pi U_2)$

However it includes square root, logarithm, trigonometric functions(sin/cos), which are costly and complex, might includes errors after calculation. J. Bell proposed a more robust and fast way of sampling, which is called polar form.

Unlike original form, polar form is a rejection sampling.

Draw samples ( $x_1, x_2$ ) from uniform distribution(-1, 1)
Calculate $s \, = x_1 * x_1 \, + x_2 * x_2$
If $s \lt 1$ , go to 4, else go to 1
Get samples $y_1, y_2$ by calculating: $y_1 = x_1 \, \sqrt{\frac{-2 \ln s}{s}}$ $y_2 = x_2 \, \sqrt{\frac{-2 \ln s}{s}}$

Python implementation

import random as rd
import math as ma
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

data = []
for i in range(0, 500):
    while True:
        x1 = rd.uniform(-1,1);
        x2 = rd.uniform(-1,1);
        s = x1 ** 2 + x2 **2
        if s < 1.0:
            break
    w = ma.sqrt((-2.0 *  ma.log(s)) / s)
    y1 = x1 * w
    y2 = x2 * w
    data.append(y1)
    data.append(y2)

sns.distplot(data)
plt.show()

The result:

f:id:nakaly:20170213050259p:plain:w500