Stanによる半正規分布の作成
半正規分布（half-normal distribution）は、Stanではリスト１に示すように簡単に作成できる。
リスト１ 半正規分布の作成（halfnormal.py）
import stan
import arviz
as az
import pickle
StanScript = """
functions {
real
hnormal_lpdf(real y, real mu, real sgm) {
if (y > mu) {
return normal_lpdf(y|mu,
sgm) + log(2.0);
}
else {
return negative_infinity();
}
}
}
data {
real
mu;
real<lower=0> sgm;
}
parameters {
real
a;
real
b;
}
model {
a ~
normal(mu, sgm);
//target += hnormal_lpdf(b|mu,
sgm);
b ~ hnormal(mu, sgm);
}
"""
sm = stan.build(StanScript,data = {'mu':1.0, 'sgm':2.0} )
print('OK!')
fit = sm.sample(
num_samples = 10000)
print(fit)
print('Statistics...\n',
az.summary(az.from_pystan(posterior
= fit, posterior_model = sm)))
with open('fita.pkl',
'wb') as f:
pickle.dump(fit['a'],
f)
with open('fitb.pkl',
'wb') as f:
pickle.dump(fit['b'],
f)
正規分布の平均より右側を取り出し、関数値を２倍（スクリプトではlog(2.0)を足している。これは密度関数の対数値であることによる）にして半正規分布の密度関数としている。
リスト１のスクリプトファイルhalfnormal.pyを次のコマンドで実行する。スクリプトファイルは、圧縮ファイルscriptfiles.zipにまとめたので、自由にダウンロードされたい。
(py310) xxxxx…xxx/halfnormalstan/scriptfiles$ python halfnormal.py
実行環境は、Windows上にインストールされたWSL（Windows Subsystem for Linux）上のUbuntuにインストールしたAnacondaのPython 3.10仮想環境である。Python 3.10仮想環境にPyStan 3をインストールして使っている。
上のコマンドでリスト１のスクリプトを実行すると、まずビルドが行われ、その後、MCMCサンプリングが始まる。
MCMCサンプリングの統計量は、次のコードで書き出している。
print('Statistics...\n',
az.summary(az.from_pystan(posterior
= fit, posterior_model = sm)))
下記のような出力であるので、MCMCサンプリングに問題は無かったと言えよう。
Statistics...
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
a
0.971 2.010 -2.886 4.634 0.019 0.013 11646.0 14065.0 1.0
b
2.597 1.217 1.000 4.772 0.015 0.010 7023.0 10306.0 1.0
正規分布と半正規分布のサンプルは、次のコードで、パラメータaとbのサンプルとして、pickleでファイルに保存している。
with open('fita.pkl',
'wb') as f:
pickle.dump(fit['a'],
f)
with open('fitb.pkl',
'wb') as f:
pickle.dump(fit['b'], f)
Windowsの場合、PyStan 3のスクリプトはLinux仮想環境で行っているが、得られたMCMCサンプルの描画は、Windows上で、リスト２のスクリプトによって行った。リスト２のファイルdrawab.pyも圧縮ファイルscriptfiles.zipにまとめた。
リスト２ リスト１のアパrメータaとbのサンプルを読み込んで表示するスクリプト（drawab.py）
import pickle
import matplotlib.pyplot
as plt
import seaborn
as sb
with open('fita.pkl',
'rb') as f:
a_smpl
= pickle.load(f)
with open('fitb.pkl',
'rb') as f:
b_smpl
= pickle.load(f)
sb.kdeplot(a_smpl[0])
sb.kdeplot(b_smpl[0])
plt.show()
plt.hist(a_smpl[0],
density = True, alpha = 0.5)
plt.hist(b_smpl[0],
density = True, alpha = 0.5)
plt.show()
リスト２のスクリプトを、Windows上にインストールしたPython 3.11 IDLEにより実行すると、次のコードによりリスト１のスクリプトの実行で保存されたファイルが読み込まれる。
with open('fita.pkl',
'rb') as f:
a_smpl
= pickle.load(f)
with open('fitb.pkl',
'rb') as f:
b_smpl
= pickle.load(f)
次のコードを実行すると、KDEによるサンプルの密度関数が表示される（図１）。
sb.kdeplot(a_smpl[0])
sb.kdeplot(b_smpl[0])
plt.show()
図１
図１では、半正規分布の左側が急峻な崖になっているが、裾は曲線で描かれている。これは、KDEによる描画の特徴である。
次のコードにより、ヒストグラムを描くと、半正規分布の左側が垂直な直線で描かれていることが分かる（図２）。
plt.hist(a_smpl[0],
density = True, alpha = 0.5)
plt.hist(b_smpl[0],
density = True, alpha = 0.5)
plt.show()
図２