Up

２項分布のベイズ分析

CmdStanPyによるStanの利用

拙著「いまさら聞けないPythonでデータ分析」では、PyStanを利用してStanによるベイズ分析を行った。ここでは、CmdStanPyによりStanによるベイズ分析を行う。

２項分布によるベイズ分析を行うためのStanスクリプトをリスト１のように用意した。ファイルはすべてbinfiles.zipまとめた。

リスト１　２項分布による分析のためのStanスクリプト（bin_model.stan）

data {

int k;

int n;

}

parameters {

real<lower=0, upper=1> p;

}

model {

p ~ beta(1,1);

k ~ binomial(n, p);

}

リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプトをリスト２のように用意した。

リスト２　リスト１のStanスクリプト（bin_model.stan）を用いてベイズ分析を行うPythonスクリプト（bin_model.py）

from cmdstanpy import CmdStanModel

import arviz as az

import matplotlib.pyplot as plt

model = CmdStanModel(stan_file="bin_model.stan") # コンパイル

fit = model.sample(data={'n':10, 'k':7}) # サンプリング

print(fit.summary())

inf_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

az.plot_trace(inf_data, figsize=(10,3)) # トレースプロット

plt.tight_layout()

plt.show()

次のコード

model = CmdStanModel(stan_file="bin_model.stan")

により、Stanスクリプトbin_model.stanのコンパイルを行い、次のコード

fit = model.sample(data={'n':10, 'k':7})

によって、データ{'n':10, 'k':7}を設定して、MCMCサンプリングを実行している。

サンプリングの統計を次のコードで出力している。

print(fit.summary())

Arviz（az）の関数を使ってトレース図を描くために、次のコードでサンプリングの結果fitをArviz用のデータの型InferrenceDataに変換している。

inf_data = az.from_cmdstanpy(fit)

変換したデータinf_dataを使って、次のコードでトレース図を描いている。

az.plot_trace(inf_data, figsize=(10,3))

リスト２のスクリプトを次のコマンドで実行してみる。

(stan) *****/binfiles$ python bin_model.py

コンパイルに多少の時間が掛かるが、しばらく待つと図１のトレース図が表示される。

図１

図１のWindowを、右上のX印のアイコンをクリックして閉じると、スクリプトの実行終了となる。

終了後の端末には、print(fit.summary())の出力が表示されている。

図２

Home