PyStanでベイズ推論+MCMC法によるデータ分析②(MCMCの結果と評価)
はじめに
こんにちはdeepblueでインターン生として働いている東です。
この記事は連載2記事目です。まだ読んでいない方は読んでいただけると嬉しいです。
PyStanでベイズ推論+MCMC法によるデータ分析③(単回帰モデル)
前回の記事と今回の概要
前回はPyStanの基本的な使い方、乱数の生成と収束確認を行いました。
この記事では参考図書第2部第5章を実装していきます。
本記事はMCMCによって得られた事後分布に従う乱数の評価、予測分布の考え方についてです。
今回もこちらの本を参考をしています。
トレースプロットと事後分布の可視化
前回とStanデータを引き続き使用しています。フォントサイズなど工夫しました。
fit = pystan.stan(file=stan,data=input_data, iter=2000, chains=4, thin=1, warmup=1000 , seed=1)
mcmc_sample = fit.extract(permuted=True)
plt.style.use('default')
fig, axes = plt.subplots(nrows=2, ncols=2,figsize=(30,10))
for i,p in zip(range(2),['mu','sigma']):
axes[i,0].set_title(p)
axes[i,0].grid()
sns.distplot(mcmc_sample[p] #事後分布の出力
,bins=11
,kde_kws={'label': 'Kernel density','color':'blue'}
,label=p+'_density'
,hist=False
,ax=axes[i,0]
)
line = axes[i,0].lines[0]
x = line.get_xydata()[:,0]
y = line.get_xydata()[:,1]
axes[i,0].fill_between(x,y, color="blue", alpha=0.3) #塗りつぶす
axes[i,0].set_title(p,fontsize=30) #フォントサイズ変更
for j in range(len(mcmc_sample_false[i])): #トレースプロット
sns.lineplot(data=mcmc_sample_false[:,j,i],label='chain'+str(j),ax=axes[i,1])
axes[i,1].set_title(p,fontsize=30)
axes[i,1].grid()
fig.show()
print(fit)
やはりpythonでの可視化は少し複雑で読みにくいコードとなってしまいました。
収束のチェックとしてn_effとRhatを確認します。
・n_eff : 有効なサンプルサイズ(100以上あることが望ましいです)
・Rhat : 収束指標(1.1未満であると良いです)
また、トレースプロットを可視化することでも視覚的に収束しているか確認できます。
収束確認(自己相関の可視化)
コレログラムを描くことでも乱数の収束を確認できます。
サンプリングした乱数の自己相関を求めるためにstatsmodelsというライブラリを使用しました。
import statsmodels.api as sm
fig, axes = plt.subplots(ncols=2,nrows=4, figsize=(18,14))
for i in range(4):
x=i*1000
y=(i+1)*1000
sm.graphics.tsa.plot_acf(mcmc_sample['mu'][x:y], lags=20, ax=axes[i,0],title='mu Autocorrelation chain:'+str(i+1))
sm.graphics.tsa.plot_acf(mcmc_sample['sigma'][x:y], lags=20, ax=axes[i,1],title='sigma Autocorrelation chain:'+str(i+1))
plt.show()
平均も標準偏差も自己相関が0に収束しています。
予測チェックとモデル選択
次に動物の発見個数のモデル化をし予測しています。今回は例として正規分布を仮定したときととポアソン分布を仮定したときの比較をしていきます。
データの読み込みをしstanファイルにデータを受け渡すためdict形式にいつもどおり変換します。
#データの読み込み
stan_normal=path+'2-5-1-normal-dist.stan'
stan_poisson=path+'2-5-2-poisson-dist.stan'
animal=pd.read_csv(path+'2-5-1-animal-num.csv')
#インプットするデータの作成
input_data = {'N': len(animal), 'animal_num':animal.animal_num.values}正規分布を仮定
#正規表現を仮定したモデル
fit_normal = pystan.stan(file=stan_normal,data=input_data, iter=2000, chains=4, thin=1, warmup=1000 , seed=1)
mcmc_sample_normal = fit_normal.extract(permuted=True)
#可視化
fig, axes = plt.subplots(ncols=3,nrows=2, figsize=(18,6))
axes[0,0].hist(animal.animal_num , bins=4 , rwidth=0.15 ,color='black')
axes[0,0].set_xlim([-2,4])
axes[0,0].set_title('animal_num')
for i in range(1,3):
axes[0,i].hist(mcmc_sample_normal['pred'][i-1] , bins=20 , rwidth=0.6 ,color='blue')
axes[0,i].set_ylim([0,100])
axes[0,i].set_xlim([-2,4])
axes[0,i].set_title('pred'+str(i-1))
for i in range(0,3):
axes[1,i].hist(mcmc_sample_normal['pred'][i+2] , bins=20 , rwidth=0.6 ,color='blue')
axes[1,i].set_ylim([0,100])
axes[1,i].set_xlim([-2,4])
axes[1,i].set_title('pred'+str(i+2))
左上の黒のグラフが観測値、青のグラフpred0から5は事後予測分布です。予測値の個体数が負の値になっており、正規分布を仮定したモデルでこのデータに対する予測をするのは無理がありそうです。
ポアソン分布を仮定
#ポアソン分布を仮定したモデル
fit_poisson = pystan.stan(file=stan_poisson,data=input_data, iter=2000, chains=4, thin=1, warmup=1000 , seed=1)
mcmc_sample_poisson = fit_poisson.extract(permuted=True)
#可視化
fig, axes = plt.subplots(ncols=3,nrows=2, figsize=(18,6))
axes[0,0].hist(animal.animal_num , bins=4 , rwidth=0.15 ,color='black')
axes[0,0].set_xlim([-2,4])
axes[0,0].set_title('animal_num')
for i in range(1,3):
axes[0,i].hist(mcmc_sample_poisson['pred'][i-1] , bins=20 , rwidth=0.6 ,color='blue')
axes[0,i].set_ylim([0,100])
axes[0,i].set_xlim([-2,4])
axes[0,i].set_title('pred'+str(i-1))
for i in range(0,3):
axes[1,i].hist(mcmc_sample_poisson['pred'][i+2] , bins=20 , rwidth=0.6 ,color='blue')
axes[1,i].set_ylim([0,100])
axes[1,i].set_xlim([-2,4])
axes[1,i].set_title('pred'+str(i+2))
左上の黒のグラフが観測値、青のグラフpred0から5は事後予測分布です。ポアソン分布を仮定したモデリングの方が正しく予測てきているように見えます。
まとめ
いかがだったでしょうか?可視化がメインになってしまいましたが、少しでもMCMCの理解が深まれば幸いです。次回は単回帰モデルについて実装していきます。
良かったら読んでみてください。
