事後分布の分布
PyStanによるシミュレーション
ベイズ分析における事後分布はデータに依存して決まり、データは母集団からランダムにサンプリングされる。したがって、母集団が固定されていても、データはランダムなサンプルとしてサンプリングごとに確率的に変動し、データに依存している事後分布もサンプリングごとに変動する。これをシミュレーションによって調べてみた。母集団として、２次元多変量正規分布を想定した。各変量の平均は０、分散は１とし、相関係数
とした。この母集団から
個のサンプルをとり、ベイズ分析により
の事後分布を求めた。このサンプリングとベイズ分析を20回行い、20個の事後分布を得たものを図示したものが図１である。
図１ サンプル数100でサンプリングを20回行った。
ベイズ分析用のStanスクリプトをリスト１に示す。
リスト１ ２変量正規分布分析のためのStanスクリプト（ファイル名はリスト２での使用に合わせてstanBiNormal.stanとした）。
data{
int<lower=2>
N;
vector[2]
XY[N];
}
parameters{
real
mux;
real
muy;
real<lower=0>
sigmax;
real<lower=0>
sigmay;
real<lower
= -1, upper = 1> rho;
}
transformed parameters {
vector[2]
mu;
cov_matrix[2]
Sgm;
mu[1]
= mux;
mu[2]
= muy;
Sgm[1][1]
= sigmax * sigmax;
Sgm[1][2]
= rho * sigmax * sigmay;
Sgm[2][1]
= rho * sigmax * sigmay;
Sgm[2][2]
= sigmay * sigmay;
}
model{
for
(i in 1:N)
XY[i]
~ multi_normal(mu, Sgm);
}
このStanスクリプトをリスト２のPythonスクリプトで実行すると、図１のグラフを得る。Stanスクリプトのコンパイルは最初に行ってpickleによってファイルに保存し、MCMCによる事後分布からのサンプリングはpickleで読み込んだコンパイル済みのコードを使用して計算の高速化を図っている。PythonスクリプトはUbuntuで実行した。
リスト２ リスト１のStanスクリプトを使用して、事後分布の変動を調べるPythonスクリプト
import numpy as np
import numpy.random as nprn
import pystan
from pystan import StanModel
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
import time
rho = 0.5
n = 100
sm = StanModel(file =
'stanBiNormal.stan')
with open('model.pkl', 'wb') as f:
pickle.dump(sm,
f)
for i in range(20):
sm
= pickle.load(open('model.pkl', 'rb'))
xy
= nprn.multivariate_normal([0, 0], [[1, rho], [rho, 1]], n)
Data
= {'N': n, 'XY': xy}
print('\n\nSampling
started...i = {}\n'.format(i))
time.sleep(1)
def
f_init():
return
dict(mux = 0.0, muy = 0.0, rho = 0.0,
sigmax
= 1.0, sigmay = 1)
fit
= sm.sampling(data = Data, init = f_init)
DistRho
= fit['rho']
sb.kdeplot(DistRho,
clip = [-1, 1], alpha = 0.5, color = 'b')
plt.hlines([0], -1, 1, 'k')
plt.plot([rho], [0], 'ro', markersize
= 10)
plt.title('Distribution of \nPosterior
distributions of $\\rho$' +
',
n = {}'.format(n))
plt.show()
図１のグラフ（MCMCによるサンプリングの確率的変動のため、グラフは実行ごとに若干異なる）を見ると、事後分布とパラメータの真値（グラフでは赤小円で示されている）との関係は満足のいくものと言える。
母集団からのサンプルデータの数を個から
個に減らしてシミュレーションを行うと、図２の結果となった。
図２ サンプル数３でサンプリングを20回行った。
図２の事後分布のサンプリングごとの変動を見ると、サンプル数個では相関係数の推定は無理であると言える。