1. Generating Normally Distributed Values

This notebook shows how to generate values from a normal (Gaussian) distribution. There are three main algorithms that are covered.

  • Box-Muller transform

  • Marsaglia-Tsang

  • Ziggurat

The Box-Muller and Marsaglia-Tsang algorithms are very similar in shape and form. Box-Muller operates in polar coordinates while Marsaglia-Tsang operates in Cartesian coordinates. All three algorithms sample for points and test for rejection/acceptance. Box-Muller and Marsaglia-Tsang sample for points within the unit square while Ziggurat samples based on stacked rectangles covering half of the normal distribution. The Ziggurat algorithm is the fastest, followed by Marsaglia-Tsang and then Box-Muller.

We simply implement all three algorithms in this notebook and compare the outputs to numpy.

1.1. Box-Mueller method

[1]:
%matplotlib inline
import sys
import time
import ctypes
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

np.random.seed(37)
sns.set(color_codes=True)
num_samples = 10000

def box_muller(m, std):
    epsilon = sys.float_info.min
    two_pi = 2.0 * np.pi

    u1 = 0.0
    u2 = 0.0

    while True:
        u1 = np.random.rand()
        u2 = np.random.rand()
        if u1 > epsilon:
            break

    z0 = np.sqrt(-2.0 * np.log(u1)) * np.cos(two_pi * u2)
    z1 = np.sqrt(-2.0 * np.log(u1)) * np.sin(two_pi * u2)

    return m + std * z0

1.2. Marsaglia-Tsang method

[2]:
def marsaglia_tsang(m, std):
    u = 0.0
    v = 0.0
    s = 0.0

    while True:
        u = np.random.rand() * 2.0 - 1.0
        v = np.random.rand() * 2.0 - 1.0
        s = np.power(u, 2.0) + np.power(v, 2.0)
        if 0.0 < s and s <= 1.0:
            break

    z = np.sqrt(-2.0 * np.log(s) / s)
    return m + std * u * z

1.3. Ziggurat method

Although the Ziggurat algorithm is the fastest, it is not the easiest to understand or implement. The implementation here is a translation from JavaScript (which was itself a translation from C, dead link). Note all the low-level bitwise operations.

[3]:
class Ziggurat(object):
    def __init__(self, mean=0.0, std=1.0, seed=37):
        self.mean = mean
        self.std = std
        self.wn = np.array([0.0 for i in range(128)], dtype=np.float)
        self.fn = np.array([0.0 for i in range(128)], dtype=np.float)
        self.kn = np.array([0.0 for i in range(128)], dtype=np.float)

        self.jsr = 123456789 ^ (seed if seed is not None else int(time.time()))

        m1 = 2147483648.0
        dn = 3.442619855899
        tn = dn
        vn = 9.91256303526217e-3

        q = vn / np.exp(-0.5 * dn * dn)
        self.kn[0] = np.floor((dn / q) * m1)
        self.kn[1] = 0

        self.wn[0] = q / m1
        self.wn[127] = dn / m1

        self.fn[0] = 1.0
        self.fn[127] = np.exp(-0.5 * dn * dn)

        i = 126
        while i >= 1:
            dn = np.sqrt(-2.0 * np.log(vn / dn + np.exp(-0.5 * dn * dn)))
            self.kn[i+1] = np.floor((dn / tn) * m1)
            tn = dn
            self.fn[i] = np.exp(-0.5 * dn * dn)
            self.wn[i] = dn / m1

            i -= 1

    def shr3(self):
        def rshift(val, n):
            if val >= 0:
                return ctypes.c_int(val >> n).value
            else:
                return ctypes.c_int((val+0x100000000) >> n).value

        jz = self.jsr
        jzr = self.jsr
        jzr = ctypes.c_int(jzr ^ (jzr << 13)).value
        jzr = ctypes.c_int(jzr ^ rshift(jzr, 17)).value
        jzr = ctypes.c_int(jzr ^ (jzr << 5)).value
        self.jsr = jzr
        return ctypes.c_int((jz+jzr) | 0).value

    def uni(self):
        return 0.5 * (1.0 + self.shr3() / -np.power(2, 31))

    def uni_safe(self):
        while True:
            u = self.uni()
            if 0 != u:
                return u

    def rnorm(self):
        hz = self.shr3()
        iz = hz & 127
        return self.mean + self.std * (hz * self.wn[iz] if np.abs(hz) < self.kn[iz] else self.nfix(hz, iz))

    def nfix(self, hz, iz):
        r = 3.442619855899
        r1 = 1.0 / r
        x = 0
        y = 0

        while True:
            x = hz * self.wn[iz]
            if iz == 0:
                x = -np.log(self.uni_safe()) * r1
                y = -np.log(self.uni_safe())

                while True:
                    x = -np.log(self.uni_safe()) * r1
                    y = -np.log(self.uni_safe())
                    if y + y >= x * x:
                        break
                return r + x if hz > 0 else -r - x

            if self.fn[iz] + self.uni() * (self.fn[iz-1] - self.fn[iz]) < np.exp(-0.5 * x * x):
                return x

            hz = self.shr3()
            iz = hz & 127

            if np.abs(hz) < self.kn[iz]:
                return hz * self.wn[iz]

1.4. Sampling

[4]:
num_samples = 10000

mean = 0
std = 1
zigg = Ziggurat(mean=mean, std=std)

np_samples = np.random.normal(mean, std, num_samples)
bm_samples = np.array([box_muller(mean, std) for _ in range(num_samples)])
mt_samples = np.array([marsaglia_tsang(mean, std) for _ in range(num_samples)])
zi_samples = np.array([zigg.rnorm() for _ in range(num_samples)])

1.5. Visualize the densities

[5]:
fig, ax = plt.subplots(2, 2, figsize=(20, 20), sharex=False, sharey=False)

sns.kdeplot(np_samples, ax=ax[0, 0])
sns.kdeplot(bm_samples, ax=ax[0, 1])
sns.kdeplot(mt_samples, ax=ax[1, 0])
sns.kdeplot(zi_samples, ax=ax[1, 1])

ax[0, 0].set_title('Numpy')
ax[0, 1].set_title('Box-Muller')
ax[1, 0].set_title('Marsaglia-Tsang')
ax[1, 1].set_title('Ziggurat')
[5]:
Text(0.5, 1.0, 'Ziggurat')
_images/generate-gaussian-distributed-values_10_1.png