この記事ではBox-Mullerの方法で$N(0,1)$に従う変数を生成するコードを書いてみます。
Box-Mullerアルゴリズム概要
2つの独立な一様分布を$U_1, U_2$とします。このとき次の$X_1, X_2$で定義される2 変数は正規分布に従うことが知られています（って結果だけかヲイ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}
→右辺の$U_1, U_2$の「混ざり具合」に注意ですねー。
（ココを間違えると、変な感じの分布になりますw）
準備
今回はnumpy,seabornを使うので、準備しときましょ。
## preparation
import numpy as np
import seaborn as sns
Step1: 一様分布を生成
$0\sim 1$の間に値をとる関数np.random.uniform(0,1)を使います。
独立な一様分布$U_1, U_2$を作る必要があるので、配列にもたせておきます。
一端、sample数を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: 一様であることの確認
一応、、、（ただの目視確認ですが、、）
sns.tsplot(random_list1)
→何となく大丈夫そう！
Step3: 2変数のぶっこみ - zip関数
「Box-Mullerアルゴリズム概要」に記載したとおり、$U_1, U_2$の「混ざり具合」が面倒ですが、ここはzip関数を使ってみましょうか。
zip(random_list1,random_list2))
とすれば、random_list1とrandom_list2の2つの配列の同じ位置にある要素を1つずつ取ってきて新しいリストを作ります。イメージ的には
\begin{eqnarray}
list1 &=& [a_0, a_1, a_2, \ldots] \\
list2 &=& [b_0, b_1, b_2, \ldots]
\end{eqnarray}
に対して
[[a_0,b_0], [a_1,b_1], [a_2,b_2] \ldots] ]
を返すイメージです。
Step4: 生成しまーす -map関数とか
今回はmap関数とlambda式を使いますー。
norm_list = map(lambda x: np.sqrt(-2*np.log(x[0]))*np.sin(2*np.pi*x[1]), zip(random_list1,random_list2))
→$U_1, U_2$をひっくり返して作った$X_2$も挑戦してみて下さい！
Step5 、、で結果は？
sns.distplot(norm_list)
何となく正規分布っぽいですねｗ
Step6 もっと簡単な方法？
まぁ、考えることはみんな一緒と言いますか、用意されてますよね、そりゃ。!
rnd = np.random.RandomState(0)
random_element = [rnd.randn() for i in range(N)]
sns.distplot(random_element)
