一様乱数の生成

Pythonによるモンテカルロ法入門（2014/6/20）

最初にPythonで [0.0, 1.0] 間の一様乱数を生成してみます。

一様分布

一様分布（Uniform distribution）の確率密度関数（pdf）は下の式になります。

引用：http://en.wikipedia.org/wiki/Uniform_distribution_(continuous)

numpy.randomによる一様乱数の生成

numpy.random.uniformで一様分布に従う乱数が生成できます。下のスクリプトでは、下限を0.0、上限を1.0とした一様乱数をNsim個生成し、ヒストグラムを描画しています。

# 一様乱数の生成してヒストグラムを描画
import numpy as np
import matplotlib.pyplot as plt

np.random.seed()
N = 10000
x = np.random.uniform(0.0, 1.0, N)
nbins = 50
plt.hist(x, nbins, normed=True)
plt.show()

実行すると

numpy.random.RandomStateのページを見るとnumpyの乱数生成器には、最強の乱数生成器として名高いMersenne Twisterアルゴリズムが使われているとのこと。

scipy.statsによる一様乱数の生成

乱数を生成するだけだったらnumpy.randomモジュールを使えばよいのだけれど、確率密度関数や累積分布関数を使いたい場合は、scipy.stats.uniformの方が便利そうなのでこちらも試してみます。下の例では、一様分布に従う乱数を生成するとともに真の確率密度関数も赤線で描画しています。locが一様分布の下限値のパラメータaにあたり、scaleが上限値のパラメータbにあたります。

# 上と同じことをscipy.statsで実現
import numpy as np
from scipy.stats import uniform
import matplotlib.pyplot as plt

# 一様分布に従う確率分布からランダムサンプリング
np.random.seed()
N = 10000
# [0.0, 1.0]の一様分布に従う確率変数
rv = uniform(loc=0.0, scale=1.0)
# 一様分布からサンプリング
x = rv.rvs(size=N)
nbins = 50
plt.hist(x, nbins, normed=True)

# 真のPDFを描画
x = np.linspace(rv.ppf(0), rv.ppf(1), 100)
plt.plot(x, uniform.pdf(x), 'r-', lw=2, label='uniform pdf')

plt.show()

実行すると

scipy.statsの方は一度確率変数のオブジェクトを作成し、その確率変数のオブジェクトのメソッドを使うことで下のようないろいろな値を計算できるみたい（Scipy Statistics Tutorial）。

• rvs: Random Variates
• pdf: Probability Density Function
• cdf: Cumulative Distribution Function
• sf: Survival Function (1-CDF)
• ppf: Percent Point Function (Inverse of CDF)
• isf: Inverse Survival Function (Inverse of SF)
• stats: Return mean, variance, (Fisher’s) skew, or (Fisher’s) kurtosis
• moment: non-central moments of the distribution

逆に確率変数オブジェクト（rv）を使わずに下のように各メソッドにパラメータを直接渡すこともできました。同じ分布を何度も使う場合は確率変数オブジェクトを作った方がすっきりしそう。

# 同じことをscipy.statsで実現
# 確率変数オブジェクトを使わないで
import numpy as np
from scipy.stats import uniform
import matplotlib.pyplot as plt

# 一様分布に従う確率分布からランダムサンプリング
np.random.seed()
N = 10000
loc = 0.0
scale = 1.0

# 一様分布からサンプリング
x = uniform.rvs(loc=loc, scale=scale, size=N)
nbins = 50
plt.hist(x, nbins, normed=True)

# 真のPDFを描画
x = np.linspace(uniform.ppf(0, loc=loc, scale=scale), uniform.ppf(1, loc=loc, scale=scale), 100)
plt.plot(x, uniform.pdf(x, loc=loc, scale=scale), 'r-', lw=2, label='uniform pdf')

plt.show()

いろいろな方法があって混乱しそうなので今後はscipy.statsの確率変数オブジェクトを使う最初の書き方で進めていこうかな。