はじめに

　計算機シミュレーションの理論に関する勉強のアウトプットの一環として、二次元標準正規分布に従う疑似乱数を取得するためのプログラムをPythonで書いてみました。

　筆者はPython初心者ですので、コードについては復習の意味も込めて少し細かく説明を記載しています。

　また、プログラミングも計算機シミュレーションもまだまだ勉強中ですので、誤っている部分や改善点など多くあると思います。もし、何かお気づきのことがありましたら、気兼ねなくコメント等いただけますと幸いです。

　ちなみに、多変量正規分布に従う疑似乱数を生成する関数は、numpy（numpy.random.multivariate_normal）やscipy（scipy.stats.multivariate_normal）などに用意されているので、実践的には自分で書く必要はありません。

疑似乱数生成器のメソッド

　以下が、疑似乱数の生成をするメソッドです。

two-dimensionalStandardNormalDistributionPRNG.py

import matplotlib.pyplot as plt
import numpy as np
from numpy.random import *

def twoDimensionalStandadNormalDistribution(n):
	theta = rand(n)*2*np.pi
	u = rand(n)
	R = [(np.sqrt(-2*np.log(i))) for i in u]
	X = [i*np.cos(j) for (i, j) in zip(R, theta)]
	Y = [i*np.sin(j) for (i, j) in zip(R, theta)]

	return X, Y

各コードの簡単な説明

　中身をそれぞれ説明していきます。

　メソッド内の初めの三行にて、以下のコードを実行して、$r$と$\theta$の疑似乱数を取り出します。

theta = rand(n)*2*np.pi
u = rand(n)
R = [(np.sqrt(-2*np.log(i))) for i in u]

　rand()は、一様分布に従う$0$以上、$1$未満の疑似乱数を、引数に渡した数だけ返してくれます。

　thetaには、$0$以上、$2\pi$未満の疑似乱数を返して欲しいので、rand()で生成されたリストの要素全てに2*np.piを掛けます。

　Rには、$\sqrt{-2\log{U}}$（$U$は一様分布に従う$0$以上、$1$未満の疑似乱数）を返したいので、事前にrand()にて$n$個の疑似乱数を生成し、それらすべてを$U$に代入してRへ格納します。

X = [i*np.cos(j) for (i, j) in zip(R, theta)]
Y = [i*np.sin(j) for (i, j) in zip(R, theta)]

　最後は、上記のように、$x=r\sin{\theta},\quad y=r\cos{\theta}$に則って、$r$と$\theta$それぞれを極座標から$x$と$y$の直交座標に変換します。

　これで、return X, YにてXとYを返せば、二次元標準正規分布に従う疑似乱数のベクトル2つを得ることができます。

疑似乱数の生成結果を描画して確認

　実際に、生成した結果をヒストグラムに描画して確認してみます。

X, Y = twoDimensionalStandadNormalDistribution(100000)

fig = plt.figure()

ax1 = fig.add_subplot(1,2,1)
ax1.hist(X, bins=1000)
ax1.set_title('X')
ax1.set_xlabel('x')
ax1.set_ylabel('freq')

ax2 = fig.add_subplot(1,2,2)
ax2.hist(Y, bins=1000)
ax2.set_title('Y')
ax2.set_xlabel('y')
ax2.set_ylabel('freq')

fig.show()
fig.savefig("img.png")

各コードの簡単な説明

　中身をそれぞれ説明していきます。
　
　まず初めに、複数の返り値があるので、カンマで区切りながら、それぞれ代入します。

X, Y = twoDimensionalStandadNormalDistribution(100000)

　次に、描画領域全体を構成するFigureオブジェクトを用いて、以下のようにインスタンスを作成します。

fig = plt.figure()

　次に、ヒストグラムが描画される範囲を構成するAxesオブジェクトを用いて、それぞれインスタンスを作成します。

ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

　今回は、2つのヒストグラムを横並びに描画したいので、add_subplot('列数','行数','描画する行 or 列の番号')に、1つ目は(1,2,1)、2つ目は(1,2,2)と入れておきます。

　そして、以下にて、ヒストグラムを生成し、タイトルとラベルを設定します。

ax1.hist(X, bins=1000)
ax1.set_title('X')
ax1.set_xlabel('x')
ax1.set_ylabel('freq')

ax2.hist(Y, bins=1000)
ax2.set_title('Y')
ax2.set_xlabel('y')
ax2.set_ylabel('freq')

　bins=では区間・階級の数を指定します。

　最後に、作成したヒストグラムを描画します。

fig.show()
fig.savefig("img.png")

　ax1とax2それぞれのヒストグラムはfigの描画領域内に含まれているので、fig.show()としてあげます。

　fig.savefig("img.png")では、作成したヒストグラムを保存しています。環境によっては、fig.show()が上手く動作しないこともあるようなので、確認方法の一つとして記載しています。

　("img.png")の部分では、保存したいフォルダのパスを、保存したい画像名（img.pngの前）に加えることで、保存先のフォルダを指定できます。

生成された疑似乱数の分布の例

実際に描画すると、例えば以下のようになります。

　$x,y$のヒストグラムの形と$x$軸の値から、$X,Y$それぞれが二次元標準正規分布に従う疑似乱数のリストであることが確認できます。

おわりに（感想）

　これまで理論を中心に勉強していたので、実際にプログラムを書くことで、途中の計算や実際に得られる数値に対する解像度が上がりました。

　今回は、記事の中でそれぞれのコードが実行する計算の具体的な意味や理由の説明まで書けませんでした。

　なので次は、今回の疑似乱数生成における理論に関してまとめたものを投稿したいと思います。

参照

Voss, Jochen. An Introduction to Statistical Computing: A Simulation-Based Approach. First edition. Wiley Series in Computational Statistics. Chichester, West Sussex, UK: John Wiley & Sons, Inc, 2013.