まずいろいろ必要なものの import
%matplotlibのインライン指定もしておく
import numpy as np
import matplotlib.pyplot as plt
from pymc import Normal, deterministic, MCMC, Matplot
%matplotlib inline
$x$ を生成する正規分布の平均は $\mu$.この $\mu$ の事前分布は,平均が 0 で,分散が $1/0.01$. とりあえずどんな感じの分布になるかプロットしておく
mu = np.random.normal(0.0, 0.1, 100)
plt.hist(mu)
(array([ 4., 6., 7., 20., 19., 20., 12., 7., 3., 2.]), array([-0.24981899, -0.19842765, -0.14703632, -0.09564499, -0.04425366, 0.00713767, 0.05852901, 0.10992034, 0.16131167, 0.212703 , 0.26409433]), <a list of 10 Patch objects>)
x_sample
は学習用のサンプル.
平均が 1.0 で,標準偏差が 1.0 の正規分布に従う.
まずは少数のサンプル 5 個だけでやってみる.
$\mu$ が 1.0 を中心に分布するようになれば成功
x_sample = np.random.normal(loc=1.0, scale=1.0, size=5)
生成モデルは,非常に単純: $$ \mu\sim\mathrm{Normal}(0, 0.1^2) $$ $$ x\sim\mathrm{Normal}(\mu, 1.0^2) $$
※ NumPy の np.random.normal は平均と標準偏差で指定するが,PyMC の正規分布は,標準偏差ではなく精度(分散の逆数・逆行列)で指定する
$x$ は,平均が $\mu$ で,分散が $1/0.01$ としておく.この $x$ については標本が観測される(グラフィカルモデルでいえば黒丸)ので, value=x_sample
と observed=True
を指定しておく.
mu = Normal('mu', 0, 1 / (0.1 ** 2))
x = Normal('x', mu=mu, tau=1/(1.0**2), value=x_sample, observed=True)
マルコフ連鎖のモデルを指定.超初心者なので,最適化パラメータの意味とかは分かってないが,とりあえず, input
に関連する変数を入れておいて,適当な回数だけ反復してみる
M = MCMC(input=[mu, x])
M.sample(iter=10000)
[-----------------100%-----------------] 10000 of 10000 complete in 0.5 sec
左上はサンプリング中の $\mu$ の値,左下は自己相関で収束してるかどうかを見る.右が $\mu$ のデータを観測したあとのパラメータの事後分布
Matplot.plot(mu)
Plotting mu
今度は同じモデルで,サンプル数をガーンと1000個まで増やしてみる.
x_sample = np.random.normal(loc=1.0, scale=1.0, size=1000)
mu = Normal('mu', 0, 1 / (0.1 ** 2))
x = Normal('x', mu=mu, tau=1/(1.0**2), value=x_sample, observed=True)
M = MCMC(input=[mu, x])
M.sample(iter=10000)
[-----------------100%-----------------] 10000 of 10000 complete in 0.7 sec
すると $\mu$ の事後分布がずっと1.0の周囲で尖った分布になってる.すなわち,事前分布の影響が減ってデータに合わせた分布になってることが分かる.
Matplot.plot(mu)
Plotting mu