PyStan 3のためのStan入門
Stanの解説はウェブサイト
https://mc-stan.org/docs/reference-manual/
などで入手できる。
Stanの基本構文は
<program> ::= [<function_block>] [<data_block>] [<transformed_data_block>]
[<parameters_block>] [<transformed_parameters_block>]
[<model_block>] [<generated_quantities_block>] EOF
である（https://mc-stan.org/docs/reference-manual/syntax.html）。
すなわち、Stanプログラムは、
［functionブロック］ ［dataブロック］ ［transformed dataブロック］
［parametersブロック］ ［transformed parametersブロック］
［modelブロック］ ［generated quantitiesブロック］
からなり、これらの順に記述する。必要がないブロックは省略する。
以下、簡単な例で説明をすすめる。なお、StanによるＭＣＭＣサンプリングはMiniconda/Ubuntu24/WSL/Windows10のPython3.12仮想環境で行い、ＭＣＭＣサンプリングで得られた事後分布の描画などは、Windows10にインストールしたUbuntuではできなかったので、Windows10にインストールしたPython3.13で行った。Linuxで実行する場合は、同じ開発環境で続けて実行すればよい。WindowsにUbuntu（Linux）をインストールする方法は、このウェブサイトで説明した。WindowsとUbuntuのソフトを並行して利用できる。
ベイズ分析の解説書としてGelman et al. 2014がある。
ファイルはまとめて圧縮ファイルintropystan3files.zipとした。
以下に、2項分布、正規分布、２変量正規分布を例として説明する。
コイン投げをＮ試行行い、表の出た試行数がｋであったときの２項分布によるStanスクリプトとPythonスクリプトをリスト1.1のように書いてみた。
import stan # pystan3がimportされる
import pickle # データの保存用モジュール
StanCode = """
data {
int
N; // 2項分布の総試行数
int
k; // 生起回数
}
parameters {
real<lower=0,upper=1> p; // 生起確率
}
model {
k ~
binomial(N, p);
// 2項分布のモデル
}
"""
sm = stan.build(StanCode,
data={'N':10, 'k':7})
fit = sm.sample()
df_fit = fit.to_frame()
# PandasのDataFrameに変換
print(df_fit['p'].describe())
# DataFameの関数describeの使用
with open('param_p.pkl', 'wb') as f:
pickle.dump(df_fit['p'],
f)
# パラメータpのMCMCサンプルの保存
PyStan3のインストールは
pip install
pystan
であるが、importは
import stan
である。
Stanスクリプトは、変数StanCodeで指し示されている文字列で与えられている。上のStanスクリプトでは、dataブロック、parametersブロックとmodelブロックが使われており、各ブロックの内容は、{と}で囲まれた部分である。各ブロックは文からなるが、各文の終わりにはセミコロン;を置く。{や;は、半角文字であることに注意。Stanスクリプトは、コメントを除き、半角英数字で書く（https://mc-stan.org/docs/reference-manual/encoding.html）。１行に、複数の文を書くこともできる。行において、/を2本書くと「//」、その右側はコメントになる。「/*」と「*/」で囲まれた部分もコメントとして扱われ、これは複数行に亘ることができる。
「data {}」にデータを記述する。リスト1.1では、整数値のデータNとkが宣言されている。具体的な値は、stan.buildで与えられている。
Stanではいろいろな型が用意されているが（https://mc-stan.org/docs/reference-manual/types.html）、リスト1.1では整数を表すint型が用いられている。変数名は、Nとkである。データを表す変数（識別子）の宣言の形式は
型 識別子
である。
「parameters {}」に、確率モデルのパラメータとしてpが宣言されている。この変数は、確率なので０以上１以下の値を取る。変数の最小値と最大値を指定した変数宣言の形式
real<lower=0,upper=1>
p;
となっている。変数の型名に
<lower=値1,upper=値２>
の形式で、最小値と最大値を設定する。最小値あるいは最大値の制約が無いときは省く。
パラメータ変数に最小値と最大値を設定すると、その間の一様分布が事前分布として設定されたと見なされる。「modelブロック」において、pの事前分布を
p ~ beta(1,
1);
と設定すると、実行時に
Messages from stanc:
Warning: The parameter p has 2 priors.
と、警告が出される。一様分布は確率密度関数が定数値であるので、事前分布が２つ重なって設定されてもMCMCへの影響はない。
「modelブロック」に確率モデルが設定される。
k ~
binomial(N, p);
２項分布についてはhttps://mc-stan.org/docs/functions-reference/bounded_discrete_distributions.htmlに説明がある。
Stanスクリプトをデータとともに次のコードでビルドする（https://pystan.readthedocs.io/en/latest/reference.html）。
sm = stan.build(StanCode,
data={'N':10, 'k':7})
第１引数にStanスクリプトを置き、引数dataにデータ値を辞書型で設定している。
ビルドしたものの関数sample()を実行してMCMCサンプリングが行われる。
fit = sm.sample()
MCMCサンプルを、次のコードでpandasのDataFrame型に変換している。
df_fit = fit.to_frame()
パラメータpのMCMCサンプルをファイルparam_p.pklに保存して、スクリプトの実行終了である。
with open('param_p.pkl', 'wb') as f:
pickle.dump(df_fit['p'], f)
ファイルparam_p.pklに保存されたパラメータｐのMCMCサンプルを読み出して、KDE (Kernel Density Estimate)曲線を描くスクリプトをリスト1.2のように用意した。
リスト1.2 事後分布の描画
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
with open('param_p.pkl', 'rb') as f:
param_p = pickle.load(f)
sb.kdeplot(param_p)
plt.show()
リスト1.2のスクリプトをWndows10にインストールしたPython 3.13で実行すると図1.1のグラフを得る。
図1.1
正規分布モデルを用いるStanスクリプト例をリスト2.1のように作成してみた。
import numpy as np
import stan
import pickle
x = [54, 73, 47, 61, 62, 53, 58, 63,
65, 58,
68, 69, 49, 63, 62,
60, 63, 64, 44, 61]
with open('x.pkl', 'wb') as f:
pickle.dump(x, f)
StanCode = """
data {
int
n;
array[n] real x;
real
min_x;
real
max_x;
}
parameters {
real
mu;
real<lower=0> sgm;
}
model {
mu ~
normal((min_x+max_x)/2, max_x-min_x); // muの事前分布
sgm
~ exponential(1/(max_x-min_x));
// sgmの事前分布
x ~
normal(mu, sgm);
// 正規分布モデル
}
"""
sm = stan.build(StanCode,
data={'n':len(x), 'x':x,
'min_x':np.min(x), 'max_x':np.max(x)})
fit = sm.sample()
df_fit = fit.to_frame()
mus = df_fit['mu']
# muの事後分布サンプル
sgms = df_fit['sgm']
# sgmの事後分布サンプル
print(mus.describe())
print(sgms.describe())
with open('mus.pkl', 'wb') as f:
pickle.dump(mus, f)
with open('sgms.pkl', 'wb') as f:
pickle.dump(sgms, f)
dataブロックでは、データ数n、データの1次元配列x、最小値min_xと最大値max_xが宣言されている。
配列は、arrayで始まる形式で
array[要素の数] 要素の型 配列名
の形式で宣言される（https://mc-stan.org/docs/stan-users-guide/matrices-arrays.html）。
パラメータsgmは、標準偏差を表すので非負の制約を置いて
real<lower=0> sgm;
と宣言している。
リスト2.1を実行すると、次の警告が提示される。
Warning: The parameter sgm has no
priors. This means either no prior is
provided, or the prior(s)
depend on data variables. In the later case,
this may be a false
positive.
Warning: The parameter mu has no
priors. This means either no prior is
provided, or the prior(s)
depend on data variables. In the later case,
this may be a false
positive.
これは、「false positive」の場合で、事前分布がデータ変数をパラメータとしている（https://mc-stan.org/docs/functions-reference/unbounded_continuous_distributions.html；https://mc-stan.org/docs/functions-reference/positive_continuous_distributions.html#exponential-distribution）。
mu ~ normal((min_x+max_x)/2,
max_x-min_x); // muの事前分布
sgm ~
exponential(1/(max_x-min_x));
// sgmの事前分布
事前分布として一様分布に近いものを選んでいる（Gelman et al., 2021）。
データxの分布が
x ~ normal(mu, sgm);
で与えられているが、これは次の文
for (i in 1:n)
x[i] ~ normal(mu, sgm);
をベクトル化したものである。
MCMCサンプリングにより得られたサンプルは、パラメータmuの事後分布のサンプルmusとパラメータsgmの事後分布のサンプルsgmsが取り出される。
mus = df_fit['mu']
# muの事後分布サンプル
sgms = df_fit['sgm']
# sgmの事後分布サンプル
それらは、それぞれファイルmus.pklとsgms.pklに保存される。
with open('mus.pkl', 'wb') as f:
pickle.dump(mus, f)
with open('sgms.pkl', 'wb') as f:
pickle.dump(sgms, f)
ファイルmus.pklとsgms.pklに保存されたMCMCサンプルを読み込んで分析するスクリプトをリスト2.2のように用意した。
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
import numpy as np
import scipy.stats as ss
with open('x.pkl', 'rb') as f:
x = pickle.load(f)
with open('mus.pkl', 'rb') as f:
mus = pickle.load(f)
with open('sgms.pkl', 'rb') as f:
sgms = pickle.load(f)
sb.kdeplot(mus)
plt.show()
sb.kdeplot(sgms)
plt.show()
#
# データと分析結果のグラフの表示
#
mu = np.median(mus) # 中央値を推定値とする
sgm = np.median(sgms) # 中央値を推定値とする
x_min = np.min(x)
x_max = np.max(x)
x_coord = np.linspace(x_min, x_max,
1000)
y_coord = ss.norm.pdf(x_coord, loc=mu,
scale=sgm) # 正規分布の値
plt.hist(x, density=True)
plt.plot(x_coord, y_coord, lw=3,
label=f'mu={mu:.1f}\nsgm={sgm:.1f}')
plt.xlabel('x', fontsize=18)
plt.ylabel('Density', fontsize=14)
plt.legend()
plt.show()
リスト2.2のスクリプトを、リスト2.1のスクリプトの次に実行する。リスト2.1はMiniconda/Ubuntu/WSL/WindowsのPython 3.12仮想環境で実行したが、リスト2.2のスクリプトはWindows上のPython 3.13で実行した。まずパラメータmuの事後分布のグラフが表示される（図2.1）。
図2.1
図2.1のWindowsを閉じると、パラメータsgmの事後分布のグラフが表示される（図2.2）。
図2.2
図2.2のWindowを閉じると、図2.3のグラフが表示される。
図2.3
データのヒストグラムと事後分布の点推定値に基づく正規分布関数のグラフである。点推定値は中央値を採用している（Gelman et al., 2021）。
図2.3のWindowを閉じると、スクリプトの実行終了である。
2変量正規分布を用いたStanスクリプト例をリスト3.1に示す。
import numpy as np
import stan
import pickle
xy = [[25, 40], [57, 65], [60, 68],
[38, 50], [55, 65],
[45, 49], [36,
55], [42, 54], [54, 61], [44, 59],
[33, 73], [47,
56], [39, 54], [44, 61], [58, 63],
[53, 58], [32,
53], [52, 70], [55, 64], [21, 43]]
xy = np.array(xy)
with open('xy.pkl', 'wb') as f:
pickle.dump(xy, f)
StanCode = """
data {
int
n;
array[n] vector[2] xy;
real
xmin;
real
xmax;
real
ymin;
real
ymax;
}
parameters {
/*
再スケール前のサンプル値
*/
real<lower=0, upper=1> vrho;
real<lower=0, upper=1> vxmu;
real<lower=0, upper=1> vymu;
real<lower=0, upper=1> vxsgm;
real<lower=0, upper=1> vysgm;
}
transformed parameters {
/*
データに合わせてパラメータのスケーリングを行う
*/
real
rho;
real
xmu;
real
ymu;
real
xsgm;
real
ysgm;
rho
= -0.9999 + vrho*1.9998;
xmu
= xmin + (xmax - xmin) * vxmu;
ymu
= ymin + (ymax - ymin) * vymu;
xsgm
= 0.0001 + (xmax - xmin) * vxsgm;
ysgm
= 0.0001 + (ymax - ymin) * vysgm;
}
model {
/*
スケーリング前のパラメータの事前分布
*/
vrho
~ uniform(0,1);
vxmu
~ uniform(0,1);
vymu
~ uniform(0,1);
vxsgm ~ uniform(0,1);
vysgm ~ uniform(0,1);
/*
２変量正規分布
*/
xy ~
multi_normal([xmu, ymu],
[[xsgm^2, xsgm*ysgm*rho],
[ysgm*xsgm*rho,
ysgm^2]]);
}
"""
sm = stan.build(StanCode,
data={'n':len(xy), 'xy':xy,
'xmin':np.min(xy.T[0]), 'xmax':np.max(xy.T[1]),
'ymin':np.min(xy.T[1]), 'ymax':np.max(xy.T[1])})
fit = sm.sample()
df_fit = fit.to_frame() # MCMCサンプルのDataFrame型への変換
xmus = df_fit['xmu']
xsgms = df_fit['xsgm']
ymus = df_fit['ymu']
ysgms = df_fit['ysgm']
rhos = df_fit['rho']
print(xmus.describe())
print(xsgms.describe())
print(rhos.describe())
"""
MCMCサンプルの保存
"""
with open('rhos.pkl', 'wb') as f:
pickle.dump(rhos, f)
with open('xmus.pkl', 'wb') as f:
pickle.dump(xmus, f)
with open('xsgms.pkl', 'wb') as f:
pickle.dump(xsgms, f)
with open('ymus.pkl', 'wb') as f:
pickle.dump(ymus, f)
with open('ysgms.pkl', 'wb') as f:
pickle.dump(ysgms, f)
dataブロックにおいて、２変量のデータを
array[n] vector[2] xy;
と宣言している。
２変量正規分布（https://mc-stan.org/docs/functions-reference/distributions_over_unbounded_vectors.html）のパラメータは、［0, 1］上の一様分布を事前分布とするものを再スケールしている。これは、事前分布において変域を大きくとる
vxmu ~
uniform(0,100);
と、以下のような警告が出るからである。
Warning in
'/tmp/httpstan_nexxzom0/model_vtzo6piu.stan', line 31, column 25: Argument
100 suggests there may be
parameters that are not unit scale; consider
rescaling with a multiplier
(see manual section 22.12).
再スケーリングにおいて、rhoの範囲をー１と＋１を含まないように
rho = -0.9999 + vrho*1.9998;
としているのは、分散共分散行列が非特異行列になるのを避けるためである。
リスト3.1のスクリプトは、MCMCサンプルをファイルに保存して、実行が終了する。
リスト3.1のスクリプトにより保存されたMCMCサンプルは、Windows上で実行されるリスト3.2のスクリプトにより読み込んで分析される。
import numpy as np
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
import scipy.stats as ss
"""
MCMCサンプルの読み込み
"""
with open('xy.pkl', 'rb') as f:
xy = pickle.load(f)
with open('rhos.pkl', 'rb') as f:
rhos = pickle.load(f)
with open('xmus.pkl', 'rb') as f:
xmus = pickle.load(f)
with open('xsgms.pkl', 'rb') as f:
xsgms = pickle.load(f)
with open('ymus.pkl', 'rb') as f:
ymus = pickle.load(f)
with open('ysgms.pkl', 'rb') as f:
ysgms = pickle.load(f)
"""
事後分布の描画
"""
sb.kdeplot(rhos)
plt.show()
sb.kdeplot(xmus)
plt.show()
sb.kdeplot(ymus)
plt.show()
sb.kdeplot(xsgms)
plt.show()
sb.kdeplot(ysgms)
plt.show()
"""
データのプロットと
推定された２変量正規分布のグラフ
"""
rho = np.median(rhos)
xmu = np.median(xmus)
xsgm = np.median(xsgms)
ymu = np.median(ymus)
ysgm = np.median(ysgms)
mins = np.min(xy, axis=0) -
np.array([2, 2])
maxs = np.max(xy, axis=0) +
np.array([2, 2])
x_range = np.linspace(mins[0],
maxs[0], 100)
y_range = np.linspace(mins[1],
maxs[1], 100)
X, Y = np.meshgrid(x_range, y_range)
xypos = np.dstack((X, Y))
Z = ss.multivariate_normal.pdf(xypos,
[xmu, ymu],
[[xsgm**2,
xsgm*ysgm*rho],
[ysgm*xsgm*rho, ysgm**2]])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xy.T[0], xy.T[1], 'o') # データの散布図
ax.contour(x_range, y_range, Z) # 正規分布の等高線図
plt.title(fr'$\rho$={rho:.3f}',
fontsize=20)
plt.show()
リスト3.2のスクリプトを実行すると、まずパラメータrhoの事後分布のグラフが表示される（図3.1）。
図3.1
図3.1のWindowを閉じると、パラメータxmuの事後分布のグラフが表示される（図3.2）。
図3.2
図3.2のWindowを閉じると、パラメータymuの事後分布のグラフが表示される（図3.3）。
図3.3
図3.3のWindowを閉じると、xsgmの事後分布のグラフが表示される（図3.4）。
図3.4
図3.4のWindowを閉じると、ysgmの事後分布のグラフが表示っされる（図3.5）。
図3.5
図３.5のWindowを閉じると、図3.6のグラフが表示される。
図3.6
正規分布の等高線は、各パラメータの事後分布の中央値を点推定値として描かれている（Gelman et al., 2021）。
等高線は、次のコード
ax.contour(x_range, y_range, Z)
により描いている（https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.contour.html）。
図3.6のWindowを閉じると、スクリプトの実行終了である。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press. https://sites.stat.columbia.edu/gelman/book/BDA3.pdf
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
岡本安晴（2019）「いまさら聞けないPythonでデータ分析」丸善出版