nykergoto’s blog

機械学習とpythonをメインに

Marsaglia法とBox-Muller法

一様分布からガウス分布を作るアルゴリズムとして有名なBox-Muller法というのがあります.
式としては、{z_1,z_2\in(0,1)}とした時に
{\displaystyle
y_1 = \sqrt{\mathstrut -2\ln z_1}\cos 2 \pi z_2 \\
y_2 = \sqrt{\mathstrut -2\ln z_1}\sin 2 \pi z_2
}
という変換を行う.これがBox-Muller
ボックス=ミュラー法 - Wikipedia


それと似たもので、一様分布する変数の範囲がちょっと変わって{z_1,z_2\in(-1,1)}としてやって、更に
{
z_1^2+z_2^2 \leq 1
}
を満たすものだけを取り出して(要するにそうじゃないやつを棄却して)
{
y_1 = z_1(-2 \ln r^2 / r^2)^{1/2}\\
y_2 = z_2(-2 \ln r^2 / r^2)^{1/2}
}
という変換をしてやるのが、Marsaglia
Marsaglia polar method - Wikipedia, the free encyclopedia

なんか読むとBox-Mullerのほうがcosとsinの計算を並列でできる?から早いとかなんとか.まあ棄却する部分があるMarsagliaのほうが遅そうやなぁという雰囲気はあるよね.

というわけで実装

import numpy as np
import math as math
import time
import matplotlib.pyplot as plt

pears = 1000
##compare marsaglia with box-muller
def boxmuller():
    """
    using box-muller method
    return:points random sampling following N(0,1)
    """
    points = []
    for i in range(pears):
        z1 = np.random.uniform(0,1)
        z2 = np.random.uniform(0,1)
        y1 = np.power((-2 * math.log(z1)),0.5) * math.cos(2 * np.pi * z2)
        y2 = np.power((-2 * math.log(z1)),0.5) * math.sin(2 * np.pi * z2)
        points.append([y1,y2])
    return points

def marsaglia():
    """
    marglia method
    return:points random sampling following N(0,1)
    """
    points = []
    for i in range(pears):
        z1,z2 = 0.,0.
        r = 0.
        while 1:
            z1 = np.random.uniform(-1,1)
            z2 = np.random.uniform(-1,1)
            r = z1 * z1 + z2 * z2
            if r <= 1:
                break
        y1 = z1 * np.power(-2. * math.log(r) / r,0.5)
        y2 = z2 * np.power(-2. * math.log(r) / r,0.5)
        points.append([y1,y2])
    return points

methods = [boxmuller,marsaglia]
for i in methods:
    start_time = time.time()
    for k in range(10):
        i()
    end_time = time.time()
    print (end_time-start_time)/100,i

points = np.array(marsaglia())
plt.plot(points[:,0],points[:,1],"ro",alpha=0.5)
points = np.array(boxmuller())
plt.plot(points[:,0],points[:,1],"go",alpha=0.5)
plt.show()

結果

0.00240999937057 <function boxmuller at 0x04357D70>
0.00396000146866 <function marsaglia at 0x04357DB0>

やっぱboxmullerのほうが早いみたい.

f:id:dette:20150802093017p:plain

どっちも標準正規分布にはなってるみたい.