PyStan3のbuild関数とsample関数の備忘録
PyStan3を用いてのMCMCサンプリングはbuild関数とsample関数を用いるとできる。簡単な例をリスト１に示す。
なお、PyStan3用Stanの備忘録も用意した。
PyStan3のインストールは
pip install
pystan
であるが、importは
import stan
である。
Pythonスクリプトの実行は、Anaconda/ubuntu22/WSL/Windowsに設定したPython3.12仮想環境において行った。ただし、scipy1.12.0を利用している。このウェブサイトを参照されたい。
import stan
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);
}
"""
sm = stan.build(binomial_model,
data={'n':10, 'k':7})
fit = sm.sample()
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
文字列binomial_modelとして記述されたStanモデルと、辞書型で与えられたデータdata={'n':10, 'k':7}を引数として関数stan.buildを呼び出してビルドを行い、作成したモデルsmの関数sm.sampleを呼び出して、MCMCサンプリングを行っている。MCMCサンプリングの結果は、fitによって表されている。
fitの内容は、arvizの関数from_pystanによってInferenceData型のオブジェクトに変換した後、arvizの関数summaryによってMCMCサンプリングの統計量が算出されるが、これらを次のスクリプトで行っている。
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
端末には、以下のように出力されている。
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.669 0.132 0.426 0.904 0.004 0.003 1293.0 1444.0 1.0
関数buildとsampleにリスト１における例以外のパラメータも設定した例は、リスト４に示している。
Stanモデルが別のファイルに記述されている場合
リスト１では、StanモデルがPythonスクリプト内に記述されているが、これを別ファイルとして用意することもできる。
リスト２は、Stanモデルを別ファイルbinomial_pstan.stanとして用意したものである。
リスト２ Stanモデル（binomial_pstan.stan）
data {
int n;
int k;
}
parameters {
real<lower=0.0,
upper=1.0> p;
}
model {
k ~ binomial(n, p);
}
リスト２のファイルを読み込んで利用するPythonスクリプトを、リスト３に示す。
リスト３ リスト２のファイルbinomial_pstan.stanを読み込んで利用するPythonスクリプト
import stan
import arviz as az
with open('binomial_pstan.stan', 'r')
as f_stan:
sm =
stan.build(f_stan.read(), data={'n':10, 'k':7})
fit = sm.sample()
print(fit)
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
Stanモデルのファイルbinomial_pstan.stanをf_stan.read()で読み込んで文字列として、stan.buildの第1パラメータとしている。
with open('binomial_pstan.stan', 'r')
as f_stan:
sm =
stan.build(f_stan.read(), data={'n':10, 'k':7})
Stanモデルを与える文字列が、ファイルの読み込みによって与えられていることを除いて、リスト１と同じである。
関数buildとsampleのよく用いる可能性のあるパラメータの使用例を、リスト４に示す。
import stan
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);
}
"""
sm = stan.build(binomial_model,
data={'n':10, 'k':7}, random_seed=12345)
fit = sm.sample(num_chains=4,
init=[{'p':0.4}, {'p':0.5},{'p':0.6}, {'p':0.7}],
num_samples
= 2000,
num_warmup = 1000)
print(fit)
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
関数buildにおけるパラメータ
random_seed=12345
は、乱数の初期値（正整数）を設定するものである。スクリプトの実行を繰り返すとき、同じMCMCサンプルを得たいときに設定する。
関数sampleにおけるパラメータ
num_chains=4,
は、MCMCサンプリングにおけるマルコフ鎖の本数を指定するものである。
パラメータ
init=[{'p':0.4}, {'p':0.5},{'p':0.6}, {'p':0.7}]
は、マルコフ鎖の初期値を指定するものである。マルコフ鎖の本数だけ設定する。マルコフ鎖がwarm up時に安定しないとき、適切な初期値を設定すると安定することがある。
パラメータ
num_samples = 2000,
num_warmup = 1000
は、MCMCサンプル数とwarm up試行数を設定するものである。この指定は、各マルコフ鎖に対するものであり、warm up試行数のサンプリング後、num_samples数のサンプリングが行われる。MCMCサンプリングが安定しないときはnum_warmup数を増やして安定を期待することができる。事後分布の推定精度を上げるときは、num_samplesの値を大きくとることになる。
リスト４のスクリプトを実行すると、端末には以下の出力が表示される。
<stan.Fit>
Parameters:
p: ()
Draws: 8000
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.67 0.134 0.423 0.907 0.002 0.002 3140.0 3622.0 1.0