離散的な乱数の生成 - 人工知能に関する断創録

Pythonによるモンテカルロ法入門（2014/6/20）のつづき。今回は、2.2.2節の離散分布に従う乱数を生成してみます。逆変換法の原理とほとんど同じだけど離散分布なので生成される乱数が自然数のみなのが違います。

$X \sim P_\theta$ を生成することを考える。ここで、 $P_\theta$ の台（確率変数の取る値の集合）を自然数とする。
$p_0 = P_\theta (X \le 0), \hspace{5mm} p_1 = P_\theta (X \le 1), \hspace{5mm} p_2 = P_\theta (X \le 2), \cdots$ を計算する。
$U \sim U_{[0,1]}$ を生成し、 $X=k \hspace{5mm} if \hspace{3mm} p_{k-1} < U < p_k$ とする。

これで、一様乱数Uから離散分布 $P_\theta$ に従う乱数Xが1つ生成されます。

ポアソン分布の例

まずはポアソン分布に従う乱数を生成してみます。

# ポアソン分布に従う乱数を累積分布関数から生成
import numpy as np
from scipy.stats import uniform, poisson
import matplotlib.pyplot as plt

# ポアソン分布のパラメータ
lam = 4

# この値まで確率値を計算
K = 20

# サンプリング数
N = 10000

rv = poisson(mu=lam)

# cdfの表を計算
t = np.arange(K)
prob = rv.cdf(t)

X = []
for i in range(N):
    u = uniform.rvs(loc=0, scale=1, size=1)
    # prob < uはcdfがuより小さいときTRUEを返す
    # TRUEは1と解釈されるのでsum()でTRUEの数をカウントして
    # インデックスを求めている
    X.append(np.sum(prob < u))

# ポアソン分布に従う乱数の分布を描画
# hist()のnormed=Trueはバーの積分が1になる確率密度関数になるため離散分布では使えない
# 離散分布ではバーの高さの合計が1になる確率質量関数にする必要がある
# http://stackoverflow.com/questions/3866520/plotting-histograms-whose-bar-heights-sum-to-1-in-matplotlib
plt.figure(1)
nbins = np.arange(-0.5, K, 1.0)
weights = np.ones_like(X) / float(len(X))
plt.hist(X, nbins, weights=weights)
plt.plot(t, rv.pmf(t), 'ro-', lw=1)
plt.xlim((0, K))

plt.show()

実行すると

f:id:aidiary:20140707213427p:plain

離散分布なので本来は赤丸のところにのみ値があります。ここでヒストグラムを書くときに要注意。今までは

plt.hist(X, nbins, normed=True)

のような感じで、Xに生成した乱数のリスト、nbinsにビンの数（ビンの幅）、normed=Trueで積分が1になるように正規化して表示してました。ただこれは連続分布にしか使えません。離散分布の場合は積分が1になるのではなく、バーの高さの和が1になるため下のように書く必要があります。

nbins = np.arange(-0.5, K, 1.0)
weights = np.ones_like(X) / float(len(X))
plt.hist(X, nbins, weights=weights)

参考にしたのはここ。

http://stackoverflow.com/questions/3866520/plotting-histograms-whose-bar-heights-sum-to-1-in-matplotlib

二項分布の例

もう一つ二項分布も試してみます。

# 二項分布に従う乱数を累積分布関数から生成
import numpy as np
from scipy.stats import uniform, binom
import matplotlib.pyplot as plt

# 二項分布のパラメータ
n = 20
p = 0.5

# この値まで確率値を計算
K = 40

# サンプリング数
N = 10000

rv = binom(n=n, p=p)

# cdfの表を計算
t = np.arange(K)
prob = rv.cdf(t)

X = []
for i in range(N):
    u = uniform.rvs(loc=0, scale=1, size=1)
    # prob < uはcdfがuより小さいときTRUEを返す
    # TRUEは1と解釈されるのでsum()でTRUEの数をカウントして
    # インデックスを求めている
    X.append(np.sum(prob < u))

# 二項分布に従う乱数の分布を描画
# hist()のnormed=Trueはバーの積分が1になる確率密度関数になるため離散分布では使えない
# 離散分布ではバーの高さの合計が1になる確率質量関数にする必要がある
# http://stackoverflow.com/questions/3866520/plotting-histograms-whose-bar-heights-sum-to-1-in-matplotlib
plt.figure(1)
nbins = np.arange(-0.5, K, 1.0)
weights = np.ones_like(X) / float(len(X))
plt.hist(X, nbins, weights=weights)
plt.plot(t, rv.pmf(t), 'ro-', lw=1)

plt.show()

f:id:aidiary:20140707213658p:plain

ちゃんと二項分布に従っています。