PyStan 3 試用記
2022.10, 2023.03
PyStan 3を使ってみた（2022.10.31, 2023.03.18）。Anacondaをインストールして、Python 3.10の環境を用意し、PyStan 3をコマンド
pip install
pystan
でインストールした。Pythonのスクリプト内でPyStan 3をimportするときは
import stan
である。
以下に、２項分布、単回帰モデル、2変量正規分布の例を示す。スクリプトファイルは、圧縮ファイルsample_script.zipにまとめたので、ダウンロード・展開して自由に利用されたい。
PyStan 3の使用例は、このウェブサイトも参照されたい。
PyStan 3は直接Windows上で実行できない。Windowsの場合は、WSL (Windows Subsystem for Linux) をWindowsにインストールして、WSL上でPyStan 3を利用する。WSLについては、このサイトで説明している。ここでの説明とは現在は少し異なった手順のようであるが、参考になると思う。また、WSLのインストールがうまくいかなくても、日にちをおいて試みれば、うまくいくことがある。WSLの使えるWindowsの範囲が広がっているようである。
PyStan 3をインストールするとき、あるいはPyStan 3実行するとき、インストールできない、あるいは実行できないときは、C++コンパイラがインストールされていない可能性がある。コマンド「g++」を実行してみて、「g++」が見つかりません、コマンド「sudo apt install g++」でインストールして下さいというようなメッセージが表示されたら、その通り実行してみる。
「g++」のインストールが済めば、PyStan 3のインストールができて、PyStan 3スクリプトの実行ができるはずである。
まず、Ubuntuでの利用について説明する。Windowsの場合については、このウェブページの後半で説明する。
PyStanについては、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。
Ubuntuの場合
試用したスクリプトは以下のようである。
リスト１ 2項分布の例
import stan
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
binomial_model = """
data {
int
N;
int
k;
}
parameters {
real<lower = 0.0, upper = 1.0> p;
}
model {
k ~
binomial(N, p);
}
"""
binomial_data = {'N':10, 'k':7}
sm = stan.build(binomial_model, data =
binomial_data)
fit = sm.sample()
print('fit =\n', fit, '\n')
p_smpls = fit['p']
print('p_smpls =\n', p_smpls)
sb.kdeplot(p_smpls[0])
plt.savefig('FigKDE.png')
plt.show()
plt.close()
df = fit.to_frame()
sb.kdeplot(df['p'])
plt.savefig('FigPD.png')
plt.show()
plt.close()
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
az.plot_trace(inference_data)
plt.savefig('FigTrace.png')
plt.show()
plt.close()
リスト１のスクリプトを実行すると、以下の端末画面への出力と描画グラフの保存ファイルが得られる。
fit =
<stan.Fit>
Parameters:
p: ()
Draws: 4000
p_smpls =
[[0.62546345 0.68693505 0.64059199 ...
0.76143154 0.69859012 0.78014042]]
FigKDE.png
FigPD.png
FigTrace.png
リスト２のスクリプトを用意した。
リスト２ 単回帰モデルの例
import stan
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
Xs = ss.uniform.rvs(loc = 10, scale =
20, size = 100)
Es = ss.norm.rvs(loc = 0, scale = 1,
size = 100)
Ys = 10 + 0.5 * Xs + Es
plt.scatter(Xs, Ys)
plt.savefig('FigXY.png')
plt.show()
plt.close()
simple_reg_model = """
data {
int
N;
vector[N] X;
vector[N] Y;
}
parameters {
real
a;
real
b;
real<lower = 0> sgm;
}
model {
a ~
normal(0, 10);
b ~
normal(0, 10);
sgm
~ exponential(0.1);
Y ~ normal(a + b*X,
sgm);
}
"""
simple_reg_data = {'N':len(Xs),
'X':Xs, 'Y':Ys}
sm = stan.build(simple_reg_model, data
= simple_reg_data)
fit = sm.sample()
print('fit =\n', fit, '\n')
infer_data = az.from_pystan(posterior=fit,
posterior_model=sm)
smry = az.summary(infer_data)
print(smry)
az.plot_trace(infer_data)
plt.savefig('FigTrace.png')
print('FigTrace.png was saved.')
plt.show()
plt.close()
a_smpls = fit['a']
b_smpls = fit['b']
plt.scatter(a_smpls, b_smpls)
plt.xlabel('a')
plt.ylabel('b')
plt.savefig('Fig_ab.png')
plt.show()
plt.close()
リスト２のスクリプトを実行すると、端末画面上に以下の出力を得るとともに、グラフがファイルに保存される。
fit =
<stan.Fit>
Parameters:
a: ()
b: ()
sgm: ()
Draws: 4000
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
a 9.913 0.362 9.211 10.556 0.010 0.007 1272.0 1507.0 1.0
b 0.502 0.018 0.469 0.535 0.000 0.000 1273.0 1509.0 1.0
sgm 0.996 0.074 0.860 1.128 0.002 0.001 1856.0 1865.0 1.0
FigXY.png
FigTrace.png
Fig_ab.png
2変量正規分布を用いたスクリプト例をリスト３のように用意した。
リスト３ 2変量正規分の例
import stan
import scipy.stats as ss
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
X = ss.multivariate_normal.rvs([0,0],[[1,
0.8],[0.8,1]], size = 100)
print(X)
plt.scatter(X.T[0], X.T[1])
plt.savefig('FigXY.png')
plt.show()
plt.close()
binormal_model = """
data {
int
N;
vector[2] X[N];
}
parameters {
vector[2]
mu;
vector[2] sgm;
real<lower = -1, upper = 1> rho;
}
transformed parameters {
cov_matrix[2] cov;
cov[1][1] = sgm[1]*sgm[1];
cov[1][2] = rho*sgm[1]*sgm[2];
cov[2][1] = rho*sgm[1]*sgm[2];
cov[2][2] = sgm[2]*sgm[2];
}
model {
mu[1] ~ normal(0, 10);
mu[2] ~ normal(0, 10);
sgm
~ exponential(0.1);
rho
~ uniform(-1,1);
X ~
multi_normal(mu, cov);
}
"""
binormal_data = {'N':len(X), 'X':X}
sm = stan.build(binormal_model, data =
binormal_data)
fit = sm.sample()
print('fit =\n', fit, '\n')
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
az.plot_trace(inference_data)
plt.savefig('FigTrace.png')
print('FigTrace.png was saved.')
plt.show()
plt.close()
mu_smpls = fit['mu']
print('mu_smpls =\n', mu_smpls)
df = fit.to_frame()
print('df =\n', df, '\n\n')
for v in zip(df['mu.1'][:5], df['mu.2'][:5],
df['cov.2.2'][:5]):
print(v)
sb.kdeplot(df['rho'])
plt.savefig('FigRho.png')
plt.show()
plt.close()
リスト３のスクリプトを実行すると、端末画面上に以下の出力を得るとともに、グラフがファイルに保存される。
fit =
<stan.Fit>
Parameters:
mu: (2,)
sgm: (2,)
rho: ()
cov: (2, 2)
Draws: 4000
Summary =
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu[0] 0.187 0.105 -0.007 0.383 0.002 0.001 2472.0 2695.0 1.01
mu[1] 0.210 0.097 0.024 0.388 0.002 0.001 2498.0 2694.0 1.00
sgm[0] 1.065 0.076 0.928 1.210 0.002 0.001 2180.0 2205.0 1.00
sgm[1] 0.984 0.069 0.853 1.112 0.001 0.001 2366.0 2355.0 1.00
rho 0.770 0.040 0.691 0.839 0.001 0.001 2295.0 2620.0 1.00
cov[0,0] 1.140 0.164 0.862 1.463 0.004 0.003 2180.0 2205.0 1.00
cov[0,1] 0.813 0.132 0.582 1.062 0.003 0.002 1880.0 2135.0 1.00
cov[1,0] 0.813 0.132 0.582 1.062 0.003 0.002 1880.0 2135.0 1.00
cov[1,1] 0.974 0.138 0.721 1.229 0.003 0.002 2366.0 2355.0 1.00
FigTrace.png was saved.
mu_smpls =
[[0.18040335 0.16903765 0.14978146 ...
0.43856426 0.01501792 0.24142373]
[0.20285295 0.15172237 0.22550795 ...
0.43964007 0.14393881 0.22252268]]
df =
parameters lp__ accept_stat__ stepsize__ treedepth__ ... cov.1.1 cov.2.1 cov.1.2 cov.2.2
draws
...
0
-58.263544
0.827678
0.458740
3.0 ... 1.133158 0.922741 0.922741 1.122094
1
-64.044667
0.562658
0.457346
2.0 ... 1.353705 0.763930 0.763930 0.725858
2
-57.713119
0.988581
0.430489
2.0 ...
1.037119 0.775723 0.775723 0.910998
3
-58.287036
0.957960
0.466302
3.0 ... 1.052202 0.752236 0.752236 0.959562
4
-57.652823
0.958025
0.458740
3.0 ... 1.041470 0.715984 0.715984 0.926422
...
...
...
...
... ... ... ... ... ...
3995
-59.680630
0.975303
0.466302
3.0 ... 1.104269 0.755304 0.755304 1.030205
3996
-57.830879 0.998800 0.458740
4.0 ... 0.920630 0.625547 0.625547 0.809076
3997
-60.108280
0.971335
0.457346
3.0 ... 1.202069 0.903426 0.903426 1.077758
3998
-58.572426
0.983354
0.430489
3.0 ... 1.125427 0.783323 0.783323 0.908921
3999
-58.904657
0.998236
0.466302
3.0 ... 1.055764 0.672293 0.672293 0.739689
[4000 rows x 16 columns]
(0.18040335026726861,
0.20285294781077376, 1.122094106032836)
(0.16903765146625016, 0.15172236853563245,
0.7258576891080232)
(0.1497814603342748,
0.22550794533864094, 0.9109977685363938)
(0.22654457223130214,
0.3351371542279542, 0.959561947061149)
(0.2556725357232705,
0.2796258006706391, 0.9264216126600003)
FigXY.png
FigTrace.png
FigRho.png
Windowsの場合のPyStan 3の利用
PyStan 3は直接Windows上で実行できないので、WSL (Windows Subsystem for Linux)をインストールして、WSL上のUbuntuでPyStan 3を利用する。WSLについては、このウェブサイトを参照されたい。
Ubuntu/WSLにおいてAnacondaをインストールしてPython 3.10の仮想環境を用意する。Python 3.10仮想環境において、コマンド
pip install
pystan
により、PyStan 3をインストールする。
Pythonスクリプトにおいて、PyStan 3をimportするときは
import stan
である。
ベイズ分析は、PyStan 3によるMCMCサンプリング（ステップ１）と、MCMCサンプリングの結果に基づく分析（ステップ２）に分かれる。
ステップ１は、PyStan 3の実行のために、仮想環境Python3.10/Anaconda/WSL/Windowsで行い、ＭＣＭＣサンプリングの結果をファイルに保存する。ステップ２において、保存されたＭＣＭＣサンプリングのデータを読み込んで。描画などの分析作業をWindows上で行う。ただし、現在（2023.08）では、仮想環境Python3.10/Anaconda/WSL/Windowsでもステップ２の実行が可能である。
PyStan ３によるMCMCサンプリングを仮想環境Python 3.10/Anaconda/Ubuntu/WSL/Windows 10で行った。モデルは、2変量正規分布である。スクリプトをリストW1に示す。
なお、最近（2023.04.13）インストールしたUbuntu/WSLでは、stan.buildの実行時にcompileエラーがでた。これは、UbuntuにC++コンパイラーがインストールされていないことに依るようである。
コマンド
g++
を実行すると、g++がインストールされていないというメッセージが表示された。g++をインストールするためにコマンド
sudo apt
install g++
を実行してインストールした。その後、Python 3.10の仮想環境でstan.buildを実行すると、うまくいった。
リストW1 PyStan 3によるMCMCサンプリング（２変量正規分布）
import stan
import scipy.stats as ss
import arviz as az
import pickle
X =
ss.multivariate_normal.rvs([0,0],[[1, 0.8],[0.8,1]], size = 100)
print(X)
binormal_model = """
data {
int
N;
vector[2] X[N];
}
parameters {
vector[2] mu;
vector[2] sgm;
real<lower = -1, upper = 1> rho;
}
transformed parameters {
cov_matrix[2] cov;
cov[1][1] = sgm[1]*sgm[1];
cov[1][2] = rho*sgm[1]*sgm[2];
cov[2][1] = rho*sgm[1]*sgm[2];
cov[2][2] = sgm[2]*sgm[2];
}
model {
mu[1] ~ normal(0, 10);
mu[2] ~ normal(0, 10);
sgm
~ exponential(0.1);
rho
~ uniform(-1,1);
X ~
multi_normal(mu, cov);
}
"""
binormal_data = {'N':len(X), 'X':X}
sm = stan.build(binormal_model, data =
binormal_data)
fit = sm.sample()
inference_data =
az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(inference_data)
print('Summary =\n',smry)
with open('i_data.pkl', 'wb') as f:
pickle.dump(inference_data,
f)
print('i_data.pkl was saved.')
df = fit.to_frame()
with open('d_frame.pkl', 'wb') as f:
pickle.dump(df, f)
print('d_frame.pkl was saved.')
print('df =\n', df, '\n\n')
for v in zip(df['mu.1'][:5],
df['mu.2'][:5], df['cov.2.2'][:5]):
print(v)
print('End')
リストＷ１のスクリプトを実行すると、PyStan 3によるＭＣＭＣサンプリングが行われ、MCMCサンプリングの結果を、InferenceData型にまとめたものがファイルi_data.pklに、pandas DataFrame型にまとめたものがファイルd_frame.pklに格納される。
次に、リストＷ２のスクリプトにより、これらのファイルからＭＣＭＣサンプリングのデータが読み込まれて描画などの処理が行われる。
ステップ１で作成されたＭＣＭＣサンプリングのデータのファイルを、Python 3.11/Windowsで読み込み、描画などを行った。スクリプトをリストＷ２に示す。
リストＷ２ ＭＣＭＣサンプリングデータの描画など
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
import pickle
with open('i_data.pkl', 'rb') as f:
inference_data =
pickle.load(f)
smry = az.summary(inference_data)
print('Summary =\n',smry)
az.plot_trace(inference_data)
plt.tight_layout()
plt.show()
plt.close()
with open('d_frame.pkl', 'rb') as f:
df = pickle.load(f)
print('df =\n', df, '\n\n')
for v in zip(df['mu.1'][:5],
df['mu.2'][:5], df['cov.2.2'][:5]):
print(v)
sb.kdeplot(df['mu.1'], label = 'mu.1')
sb.kdeplot(df['mu.2'], label = 'mu.2')
plt.xlabel('mu')
plt.legend()
plt.show()
print('End')
リストW2のスクリプトを実行すると、まず図W1のグラフが表示される。
図W1
図W1のグラフを閉じると、次に図W2のグラフが表示される。
図W2
端末には、以下の出力などが表示されている。
Summary =
mean
sd hdi_3% hdi_97% ...
mcse_sd ess_bulk ess_tail r_hat
mu[0] 0.111 0.107 -0.093 0.303 ... 0.002 2183.0 2089.0 1.0
mu[1] -0.009 0.101 -0.202 0.176 ... 0.002 2232.0 2328.0 1.0
sgm[0] 1.046 0.074 0.912 1.187 ... 0.001 1844.0 2132.0 1.0
sgm[1] 0.989 0.070 0.861 1.122 ... 0.001 2008.0 2336.0 1.0
rho
0.833 0.031 0.771 0.885 ... 0.000 2251.0 2075.0 1.0
cov[0, 0] 1.100 0.158 0.828 1.404 ... 0.003 1844.0 2132.0 1.0
cov[0, 1] 0.867 0.137 0.618 1.124 ... 0.002 1713.0 1974.0 1.0
cov[1, 0] 0.867 0.137 0.618 1.124 ... 0.002 1713.0 1974.0 1.0
cov[1, 1] 0.983 0.141 0.729 1.247 ... 0.002 2008.0 2336.0 1.0