PyStan3用Stanメモ
PyStan3用Stanにおける事前分布（一様分布）、事前分布（分布のスケーリング）、配列について備忘録を用意した。
PyStan3の関数buildとsampleについての備忘録ウェブサイトも用意している。
パラメータの変域に制約を設定すると、一様分布として扱われるようである。一様分布は、Gelmanら（2021）では、分析の第１ステップで用いられる事前分布と説明されている（p.124）。
リスト1.1に示すスクリプトは、パラメータpの変域のみを設定するものである。
リスト1.1 事前分布を明示しない例
import stan
import arviz as az
model_str = """
data {
int
n;
int
k;
}
parameters {
real<lower=0.0, upper=1.0> p;
}
model {
k ~
binomial(n, p);
}
"""
sm = stan.build(model_str,
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))
リスト1.1を実行すると、下記の結果を得る。
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$ python uni_0.py
Building: found in cache, done.
Sampling: 100% (8000/8000), done.
Messages received during sampling:
Gradient evaluation took 1.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.18
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 1.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.18
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 1.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.15
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 1.7e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.17
seconds.
Adjust your expectations accordingly!
<stan.Fit>
Parameters:
p: ()
Draws: 4000
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.667 0.13 0.422 0.889 0.003 0.002 1416.0 1534.0 1.0
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$
リスト1.2は、ｐの事前分布として一様分布を明示的に設定したものである。
リスト1.2 事前分布の明示的設定
import stan
import arviz as az
model_str = """
data {
int
n;
int
k;
}
parameters {
real<lower=0.0, upper=1.0> p;
}
model {
p ~
uniform(0, 1);
// 事前分布の設定
k ~
binomial(n, p);
}
"""
sm = stan.build(model_str,
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))
リスト1.2を実行すると、以下の結果を得る。
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/www/****/samplefiles$ python uni_1.py
Building: found in cache, done.
Messages from stanc:
Warning: The parameter p has 2 priors.
Sampling: 100% (8000/8000), done.
Messages received during sampling:
Gradient evaluation took 3.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.32
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2.3e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.23
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 4.7e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.47
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2.2e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.22
seconds.
Adjust your expectations accordingly!
<stan.Fit>
Parameters:
p: ()
Draws: 4000
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
p
0.666 0.13 0.434 0.91 0.004 0.003 1220.0 1543.0 1.01
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$
上の出力の４行目に次の警告メッセージが出力されている。
Warning: The parameter p has 2 priors.
モデル宣言部での
p ~ uniform(0, 1);
による事前分布の設定に加えて、パラメータ宣言部の
real<lower=0.0, upper=1.0> p;
により、一様分布が暗黙に設定されているようである。
この暗黙の設定は、常に行われているのではないようである。リスト3.1の場合は、パラメータrhoに２つの事前分布が設定されているという警告は出ない。
ちなみに、一様分布の確率密度関数値は定数であるので、微分すると０になり消える。したがって、一様分布は、重複して設定しても対数事後分布（log-posterior density; Gelmanら, 2014, p.301）の微分には影響しないが、モデルの適合度を算出するために対数尤度を求めるときは注意が必要である。
引用文献
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University Press.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, third edition. CRC Press.
パラメータの事前分布の変域が広いと警告が出される。
リスト2.1では、パラメータmuとsgmの変域を十分に広くとっている。
リスト2.1 スケーリングが考慮されていない場合
import stan
import scipy.stats as ss
import arviz as az
model_str = """
data {
int
n;
array[n] real x;
}
parameters {
real
mu;
real<lower=0.0> sgm;
}
model {
mu ~
normal(0.0, 100.0);
sgm ~
exponential(0.01);
x ~
normal(mu, sgm);
}
"""
n = 100
x = ss.norm.rvs(loc=50, scale=15,
size=n)
print(x[:10])
sm = stan.build(model_str,
data={'n':len(x), 'x':x})
fit = sm.sample()
print(fit)
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
リスト2.1を実行すると、以下の出力を得る。
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$ python uni_normal.py
[31.78297055 31.62761252 61.44566053
62.38655417 54.97553818 21.24494148
66.38688251 74.50928458 59.78637668 31.0916614
]
Building: found in cache, done.
Messages from stanc:
Warning in
'/tmp/httpstan_qzlrip31/model_vkggkwq5.stan', line 12, column 26: Argument
0.01 suggests there may be
parameters that are not unit scale; consider
rescaling with a multiplier
(see manual section 22.12).
Warning in
'/tmp/httpstan_qzlrip31/model_vkggkwq5.stan', line 11, column 25: Argument
100.0 suggests there may be
parameters that are not unit scale; consider
rescaling with a multiplier
(see manual section 22.12).
Sampling: 100% (8000/8000), done.
Messages received during sampling:
Gradient evaluation took 3.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.38
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 3.6e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.36
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 3.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.35
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 3.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.35
seconds.
Adjust your expectations accordingly!
<stan.Fit>
Parameters:
mu: ()
sgm: ()
Draws: 4000
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
mu 48.077 1.451 45.403 50.870 0.023 0.016 3907.0 2363.0 1.0
sgm 14.330 1.000 12.564 16.277 0.016 0.012 3799.0 2615.0 1.0
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$
警告メッセージ
Warning in
'/tmp/httpstan_qzlrip31/model_vkggkwq5.stan', line 12, column 26: Argument
0.01 suggests there may be
parameters that are not unit scale; consider
rescaling with a multiplier
(see manual section 22.12).
Warning in
'/tmp/httpstan_qzlrip31/model_vkggkwq5.stan', line 11, column 25: Argument
100.0 suggests there may be
parameters that are not unit scale; consider
rescaling with a multiplier
(see manual section 22.12).
があるので、リスト2.1をリスト2.2のように改めてみる。
import stan
import scipy.stats as ss
import arviz as az
model_str = """
data {
int
n;
array[n] real x;
}
parameters {
real
v_mu;
real<lower=0.0> v_sgm;
}
transformed parameters {
real
mu;
real
sgm;
mu =
v_mu * 100;
// スケーリング
sgm
= 0.001 + v_sgm * 100;
// スケーリング
}
model {
v_mu
~ normal(0.0, 10.0);
// 変域を抑えている
v_sgm ~ exponential(0.1); // 変域を抑えている
x ~
normal(mu, sgm);
}
"""
n = 100
x = ss.norm.rvs(loc=50, scale=15,
size=n)
print(x[:10])
sm = stan.build(model_str,
data={'n':len(x), 'x':x})
fit = sm.sample()
print(fit)
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
リスト2.2を実行すると、以下の結果を得る。
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$ python
uni_normal_scaled.py
[75.264407 36.51387646 49.78169614
41.90098387 43.11311503 68.19862408
71.14342016 51.04580772 55.36301218 74.44421151]
Building: found in cache, done.
Sampling: 100% (8000/8000), done.
Messages received during sampling:
Gradient evaluation took 2.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.25
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.25
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2.4e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.24
seconds.
Adjust your expectations accordingly!
Gradient evaluation took 2.5e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.25
seconds.
Adjust your expectations accordingly!
<stan.Fit>
Parameters:
v_mu: ()
v_sgm: ()
mu: ()
sgm: ()
Draws: 4000
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
v_mu 0.509 0.013 0.485 0.532 0.000 0.000 3410.0 2421.0 1.0
v_sgm 0.127 0.009 0.110 0.143 0.000 0.000 2987.0 1874.0 1.0
mu 50.880 1.272 48.470 53.237 0.022 0.015 3410.0 2421.0 1.0
sgm 12.690 0.897 11.004 14.325 0.017 0.012 2987.0 1874.0 1.0
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$
警告メッセージの出力はない。
配列は、予約語arrayで始める。例を、リスト3.1に示す。
import stan
import scipy.stats as ss
import arviz as az
model_str = """
data {
int
n;
array[n]
vector[2] x;
}
parameters {
vector[2] v_mu;
vector<lower=0.0>[2] v_sgm;
real<lower=-1.0, upper=1.0> rho;
}
transformed parameters {
vector[2] mu;
vector[2] sgm;
mu =
v_mu * 100;
// スケーリング
sgm
= 0.001 + v_sgm * 10;
// スケーリング
}
model {
v_mu
~ normal(0.0, 10.0);
// 変域を抑えている
v_sgm ~ exponential(0.1); // 変域を抑えている
rho
~ uniform(-1, 1);
x ~
multi_normal(mu, [[sgm[1]^2, rho*sgm[1]*sgm[2]],
[rho*sgm[2]*sgm[1], sgm[2]^2]]);
}
"""
n = 100
x =
ss.multivariate_normal.rvs(mean=[50, 55],
cov=[[15**2, 0.5*(15**2)],
[0.5*(15**2), 15**2]],
size=n)
print(x[:10])
sm = stan.build(model_str,
data={'n':len(x), 'x':x})
fit = sm.sample()
print(fit)
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
リスト3.1のdata宣言部において、配列が
array[n] vector[2] x;
と、vector[2]型のn個の配列array[n]として宣言されている。
リスト3.1を実行すると、以下の結果を得る。
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/****/samplefiles$ python
multi_normal_scaled.py
[[41.21371655 68.42238388]
[55.49835866 66.2519432 ]
[41.99030993 57.13919303]
[43.02091292 82.15598872]
[38.46163839 44.91750327]
[63.84962356 89.32400035]
[53.83770669 60.23709107]
[29.83379027 9.78398135]
[55.15986798 41.35552297]
[36.57292359 32.33198216]]
Building: found in cache, done.
Sampling: 100% (8000/8000), done.
Messages received during sampling:
Gradient evaluation took 7.8e-05 seconds
1000 transitions using 10 leapfrog steps per transition would take 0.78
seconds.
Adjust your expectations accordingly!
Informational Message: The current Metropolis proposal is about to be
rejected because of the following issue:
Exception: multi_normal_lpdf: Covariance matrix is not symmetric.
Covariance matrix[1,2] = inf, but Covariance matrix[2,1] = inf (in
'/tmp/httpstan_y8lqar38/model_npbz37m7.stan', line 21, column 8 to line 22,
column 64)
If this warning occurs sporadically, such
as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely
ill-conditioned or misspecified.
Informational Message: The current Metropolis proposal is about to be
rejected because of the following issue:
Exception: multi_normal_lpdf: Covariance matrix is not symmetric.
Covariance matrix[1,2] = inf, but Covariance matrix[2,1] = inf (in
'/tmp/httpstan_y8lqar38/model_npbz37m7.stan', line 21, column 8 to line 22,
column 64)
If this warning occurs sporadically, such
as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely
ill-conditioned or misspecified.
Informational Message: The current Metropolis proposal is about to be
rejected because of the following issue:
Exception: multi_normal_lpdf: Covariance matrix is not symmetric.
Covariance matrix[1,2] = inf, but Covariance matrix[2,1] = inf (in
'/tmp/httpstan_y8lqar38/model_npbz37m7.stan', line 21, column 8 to line 22,
column 64)
・
・
・
If this warning occurs sporadically, such
as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely
ill-conditioned or misspecified.
<stan.Fit>
Parameters:
v_mu: (2,)
v_sgm: (2,)
rho: ()
mu: (2,)
sgm: (2,)
Draws: 4000
mean
sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
v_mu[0] 0.475 0.014 0.449 0.501 0.000 0.000 3494.0 2955.0 1.0
v_mu[1] 0.527 0.016 0.498 0.558 0.000 0.000 3775.0 2789.0 1.0
v_sgm[0] 1.387 0.102 1.197 1.574 0.002 0.001 4026.0 3247.0 1.0
v_sgm[1] 1.640 0.118 1.432 1.863 0.002 0.001 4101.0 2936.0 1.0
rho
0.462 0.082 0.305 0.602 0.001 0.001 3683.0 2837.0 1.0
mu[0] 47.540 1.386 44.874 50.118 0.023 0.017 3494.0 2955.0 1.0
mu[1] 52.684 1.633 49.787 55.840 0.027 0.019 3775.0 2789.0 1.0
sgm[0] 13.870 1.017 11.969 15.737 0.016 0.011 4026.0 3247.0 1.0
sgm[1] 16.399 1.175 14.320 18.630 0.019 0.013 4101.0 2936.0 1.0
(py312)
yasuharu@PSY-PM-PC:/mnt/d/yasuharu/www/LaCoocan/booksetc/pyda/PS2toPS3/grammar_stan/samplefiles$
警告メッセージが出されているが、4000試行中のメッセージであるので、sporadicallyであると考えられる。すなわち、問題は無いと考える。