読者です 読者をやめる 読者になる 読者になる

人工知能に関する断創録

人工知能、認知科学、心理学、ロボティクス、生物学などに興味を持っています。このブログでは人工知能のさまざまな分野について調査したことをまとめています。最近は、機械学習・Deep Learningに関する記事が多いです。



混合分布から乱数を生成

機械学習 PRML

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)という二つの分布の混合で表せるという例が教科書で紹介されています。このような混合分布から乱数を生成するには次の手順に従えばよいとのこと。

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

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

スチューデントのt分布

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

  1. カイ二乗分布 \chi_\nu^2 から乱数 y を生成し、
  2. 平均0、分散 \nu/y である正規分布から乱数 x を生成すると、
  3. 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 を生成してみます。負の二項分布はガンマ分布とポアソン分布の混合分布で表せるそうです。連続のガンマ分布と離散のポアソン分布を組み合わせるとは面白いな。

  1. ガンマ分布 Gamma(n, (1-p)/p) から乱数 y を生成し、
  2. ポアソン分布 P(y) から乱数 x を生成すると、
  3. 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()

f:id:aidiary:20140709215116p:plain

生成された乱数Xのヒストグラムは負の二項分布の理論値と一致していることがわかります。

Wikipediaを読むと確率分布間の関係などが細かに書いてあるのですがそのような理論的背景まで理解できていません。統計学をもう一回復習した方がいいかな。