Marsaglia法とBox-Muller法
一様分布からガウス分布を作るアルゴリズムとして有名なBox-Muller法というのがあります.
式としては、とした時に
という変換を行う.これがBox-Muller
ボックス=ミュラー法 - Wikipedia
それと似たもので、一様分布する変数の範囲がちょっと変わってとしてやって、更に
を満たすものだけを取り出して(要するにそうじゃないやつを棄却して)
という変換をしてやるのが、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のほうが早いみたい.
どっちも標準正規分布にはなってるみたい.