Stan事始め
PyStan 3 試用記
ベイズ分析のプラットフォームとしてStanは使い易いと思われるが、StanのPython用ライブラリがPyStanであり、現在PyStan 3がリリースされている。PyStan 3をWindowsで利用するときは、WindowsにUbuntuをインストールして、Ubuntu上で利用する。
WindowsにUbuntuとPyStan 3をインストールする方法については、このウェブサイトを参照されたい。これらのインストール手順は、簡単である。
PyStan 3の使用例として、１変量データに正規分布モデルを適用して平均と標準偏差の母数パラメータ値の事後分布を求めるものと、２変量データに多変量正規分布モデルを適用して平均、標準偏差と相関係数を求めるものを作成してみた。、まず、１変量データの場合について説明する。その後、２変量データの場合について説明する。
ベイズ分析は、データからの情報を事後分布として得るものである。分布のサンプル値の代表値としてとして平均値がよく用いられるが、これは統計分析が手計算で行われていた時代におけるものと言え、PCで統計計算が行われる現代では一考されるべき事柄である。Gelmanら（2021）は、事後分布の代表値としてサンプルの変動に対する安定性に問題のある平均値ではなく中央値を薦めている。Gelmanら（2021）は、事後分布のモードは、データと事前分布に対して最良の全体的適合を与えるとしている（p. 107）。モード（MAP推定値）が精度の観点からよい点推定値であるという報告がある（岡本、2022a, b）。MAP（maximum a posteriory）推定値を求めるときは、モード付近のKDE（kernel density estimate）密度関数の形状が、山が崩れたような形ではなく、きれいな山形でなければならないが、このためにはサンプル数を多く取る必要があり、それに伴って計算時間も長くなる。ここでは、Gelmanら（2021）に倣って、事後分布の代表値として中央値を採用した。
事前分布として無情報事前分布を採用するときは、一様分布、あるいはそれに近いものとして分散の大きい正規分布あるいは指数分布を用いることができるが、PyStan 3では分散が大きすぎるという警告が出されるようになった。特に理由がない限り、事前分布として正規分布あるいは指数分布を用いるが（cf. Gelmanら（2021）, p. 124）、大きい分散の分布を設定すると、Stanから警告が提示されるので、まず、警告の提示されない程度の分散である分布を設定して、それをデータの分散に合わせてスケールするという方法を取った。
なお、岡本安晴（著）「いまさら聞けないPythonでデータ分析―多変量解析，ベイズ統計分析（PyStan，PyMC）―」でStan (PyStan 2)の説明を行っている。
１変量正規分布モデルに基づくベイズ分析スクリプトをリスト１のように用意した。リスト１は、Anaconda/Ubuntu/WSL/WindowsのPython 3.11仮想環境にPyStan 3をインストールして開発した。スクリプトは、WindowsにインストールしたPython IDLEのエディタを用いて編集した。すなわち、スクリプトファイルの編集はWindowsにインストールされたエディタで、スクリプトファイルの実行はWSL上の端末で行った。WindowsのエディタとWSLの端末を同時に開いておけば、スクリプトファイルの編集と実行がシームレスに行える。
スクリプトファイルは、圧縮ファイルSingleFilePS3.zipをダウンロード・展開すれば利用できる。
リスト１ 1変量正規分布によるベイズ分析スクリプト（OneVar.py）
import stan
import
matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import numpy as np
OneVar_code =
"""
data {
int<lower=0> J;
array[J] real y;
real
mu_p;
real
sgm_p;
}
parameters {
real
mu_t;
real<lower=0> sigma_t;
}
transformed
parameters {
real
mu;
real
sigma;
mu =
sgm_p * mu_t + mu_p;
sigma = sgm_p * sigma_t + 1.0e-7;
}
model {
mu_t
~ normal(0.0, 10.0);
sigma_t ~ exponential(0.1);
y ~
normal(mu, sigma);
}
"""
LData = [ 44, 44, 56, 59, 44, 44, 63, 60, 34,
58,
55, 50, 32, 61, 47, 25, 38, 47, 30, 42]
Data =
np.array(LData)
mu_p = Data.mean()
sgm_p = np.std(Data)
OneVar_dat = { 'J':
len(LData), 'y': LData, 'mu_p': mu_p, 'sgm_p': sgm_p}
sm =
stan.build(OneVar_code, data = OneVar_dat)
fit = sm.sample()
i_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry =
az.summary(i_data)
print('Summary
=\n',smry)
az.plot_trace(i_data,
var_names = ['mu', 'sigma'])
plt.show()
d_fm = fit.to_frame()
q1_mu, q2_mu, q3_mu =
np.percentile(d_fm['mu'], [25, 50, 75])
plt.title('Posterior
Distribution of $\mu$' +
f'\nQ1 = {q1_mu:.2f}, Med.= {q2_mu:.2f}, Q3 = {q3_mu:.2f}')
sb.kdeplot(d_fm['mu'])
plt.xlabel('$\mu$',
fontsize = 16)
plt.show()
q1_sgm, q2_sgm,
q3_sgm = np.percentile(d_fm['sigma'], [25, 50, 75])
plt.title('Posterior
Distribution of $\sigma$' +
f'\nQ1 = {q1_sgm:.2f}, Med.= {q2_sgm:.2f}, Q3 = {q3_sgm:.2f}')
sb.kdeplot(d_fm['sigma'])
plt.xlabel('$\sigma$',
fontsize = 16)
plt.show()
リスト１では、Stanスクリプトは変数OneVar_codeの文字列で与えられている。
配列は、予約語arrayを用いて
array[J] real y;
と記述されている。
データy[i]が正規分布に従うという記述
for (I in 1:N)
y[i] ~ normal(mu, sigma);
は、ベクトル化して
y ~ normal(mu,
sigma);
と記述されている。
正規分布のパラメータmuとmuは、まず、mu_tとsigma_tを
mu_t
~ normal(0.0, 10.0);
sigma_t ~ exponential(0.1);
と生成して、それらを次式
mu =
sgm_p * mu_t + mu_p;
sigma = sgm_p * sigma_t + 1.0e-7;
でスケーリングしている。
Stanスクリプトのビルドは、次のコードで行っている。
sm =
stan.build(OneVar_code, data = OneVar_dat)
ビルドに基づいたMCMCサンプリングは次のコードで行っている。
fit = sm.sample()
MCMCサンプリングの結果をInferenceData型に変換してサンプリングの統計量を表示することを以下のコードで行っている。
i_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary
=\n',smry)
MCMCサンプリングのトレース図を次のコードで行っている。
az.plot_trace(i_data,
var_names = ['mu', 'sigma'])
plt.show()
次のコードでは、MCMCサンプリングのサンプルがPandasのDataFrame型オブジェクトとして表される。
d_fm = fit.to_frame()
MCMCサンプルの統計量やグラフがPnadasのDataFrameに対するものとして求めることができる。
リスト１のスクリプトをAnaconda/Ubuntu/WSL/WindowsのPython 3.11仮想環境において次のコマンドで実行した。
(py311)
・・・/SingleFilePS3$ python OneVar.py
MCMCサンプリングが終わると、トレース図が表示される（図１）。
図１
図１のフォームを閉じると、パラメータ
の事後分布が表示される（図２）。
図２
図２のフォームを閉じると、パラメータ
の事後分布が表示される（図３）。
図３
図３のフォームを閉じると、プログラムの実行終了である。
端末画面には、MCMCサンプリングの統計量が以下のように表示されている。
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu_t 0.002 0.254 -0.481 0.491 0.006 0.005 2272.0 2091.0 1.00
sigma_t 1.099 0.197 0.769 1.465 0.005 0.004 1586.0 1576.0 1.01
mu
46.670 2.747 41.435 51.968 0.060 0.043 2272.0 2091.0 1.00
sigma 11.903 2.134 8.325 15.864 0.056 0.041 1586.0 1576.0 1.01
Stanスクリプトファイルを独立して作成する場合
リスト１では、StanスクリプトがPythonスクリプト内において文字列として与えられている。StanスクリプトをPythonスクリプトのファイルとは別のファイルとして用意した場合は、Stanスクリプトファイルを読み込んで文字列としてからbuildの第1引数として与えればよい。
いま、Stanスクリプトをリスト２のように用意したとする。これらのファイルは、圧縮ファイルMultiFilesPS3.zipとした。ダウンロード・展開して自由に利用されたい。
リスト２ Stanスクリプトファイル（OneVar.stan）
data {
int<lower=0> J;
array[J] real y;
real
mu_p;
real
sgm_p;
}
parameters {
real
mu_t;
real<lower=0> sigma_t;
}
transformed parameters {
real
mu;
real
sigma;
mu = sgm_p * mu_t +
mu_p;
sigma = sgm_p * sigma_t + 1.0e-7;
}
model {
mu_t
~ normal(0.0, 10.0);
sigma_t ~ exponential(0.1);
y ~
normal(mu, sigma);
}
リスト２のStanスクリプトを利用するPythonスクリプトをリスト３のように作成した。
リスト３ Stanスクリプトファイルを読み込むPythonスクリプト（OneVar_mf.py）
import stan
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import numpy as np
LData = [ 44, 44, 56, 59, 44, 44, 63, 60, 34,
58,
55, 50, 32, 61, 47, 25, 38, 47, 30, 42]
Data = np.array(LData)
mu_p = Data.mean()
sgm_p = np.std(Data)
OneVar_dat = { 'J': len(LData), 'y':
LData, 'mu_p': mu_p, 'sgm_p': sgm_p}
with open('OneVar.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = OneVar_dat)
fit = sm.sample()
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary =\n',smry)
az.plot_trace(i_data, var_names =
['mu', 'sigma'])
plt.show()
d_fm = fit.to_frame()
q1_mu, q2_mu, q3_mu =
np.percentile(d_fm['mu'], [25, 50, 75])
plt.title('Posterior Distribution of
$\mu$' +
f'\nQ1 = {q1_mu:.2f}, Med.= {q2_mu:.2f}, Q3 = {q3_mu:.2f}')
sb.kdeplot(d_fm['mu'])
plt.xlabel('$\mu$', fontsize = 16)
plt.show()
q1_sgm, q2_sgm, q3_sgm =
np.percentile(d_fm['sigma'], [25, 50, 75])
plt.title('Posterior Distribution of
$\sigma$' +
f'\nQ1 = {q1_sgm:.2f}, Med.= {q2_sgm:.2f}, Q3 = {q3_sgm:.2f}')
sb.kdeplot(d_fm['sigma'])
plt.xlabel('$\sigma$', fontsize = 16)
plt.show()
リスト３では、次のコードにより、リスト２のファイルを文字列として、stan.buildの第１引数としている。
with open('OneVar.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = OneVar_dat)
上のコード以外は、リスト１と同じである。
次のコマンドで実行すると、リスト１の場合と同じように実行が進む。
(py311) ・・・/MultiFilesPS3$ python
OneVar_mf.py
多変量正規分布の簡単な場合として、２変量正規分布を取り上げる。
ベイズ分析のためのStanスクリプトをリスト４のように作成した。これらのファイルは、まとめて圧縮ファイルBiVariateNormalPS3.zipとした。ダウンロード・展開して自由に利用されたい。
リスト４ ２変量正規分布モデルのStanスクリプト（TwoVarNormal.stan）
data {
int
N;
array[N]
vector[2] Data;
real
m1;
real
m2;
real
sd1;
real
sd2;
}
parameters {
real
p_mu1;
real
p_mu2;
real<lower
= 0.0> p_sgm1;
real<lower
= 0.0> p_sgm2;
real
p_r;
}
transformed parameters {
vector[2]
mu;
real
sgm1;
real
sgm2;
matrix[2,
2] cov;
real<lower=-1, upper=1> r;
r =
(inv_logit(p_r) - 0.5) * 2;
mu[1]
= p_mu1 * sd1 + m1;
mu[2]
= p_mu2 * sd2 + m2;
sgm1
= p_sgm1 * sd1 + 0.001;
sgm2
= p_sgm2 * sd2 + 0.001;
cov[1][1]
= sgm1 * sgm1;
cov[1][2]
= sgm1 * sgm2 * r;
cov[2][1]
= sgm1 * sgm2 * r;
cov[2][2]
= sgm2 * sgm2;
}
model {
p_mu1
~ normal(0.0, 10.0);
p_mu2
~ normal(0.0, 10.0);
p_sgm1
~ exponential(0.1);
p_sgm2
~ exponential(0.1);
p_r
~ normal(0.0, 5.0);
Data
~ multi_normal(mu, cov);
}
２変量データの配列は、予約語arrayを用いて次のように宣言されている。
array[N] vector[2] Data;
Dataの個々のData[i]の分布は、for文ではなく、ベクトル化されたコードで記述されている。
Data ~ multi_normal(mu, cov);
パラメータmuとcovを設定するパラメータsgm1とsgm2は、次のスケーリングにより与えられている。
mu[1]
= p_mu1 * sd1 + m1;
mu[2]
= p_mu2 * sd2 + m2;
sgm1
= p_sgm1 * sd1 + 0.001;
sgm2
= p_sgm2 * sd2 + 0.001;
スケーリングのもとのパラメータ値は次のコードで与えられ、これらはパラメータの分布の「分散が大きすぎる」という警告を避けるように設定されている。
p_mu1
~ normal(0.0, 10.0);
p_mu2
~ normal(0.0, 10.0);
p_sgm1
~ exponential(0.1);
p_sgm2
~ exponential(0.1);
相関係数rは、次のコードによるp_rの変換値として与えられている。
r =
(inv_logit(p_r) - 0.5) * 2;
相関係数rの事前分布を区間
上の一様分布で与えると、一様分布は数値微分で問題を起こす可能性があるという警告が出されるので、正規分布に従う値p_rを変換している。
リスト４のStanスクリプトを用いて２変量正規分布モデルによるベイズ分析を行うPythonスクリプト例をリスト５のように用意した。
リスト５ ２変量正規分布モデルによるベイズ分析Pythonスクリプト（BiVariateNormal.py）
import numpy as np
import stan
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] ]
mu = np.mean(Data, axis = 0)
print('mu =\n', mu)
cov = np.cov(np.array(Data).T)
print('cov =\n', cov)
sgms = [cov[0][0]**0.5,
cov[1][1]**0.5]
print('sgm =\n', sgms)
def f_init():
return {'p_mu1': 0.0, 'p_mu2': 0.0,
'p_sgm1': 1.0, 'p_sgm2': 1.0, 'p_r': 0.0}
stan_data = {'N': len(Data), 'Data':
Data, 'm1': mu[0], 'm2': mu[1],
'sd1': sgms[0], 'sd2': sgms[1]}
with open('TwoVarNormal.stan', 'r') as
fstan:
sm =
stan.build(fstan.read(), data = stan_data)
fit = sm.sample(num_chains = 4, init = [f_init()]*4)
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary =\n',smry)
az.plot_trace(i_data, var_names =
['mu', 'sgm1', 'sgm2', 'r'])
plt.show()
d_fm = fit.to_frame()
print('d_fm keys =\n', d_fm.keys())
samples_mu1 = d_fm['mu.1']
samples_mu2 = d_fm['mu.2']
med_mu1 = np.median(samples_mu1)
med_mu2 = np.median(samples_mu2)
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_mu1, label =
f'$\mu_1$ \nmed.={med_mu1:.2f}')
sb.kdeplot(samples_mu2, label =
f'$\mu_2$ \nmed.={med_mu2:.2f}')
plt.xlabel('')
plt.legend(fontsize = 14)
plt.show()
samples_sgm1 = d_fm['sgm1']
samples_sgm2 = d_fm['sgm2']
med_sgm1 = np.median(samples_sgm1)
med_sgm2 = np.median(samples_sgm2)
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_sgm1, label =
f'$\sigma_1$ \nmed.={med_sgm1:.2f}')
sb.kdeplot(samples_sgm2, label =
f'$\sigma_2$ \nmed.={med_sgm2:.2f}')
plt.xlabel('')
plt.legend(fontsize = 14)
plt.show()
samples_r = d_fm['r']
med_r = np.median(samples_r)
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_r, label = f'r
\nmed.={med_r:.2f}')
plt.xlabel('r', fontsize = 20)
plt.legend(fontsize = 14)
plt.show()
リスト５では、MCMCサンプリング時に次のコードでパラメータの初期値を指定している。
fit = sm.sample(num_chains = 4, init = [f_init()]*4)
これは、サンプリングにおいて不適切な値が生じた旨のInformation Messageが表示される回数を減らすためである。共分散行列のような制約の強いパラメータの場合、このようなメッセージが時々出るのはサンプリングがよい証拠であると説明されている。初期値の指定は、各マルコフ鎖ごとに設定するので、マルコフ鎖の数を４と明示して、それに合わせて初期値の指定を行っている。
リスト５のスクリプトを、Anaconda/Ubuntu/WSL/WindowsのPython 3.11仮想環境において次のコマンドで実行した。
(py311) ・・・/BiVariateNormalPS3$ python
BiVariateNormal.py
実行されると、MCMCサンプリング後、サンプリングのトレース図が表示される（図４）。
図４
図４のフォームを閉じると、平均値パラメータの事後分布のグラフが表示される（図５）。
図５
図５のフォームを閉じると、標準偏差パラメータの事後分布が表示される（図６）。
図６
図６のフォームを閉じると、相関係数の事後分布のグラフが表示される（図７）。
図７
図７のフォームを閉じると、スクリプトの実行終了である。
実行終了後の端末の画面には以下の出力がある。
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p_mu1 0.002 0.101 -0.180 0.198 0.002 0.002 4141.0 2613.0 1.0
p_mu2 0.002 0.102 -0.178 0.199 0.002 0.001 4054.0 3314.0 1.0
p_sgm1 1.018 0.074 0.880 1.153 0.001 0.001 4297.0 3055.0 1.0
p_sgm2 1.016 0.074 0.886 1.158 0.001 0.001 3969.0 3012.0 1.0
p_r
0.927 0.202 0.551 1.303 0.003 0.002 3512.0 2971.0 1.0
mu[0] 203.050 1.923 199.565 206.800 0.030 0.021 4141.0 2613.0 1.0
mu[1] 99.662 1.035 97.826 101.665 0.016 0.011 4054.0 3314.0 1.0
sgm1
19.477 1.414 16.848 22.072 0.022 0.015 4297.0 3055.0 1.0
sgm2 10.342 0.750 9.011 11.785 0.012 0.009 3969.0 3012.0 1.0
cov[0,0] 381.356 55.892 283.346 486.657 0.863 0.620 4297.0 3055.0 1.0
cov[0,1] 87.290 22.444 48.330 131.030 0.394 0.289 3349.0 2627.0 1.0
cov[1,0] 87.290 22.444 48.330 131.030 0.394 0.289 3349.0 2627.0 1.0
cov[1,1] 107.519 15.736 81.189 138.891 0.253 0.181 3969.0 3012.0 1.0
r
0.429 0.082 0.270 0.574 0.001 0.001 3512.0 2971.0 1.0
d_fm keys =
Index(['lp__', 'accept_stat__',
'stepsize__', 'treedepth__', 'n_leapfrog__',
'divergent__', 'energy__', 'p_mu1', 'p_mu2', 'p_sgm1', 'p_sgm2', 'p_r',
'mu.1',
'mu.2', 'sgm1', 'sgm2', 'cov.1.1', 'cov.2.1', 'cov.1.2',
'cov.2.2',
'r'],
dtype='object',
name='parameters')
前半の統計量の出力は、コード
i_data = az.from_pystan(posterior =
fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary =\n',smry)
によるものである。
後半の「d_fm keys =」に続く部分は、コード
d_fm = fit.to_frame()
print('d_fm keys =\n', d_fm.keys())
によるものである。パラメータcov[0,1]などが、PandasのDataFrameでは'cov.1.2'などのキーで表されていることが解る。InferenceData型では添え字は０から始まり、［ ］内にコンマ，で区切って並べられるが、DataFrame型では添え字は１から始まり、ピリオド．で繋いで並べられる。
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
岡本安晴（2019）いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan, PyMC）――．丸善出版
岡本安晴（2022a）ベイズ分析における点推定値例のRMSEとbiasによる比較――実験心理学の観点からの検討――．基礎心理学研究，41，1，1-7．https://doi.org/10.14947/psychono.41.1
岡本安晴（2022b）検出力指標d′の推定における実験上の問題．基礎心理学研究、41、2、107-114．https://doi.org/10.14947/psychono.41.16