PythonでStanを使う
Python用のStan（CmdStanPy）を使ってみた。１変量データに正規分布モデルを適用して平均と標準偏差の母数パラメータ値の事後分布を求めるものと、２変量データに多変量正規分布モデルを適用して平均、標準偏差と相関係数を求めるものの２例を作成してみた。まず、１変量データの場合について説明する。その後、２変量データの場合について説明する。ファイルは、cmdstanpysamples.zipにまとめた。
CmdStanPyのインストールについては、このウェブサイト、あるいはこのウェブサイトなどで説明している。
なお、岡本安晴（著）「いまさら聞けないPythonでデータ分析―多変量解析，ベイズ統計分析（PyStan，PyMC）―」でStan(PyStan2)の説明を行っている。PyStan2のStanに対してCmdStanPyのStanは、若干の違いがある。PyStan2を使用するときは、このウェブサイト、あるいはこのウェブサイトを参照されたい。
１変量正規分布によるStanを用いたベイズ分析例を示す。
ベイズ分析のStanスクリプトをリストA.1のように作成した。
リストA.1 正規分布によるベイズ分析のStanスクリプト（OneVar.stan）
data {
int<lower=0> J;
array[J] real y;
// 要素数Jの1次元配列
}
parameters {
real<lower=-1000, upper=1000> mu;
real<lower=0.0001, upper=1000> sigma; // sigmaが０にならないようにする
}
model {
mu ~
uniform(-1000, 1000);
sigma ~ uniform(0.0001, 1000);
y ~
normal(mu, sigma);
}
データ宣言部data {…}においてデータが宣言される。
パラメータ宣言部において、正規分布の２つのパラメータ、平均muと標準偏差sigma、が宣言されている。パラメータは、十分に広い変域が設定されている。十分に広いかどうかは、分析後の事後分布で判断できる（図A.2、図A.3）。パラメータsigmaは、値0を取らないように、0より大きい値に制限されている。
モデル宣言部において、確率分布が設定されている。パラメータmuとsigmaは一様分布が設定されている。データｙは、平均mu、標準偏差sigmaの正規分布に従うという設定である。
リストA.1のStanスクリプトを用いてベイズ分析を行うPythonスクリプトを、リストA.2に示す。
リストA.2 リストA.1のStanスクリプトを用いるPythonスクリプト（OneVar.py）
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
import numpy
as np
Data = [ 44, 44, 56, 59, 44, 44, 63, 60, 34,
58,
55, 50, 32, 61, 47, 25, 38, 47, 30, 42]
OneVar_dat = { 'J' : len(Data),
'y': Data }
sm = CmdStanModel(stan_file='OneVar.stan') # ビルディング
fit = sm.sample(data=OneVar_dat)
# サンプリング
print('fit =\n', fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
az.plot_trace(inf_data)
# トレースプロット
plt.tight_layout()
plt.show()
dframe = fit.draws_pd()
# Pandas DataFrameへ変換
print('Keys...\n', dframe.keys())
sb.kdeplot(dframe['mu'])
plt.title(f"Med.
of mu = {np.median(dframe['mu']):.2f}",
fontsize=16)
plt.show()
sb.kdeplot(dframe['sigma'])
plt.title(f"Med.
of sigma = {np.median(dframe['sigma']):.2f}",
fontsize=16)
plt.show()
リストA.1とリストA.2のスクリプトファイルを同じフォルダに置き、リストA.2のPythonスクリプトファイルOneVar.pyを、CmdStanPyのインストールされた環境内で実行する。CmdStanPyのインストールについては、このウェブサイトなどで説明している。
カレントディレクトリをリストA.1とリストA.2の置かれているフォルダに移して、次のコマンドを実行する。
(stan) ****/cmdstanpysamples$ python OneVar.py
次のコードで、リストA.1のStanスクリプトがコンパイル・ビルドされる。
sm = CmdStanModel(stan_file='OneVar.stan') # ビルディング
ビルドされたモデルは、次のコードで実行され、MCMCサンプリングが行われる。
fit = sm.sample(data=OneVar_dat)
# サンプリング
サンプリングの統計が、次のコードで出力される。
print('fit =\n', fit.summary())
出力結果は、以下のようである。
fit =
Mean
MCSE StdDev MAD ...
ESS_bulk
ESS_tail
ESS_bulk/s R_hat
lp__ -50.1175 0.027823 1.09002 0.747015 ... 1733.27 2078.02 50978.5 1.00252
mu 46.6872 0.051761 2.67728 2.589010 ... 2732.62 2100.39 80371.1 1.00052
sigma 11.9363 0.048653 2.19670 1.984270 ... 2222.10 1916.64 65355.8 1.00369
次のスクリプトにより、MCMCサンプルがArvizのInferenceData型に変換されて、plot_traceによりトレース図が表示される（図A.1）。
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
az.plot_trace(inf_data)
# トレースプロット
plt.tight_layout()
plt.show()
図A.1
図A.1のWindowを閉じると、次のコードでMCMCサンプルがPandasのDataFrame型に変換される。
dframe = fit.draws_pd()
# Pandas DataFrameへ変換
DataFrameのキーが、次のコードで出力されている。
print('Keys...\n', dframe.keys())
上のコードの実行により、以下が出力される。
Keys...
Index(['chain__', 'iter__', 'draw__', 'lp__', 'accept_stat__', 'stepsize__',
'treedepth__', 'n_leapfrog__',
'divergent__', 'energy__', 'mu',
'sigma'],
dtype='str')
パラメータmuとsigmaのサンプルが、変数名の文字列をキーとして納められていることが分かる。
次のスクリプトにより、パラメータmuの事後分布のkdeplotによるグラフが表示される（図A.2）。
sb.kdeplot(dframe['mu'])
plt.title(f"Med.
of mu = {np.median(dframe['mu']):.2f}",
fontsize=16)
plt.show()
図A.2
図A.2のWindowを閉じると、次のスクリプトにより、パラメータsigmaの事後分布のグラフが描かれる（図A.3）。
sb.kdeplot(dframe['sigma'])
plt.title(f"Med.
of sigma = {np.median(dframe['sigma']):.2f}",
fontsize=16)
plt.show()
図A.3
図A.3のWindowを閉じると、スクリプトの実行終了である。
２変量正規分布によるベイズ分析のStanスクリプトを、リストB.1のように作成した。
リストB.1 ２変量正規分布によるベイズ分析（TwoVarNormal.stan）
data {
int
N;
array[N]
vector[2] Data; // ２次元ベクトルのＮ個の１次元配列
// ２変量正規分布に従うとする
}
parameters {
real
mu1;
real
mu2;
/*
共分散行列が計算可能な範囲になるように
標準偏差と相関係数の変域を制約する
*/
real<lower
= 0.0001, upper=1000.0> sgm1;
real<lower
= 0.0001, upper=1000.0> sgm2;
real<lower
= -0.9999, upper = 0.9999> r;
}
transformed parameters {
matrix[2,
2] cov;
cov[1][1] = sgm1 * sgm1;
cov[1][2] = sgm1 * sgm2 * r;
cov[2][1] = sgm1 * sgm2 * r;
cov[2][2] = sgm2 * sgm2;
}
model {
mu1
~ normal(0.0, 1000.0);
mu2
~ normal(0.0, 1000.0);
sgm1
~ uniform(0.0001, 1000.0);
sgm2
~ uniform(0.0001, 1000.0);
r
~ uniform(-0.9999, 0.9999);
for
(i in 1:N){
Data[i] ~ multi_normal([mu1,mu2], cov);
}
}
入力データDataは、２次元ベクトルのN個の配列である。
パラメータは、共分散行列の処理において値が大きすぎる、非正則であるなどの問題が生じないように、変域を制約している。
real<lower
= 0.0001, upper=1000.0> sgm1;
real<lower
= 0.0001, upper=1000.0> sgm2;
real<lower
= -0.9999, upper = 0.9999> r;
データDataの分布は、２変量正規分布に従うとして、次のようにおいている。
Data[i] ~ multi_normal([mu1,mu2], cov);
リストB.1のStanスクリプトを用いるPythonスクリプトを、リストB.2のように用意した（Sample2VarNorm.py）。
import numpy
as np
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import seaborn as sb
import arviz
as az
Data = [
[179, 121], [185, 89], [181, 106], [210,
100], [202, 104],
[197, 91], [227, 107], [208, 101], [188,
99], [189, 112],
[240, 111], [203, 103], [202, 85], [190,
95], [221, 96],
[150, 83], [205, 103], [177, 90], [221,
92], [172, 86],
[190, 87], [218, 101], [205, 100], [208,
120], [224, 98],
[199, 116], [203, 96], [213, 91], [200,
109], [190, 100],
[211, 115], [204, 85], [185, 92], [178,
86], [231, 122],
[206, 84], [197, 93], [236, 93], [213,
106], [205, 103],
[212, 102], [177, 86], [193, 96], [183,
93], [204, 87],
[203, 108], [203, 103], [180, 100], [196,
101], [167, 75],
[202, 106], [182, 98], [225, 109], [219,
115], [205, 102],
[200, 99], [219, 94], [189, 102], [220,
86], [186, 102],
[203, 101], [210, 124], [254, 103], [208,
95], [192, 106],
[176, 91], [173, 97], [195, 100], [189,
98], [220, 115],
[202, 112], [257, 108], [207, 104], [199,
96], [173, 98],
[230, 101], [230, 111], [212, 105], [196,
82], [183, 82],
[230, 103], [207, 112], [225, 101], [192,
104], [175, 83],
[223, 105], [192, 106], [202, 87], [204,
87], [184, 97],
[193, 94], [204, 109], [196, 92], [216,
109], [219, 113],
[256, 117], [210, 86], [217, 103], [206,
87], [213, 105] ]
stan_data = {'N': len(Data),
'Data': Data}
sm = CmdStanModel(stan_file='TwoVarNormal.stan') # ビルド
fit = sm.sample(data
= stan_data)
# サンプリング
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
az.plot_trace(inf_data, var_names=['mu1','mu2','sgm1','sgm2','r'],
figsize=(14,7))
# トレースプロット
plt.tight_layout()
plt.show()
dframe = fit.draws_pd()
# Pandas DataFrameへ変換
print('Keys...\n', dframe.keys())
p_pos_r = np.mean(dframe['r'] > 0.0) # Prop(r>0)の計算
sb.kdeplot(dframe['r'])
# rの事後分布
plt.xlabel('r', fontsize=16)
plt.title(f"Med.
of r = {np.median(dframe['r']):.3f}"
+
f"
P(r>0)={p_pos_r:.5f}", fontsize=18)
plt.show()
plt.figure(figsize=(14,8))
# 画面を２×２＝４の４分割にして描画
plt.subplot(221)
# 左上にmu1の事後分布を描画
sb.kdeplot(dframe['mu1'])
plt.xlabel('mu1', fontsize=14)
plt.title(f"Med.
of mu1 = {np.median(dframe['mu1']):.2f}",
fontsize=16)
plt.subplot(222)
# 右上にsgm1の事後分布を描画
sb.kdeplot(dframe['sgm1'])
plt.xlabel('sgm1', fontsize=14)
plt.title(f"Med.
of sgm1 = {np.median(dframe['sgm1']):.2f}",
fontsize=16)
plt.subplot(223)
# 左下にmu2の事後分布を描画
sb.kdeplot(dframe['mu2'])
plt.xlabel('mu2', fontsize=14)
plt.title(f"Med.
of mu2 = {np.median(dframe['mu2']):.2f}",
fontsize=16)
plt.subplot(224)
# 右下にsgm2の事後分布を描画
sb.kdeplot(dframe['sgm2'])
plt.xlabel('sgm2', fontsize=14)
plt.title(f"Med.
of sgm2 = {np.median(dframe['sgm2']):.2f}",
fontsize=16)
plt.tight_layout() # グラフの重なりを調整
plt.show()
# グラフの表示
リストB.1とリストB.2のファイルを同じフォルダに置き、CmdStanPyがインストールされている環境において、カレントディレクトリをリストB.1とB.2のファイルの置かれているフォルダに移して、次のコマンドを実行する。
(stan) ****/cmdstanpysamples$ python Sample2VarNorm.py
スクリプトファイルSample2VarNorm.pyが実行されると、次のコードでリストB.1のStanスクリプトのコンパイル・ビルドが行われる。
sm = CmdStanModel(stan_file='TwoVarNormal.stan') # ビルド
コンパイル・ビルドが済むと、次のコードでMCMCサンプリングが行われる。
fit = sm.sample(data
= stan_data)
# サンプリング
MCMCサンプリングの統計が、次のコードで出力される。
print(fit.summary())
出力は、以下のようである。
Mean
MCSE StdDev ... ESS_tail ESS_bulk/s R_hat
lp__ -613.963000 0.036220 1.559050 ... 2684.54 1855.76 1.000830
mu1
203.019000 0.029963 1.927020 ... 3333.09 4169.95 1.000450
mu2
99.640800 0.016816 1.038630 ... 2856.73 3844.25 1.001020
sgm1
19.460000 0.021169 1.398750 ... 3074.52 4404.76 1.004410
sgm2
10.345400 0.011723 0.751283 ... 3220.11 4232.36 1.001120
r
0.425591 0.001334 0.082328 ... 3055.01 3809.62 0.999764
cov[1,1] 380.649000 0.839835 55.162500 ...
3074.52
4404.76 1.004240
cov[1,2] 86.538400 0.391451 22.748100 ... 2691.65 3511.69 1.001030
cov[2,1] 86.538400 0.391451 22.748100 ... 2691.65 3511.69 1.001030
cov[2,2] 107.592000 0.248814 15.777300 ...
3220.11
4232.36 1.001210
次のスクリプトで、MCMCサンプルがArvizのInferenceData型に変換されて、Arvizのplot_traceによりトレース図が表示される（図B.1）。
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
az.plot_trace(inf_data, var_names=['mu1','mu2','sgm1','sgm2','r'],
figsize=(14,7))
# トレースプロット
plt.tight_layout()
plt.show()
図B.1
図B.1のWindowを閉じると、次のコードにより、MCMCサンプルがPandasのDataFrame型に変換される。
dframe = fit.draws_pd()
# Pandas DataFrameへ変換
DataFrame型のキーは、次のコードで出力される。
print('Keys...\n', dframe.keys())
出力は、以下のようである。
Keys...
Index(['chain__', 'iter__', 'draw__', 'lp__', 'accept_stat__', 'stepsize__',
'treedepth__', 'n_leapfrog__',
'divergent__', 'energy__', 'mu1', 'mu2',
'sgm1',
'sgm2', 'r', 'cov[1,1]', 'cov[2,1]',
'cov[1,2]', 'cov[2,2]'],
dtype='str')
パラメータ名の文字列がキーとなっていることが分かる。
次のスクリプトは、パラメータrの事後分布のグラフを描くものである（図B.2）。
p_pos_r = np.mean(dframe['r'] > 0.0) # Prop(r>0)の計算
sb.kdeplot(dframe['r'])
# rの事後分布
plt.xlabel('r', fontsize=16)
plt.title(f"Med.
of r = {np.median(dframe['r']):.3f}"
+
f"
P(r>0)={p_pos_r:.5f}", fontsize=18)
plt.show()
図B.2
図B.2のWindowを閉じると、以下のスクリプトが実行されて、平均値と標準偏差のパラメータの事後分布が表示される（図B.3）。
plt.figure(figsize=(14,8))
# 画面を２×２＝４の４分割にして描画
plt.subplot(221)
# 左上にmu1の事後分布を描画
sb.kdeplot(dframe['mu1'])
plt.xlabel('mu1', fontsize=14)
plt.title(f"Med.
of mu1 = {np.median(dframe['mu1']):.2f}",
fontsize=16)
plt.subplot(222)
# 右上にsgm1の事後分布を描画
sb.kdeplot(dframe['sgm1'])
plt.xlabel('sgm1', fontsize=14)
plt.title(f"Med.
of sgm1 = {np.median(dframe['sgm1']):.2f}",
fontsize=16)
plt.subplot(223)
# 左下にmu2の事後分布を描画
sb.kdeplot(dframe['mu2'])
plt.xlabel('mu2', fontsize=14)
plt.title(f"Med.
of mu2 = {np.median(dframe['mu2']):.2f}",
fontsize=16)
plt.subplot(224)
# 右下にsgm2の事後分布を描画
sb.kdeplot(dframe['sgm2'])
plt.xlabel('sgm2', fontsize=14)
plt.title(f"Med.
of sgm2 = {np.median(dframe['sgm2']):.2f}",
fontsize=16)
plt.tight_layout() # グラフの重なりを調整
plt.show()
# グラフの表示
図B.3
図B.3のWindowを閉じると、スクリプトの実行終了である。