混合分布から乱数を生成 - 人工知能に関する断創録

Pythonによるモンテカルロ法入門（2014/6/20）のつづき。今回は、2.2.3節の混合分布から乱数を生成する実験してみました。2.2.3節によるとある確率分布fが別の確率分布gとpの混合分布として次のように書ける場合があるそうです。

$\displaystyle f(x) = \int_Y g(x|y) p(y) dy$

たとえば、スチューデントのt分布（f）は、カイ二乗分布（p）と正規分布（g）という二つの分布の混合で表せるという例が教科書で紹介されています。このような混合分布から乱数を生成するには次の手順に従えばよいとのこと。

分布 p(y) に従う乱数 y を生成し、
分布 g(x|y) に従う乱数 x を生成すると、
x は分布 f(x) に従う乱数になっている

教科書では、スチューデントのt分布と負の二項分布の例題があるのでさっそくPythonで検証してみます。

スチューデントのt分布

ここでは、パラメータを $\nu$ としたスチューデントのt分布 $T_\nu$ に従う乱数を生成したいとします。この場合、

カイ二乗分布 $\chi_\nu^2$ から乱数 y を生成し、
平均0、分散 $\nu/y$ である正規分布から乱数 x を生成すると、
x はスチューデントのt分布 $T_\nu$ に従う乱数になっている

本当かな？

# 混合分布によってカイ二乗乱数からスチューデントのt分布に従う乱数を生成
# X|y〜N(0, ν/y) && Y〜Xν^2 => X〜t(ν)
import numpy as np
from scipy.stats import chi2, norm, t
import matplotlib.pyplot as plt

# t分布のパラメータ
nu = 2

# サンプリング数
N = 100000

# カイ二乗分布に従う乱数を生成
Y = chi2(df=nu).rvs(size=N)

# 生成したカイ二乗分布の乱数をパラメータに持つ正規分布から乱数を生成
# scaleの次元数 = sizeとなっていないとエラーになる
# scaleの各値について1つの乱数を生成してくれる？
# nu=1のときに微妙に理論値とずれるがscipy.statsの乱数生成でも同じ結果なのでOK？
X = norm(loc=0, scale=np.sqrt(nu/Y)).rvs(size=N)

# 上の書き方の意図
# X = np.array([])
# for y in Y:
#     X = np.append(X, norm(loc=0, scale=np.sqrt(nu/y)).rvs(size=1))

# scipy.statsの乱数生成を使ったとき
# X = t(df=nu).rvs(N)

# t分布の裾野は広いので表示範囲を制限するため
X = X[(X>-5) & (X<5)]

plt.figure(1)
nbins = 100
plt.hist(X, nbins, normed=True)
rv = t(df=nu)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((-5, 5))

plt.show()

f:id:aidiary:20140709215106p:plain

混合分布から生成した乱数Xのヒストグラムはスチューデントのt分布の理論値と一致していることがわかります。

スクリプトのコメントにも書きましたが、正規分布のscaleパラメータに乱数の配列を渡しても機能するようです。もっとわかりやすく書きたいならYのforループにすればよいと思います。また、 $\nu=1$ のときにヒストグラムと理論値がずれるという結果になりました。これは、scipy.statsの乱数生成器を使っても同じ結果になったのですが結局なぜかわかりません。ヒストグラムの書き方が悪いのかな？

負の二項分布

次は、パラメータ n と p を持つ負の二項分布に従う乱数 X を生成してみます。負の二項分布はガンマ分布とポアソン分布の混合分布で表せるそうです。連続のガンマ分布と離散のポアソン分布を組み合わせるとは面白いな。

ガンマ分布 $Gamma(n, (1-p)/p)$ から乱数 y を生成し、
ポアソン分布 $P(y)$ から乱数 x を生成すると、
x は負の二項分布 $Neg(n, p)$ に従う乱数になっている

# 混合分布によってガンマ乱数から負の二項分布に従う乱数を生成
# Y〜Gamma(n,(1-p)/p) && X|y〜Poisson(y) => X〜Neg(n,p)
import numpy as np
from scipy.stats import gamma, poisson, nbinom
import matplotlib.pyplot as plt

# 負の二項分布のパラメータ
n = 6
p = 0.3

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

# サンプリング数
N = 100000

# ガンマ分布に従う乱数を生成
Y = gamma(a=n, scale=(1-p)/p).rvs(size=N)

# 生成したガンマ分布の乱数をパラメータに持つポアソン分布から乱数を生成
X = poisson(mu=Y).rvs(size=N)

# scipy.statsの機能で負の二項分布のパラメータ生成
# X = nbinom(n=n, p=p).rvs(N)

plt.figure(1)
weights = np.ones_like(X) / float(len(X))
nbins = np.arange(-0.5, K, 1.0)
plt.hist(X, nbins, weights=weights)
rv = nbinom(n=n, p=p)
t = np.arange(K)
plt.plot(t, rv.pmf(t), 'r-', lw=2)
plt.xlim((0, K))
plt.show()