PythonでStanを使う
Python用のStan（PyStan）を使ってみた。１変量データに正規分布モデルを適用して平均と標準偏差の母数パラメータ値の事後分布を求めるものと、２変量データに多変量正規分布モデルを適用して平均、標準偏差と相関係数を求めるものの２例を作成してみた。、まず、１変量データの場合について説明する。その後、２変量データの場合について説明する。
PyStanのインストールについては、別のウェブサイトで説明している。
PyStan 3対応は、別のウェブサイトを用意した。
なお、岡本安晴（著）「いまさら聞けないPythonでデータ分析―多変量解析，ベイズ統計分析（PyStan，PyMC）―」でStanの説明を行っている。
１変量データの場合
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import pystan
import matplotlib.pyplot as plt
import numpy as np
OneVar_code = """
data {
int<lower=0> J;
real y[J];
}
parameters {
real<lower=-1000, upper=1000> mu;
real<lower=0, upper=1000> sigma;
}
transformed parameters {
}
model {
y ~
normal(mu, sigma);
}
"""
LData = [ 44, 44, 56, 59, 44, 44, 63, 60, 34,
58,
55, 50, 32, 61, 47, 25, 38, 47, 30, 42]
OneVar_dat = { 'J' : len(LData), 'y':
LData }
sm =
pystan.StanModel(model_code=OneVar_code)
fit = sm.sampling(data=OneVar_dat,
iter=5000, chains=4, n_jobs=1)
print(fit)
LResults = fit.extract(permuted=True)
Lmu = LResults['mu']
meanMu = np.mean(Lmu)
plt.figure()
plt.title("Posterior Distribution
of mu Mean =
{0:.7}".format(meanMu))
plt.hist(Lmu, bins = 50)
plt.show()
Lsigma = LResults['sigma']
meanSigma = np.mean(Lsigma)
plt.figure()
plt.title("Posterior Distribution
of sigma Mean =
{0:.7}".format(meanSigma))
plt.hist(Lsigma, bins = 50)
plt.show()
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Stanのコードは、文字列データとして用意する。上のコードでは、文字列データとして用意したものをOneVar_codeに設定している。データとパラメータの値は、辞書型のオブジェクトにそれぞれキーとオブジェクトのペアとして用意して、変数OneVar_datに設定している。これらの準備をして、Stanを以下のように呼び出している。Windowsで実行するときは、「n_jobs=1」の設定が必要である。
sm =
pystan.StanModel(model_code=OneVar_code)
fit =
sm.sampling(data=OneVar_dat, iter=5000, chains=4, n_jobs=1)
計算結果は変数fitに返される。上のプログラムでは、fitに返された結果を、print()で出力し、さらにnumpyで平均値を求め、matplotlibで事後分布のヒストグラムを描画している。
print文による出力は図１に示すようになっている。
図１
muの事後分布は図２のように示される。
図２
図２のウィンドウの右上のＸ印をクリックして閉じると、図３のsigmaの事後分布が表示される。
図３
図３の右上のＸ印のクリックでウィンドウは閉じられてプログラムの実行終了となる。
上のプログラムは圧縮ファイルSingleFile.zipにまとめておいた。ファイル名SingleFile.zipをクリックしてダウンロード・展開すれば利用できる。
上のプログラムでは、コードが１つのファイルにまとめて書かれている。これをStanの部分とPythonの部分に分けて用意して、Pythonの実行時にStanのコードの書かれたファイルを読み込むこともできる。上のコードを２つのファイルに分けて用意したものが以下に示すものである。
まず、Stanのコードを次のように用意する。ファイル拡張子が .stan であることに注意。名前は適当に付ける。また、上のプログラム例では、Stanのコードは文字列データとして用意されていたが、別ファイルとして用意したOneVar.stanでは、文字列データではなく、普通のテキストファイルのコードとして書かれていることに注意。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝ OneVar.stan ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
data {
int<lower=0> J;
real
y[J];
}
parameters {
real<lower=-1000, upper=1000> mu;
real<lower=0, upper=1000> sigma;
}
transformed parameters {
}
model {
y ~
normal(mu, sigma);
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のStanコードを読み込んで実行するPythonコードを次のように用意した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import pystan
import matplotlib.pyplot as plt
import numpy as np
LData = [ 44, 44, 56, 59, 44, 44, 63, 60, 34,
58,
55, 50, 32, 61, 47, 25, 38, 47, 30, 42]
OneVar_dat = { 'J' : len(LData), 'y':
LData }
sm =
pystan.StanModel(file='OneVar.stan')
fit = sm.sampling(data=OneVar_dat,
iter=5000, chains=4, n_jobs=1)
print(fit)
LResults = fit.extract(permuted=True)
Lmu = LResults['mu']
meanMu = np.mean(Lmu)
plt.figure()
plt.title("Posterior Distribution
of mu Mean =
{0:.7}".format(meanMu))
plt.hist(Lmu, bins = 50)
plt.show()
Lsigma = LResults['sigma']
meanSigma = np.mean(Lsigma)
plt.figure()
plt.title("Posterior Distribution
of sigma Mean =
{0:.7}".format(meanSigma))
plt.hist(Lsigma, bins = 50)
plt.show()
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Stanの呼び出し
sm =
pystan.StanModel(file='OneVar.stan')
において、Stanのコードが「file='OneVar.stan'」とその格納されているファイル名で指定されていることに注意。
上のPythonコードの実行結果は、１つのファイルにコードがまとめられているプログラムの実行結果を示す図１から図３までのものと同じである。StanコードとPythonコードの２つのファイルに分けられたプログラムのファイルは圧縮ファイルMultiFiles.zipにまとめた。ファイル名MultiFiles.zipをクリックダウンロード・展開すれば利用できる。
多変量正規分布の簡単な場合として、２変量正規分布について考える。
まず、２変量正規分布モデルを使うStanスクリプトを以下のように作成した。スクリプトファイルは、圧縮ファイルSampleFiles.zipにまとめているので、ダウンロード解凍すれば自由に利用できる。
＝＝＝＝＝＝＝＝＝＝＝＝ ファイル名：TwoVarNormal.stan ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
data {
int
N;
vector[2]
Data[N];
real
l1;
real
u1;
real
l2;
real
u2;
real
sgm_l1;
real
sgm_u1;
real
sgm_l2;
real
sgm_u2;
}
parameters {
real<lower
= l1, upper = u1> mu1;
real<lower
= l2, upper = u2> mu2;
real<lower
= sgm_l1, upper = sgm_u1> sgm1;
real<lower
= sgm_l2, upper = sgm_u2> sgm2;
real<lower
= -1.0, upper = 1.0> r;
}
transformed parameters {
vector[2]
mu;
matrix[2,
2] cov;
mu[1]
= mu1;
mu[2]
= mu2;
cov[1][1]
= sgm1 * sgm1;
cov[1][2]
= sgm1 * sgm2 * r;
cov[2][1]
= sgm1 * sgm2 * r;
cov[2][2]
= sgm2 * sgm2;
}
model {
mu1
~ uniform(l1, u1);
mu2
~ uniform(l2, u2);
sgm1
~ uniform(sgm_l1, sgm_u1);
sgm2
~ uniform(sgm_l2, sgm_u2);
r
~ uniform(-1.0, 1.0);
for
(i in 1:N){
Data[i]
~ multi_normal(mu, cov);
}
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のStanスクリプト（ファイル名：TawVarNormal.stan）を使用するPythonスクリプト例を以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝ ファイル名：Sample2VarNorm.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
Data = [
[179, 121], [185, 89], [181, 106], [210,
100], [202, 104],
[197, 91], [227, 107], [208, 101], [188,
99], [189, 112],
[240, 111], [203, 103], [202, 85], [190,
95], [221, 96],
[150, 83], [205, 103], [177, 90], [221,
92], [172, 86],
[190, 87], [218, 101], [205, 100], [208,
120], [224, 98],
[199, 116], [203, 96], [213, 91], [200,
109], [190, 100],
[211, 115], [204, 85], [185, 92], [178,
86], [231, 122],
[206, 84], [197, 93], [236, 93], [213,
106], [205, 103],
[212, 102], [177, 86], [193, 96], [183,
93], [204, 87],
[203, 108], [203, 103], [180, 100], [196,
101], [167, 75],
[202, 106], [182, 98], [225, 109], [219,
115], [205, 102],
[200, 99], [219, 94], [189, 102], [220,
86], [186, 102],
[203, 101], [210, 124], [254, 103], [208,
95], [192, 106],
[176, 91], [173, 97], [195, 100], [189,
98], [220, 115],
[202, 112], [257, 108], [207, 104], [199,
96], [173, 98],
[230, 101], [230, 111], [212, 105], [196,
82], [183, 82],
[230, 103], [207, 112], [225, 101], [192,
104], [175, 83],
[223, 105], [192, 106], [202, 87], [204,
87], [184, 97],
[193, 94], [204, 109], [196, 92], [216,
109], [219, 113],
[256, 117], [210, 86], [217, 103], [206,
87], [213, 105] ]
mu = np.mean(Data, axis = 0)
print('mu =\n', mu)
cov = np.cov(np.array(Data).T)
print('cov =\n', cov)
sgms = [cov[0][0]**0.5,
cov[1][1]**0.5]
print('sgm =\n', sgms)
L1 = mu[0] - 5 * sgms[0]
U1 = mu[0] + 5 * sgms[0]
L2 = mu[1] - 5 * sgms[1]
U2 = mu[1] + 5 * sgms[1]
sgm_L1 = 0.0
sgm_U1 = 10 * sgms[0]
sgm_L2 = 0.0
sgm_U2 = 10 * sgms[1]
def f_init():
return
{'mu1': mu[0], 'mu2': mu[1], 'sgm1': sgms[0], 'sgm2': sgms[1],
'r':
0.0}
sm = pystan.StanModel(file =
'TwoVarNormal.stan')
stan_data = {'N': len(Data), 'Data':
Data, 'l1': L1, 'u1': U1, 'l2': L2, 'u2': U2,
'sgm_l1':
sgm_L1, 'sgm_u1': sgm_U1,
'sgm_l2':
sgm_L2, 'sgm_u2': sgm_U2}
fit = sm.sampling(data = stan_data,
init = f_init, n_jobs = 1)
print(fit)
samples_mu1 = fit['mu1']
samples_mu2 = fit['mu2']
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_mu1, label = 'mu1')
sb.kdeplot(samples_mu2, label = 'mu2')
plt.legend(fontsize = 20)
plt.show()
samples_sgm1 = fit['sgm1']
samples_sgm2 = fit['sgm2']
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_sgm1, label =
'$\sigma_1$')
sb.kdeplot(samples_sgm2, label =
'$\sigma_2$')
plt.legend(fontsize = 20)
plt.show()
samples_r = fit['r']
plt.title('Posterior Distributions',
fontsize = 20)
sb.kdeplot(samples_r, label = 'r')
plt.xlabel('r', fontsize = 20)
plt.legend(fontsize = 20)
plt.show()
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のPythonスクリプトを、以下のように実行した。
コマンドpythonで上のPythonスクリプトSample2VarNorm.pyを実行する（図C.1）。
図C.1
データの平均値、共分散行列、標準偏差が出力され、その後、Stanスクリプトのコンパイルが始まる（図C.2）。
図C.2
Stanによるサンプリングが終わると、まずパラメータmuの事後分布が表示される（図C.3）。
図C.3
図C.3のフォームを閉じると、パラメータsgmの事後分布のグラフが表示される（図C.4）。
図C.4
図C.4のフォームを閉じると、rの事後分布のグラフが表示される（図C.5）。
図C.5
図C.5のフォームを閉じると、プログラムの実行終了となる。実行終了後の端末画面は図C.6のようになっている。
図C.6
StanによるMCMCサンプリングの統計量が表示されている。
StanによるMCMCサンプリングと統計量の計算
上の例では事後分布の平均値などの統計量は、メソッドsamplingの返り値の出力で見ている。平均値とか四分位数などは出力されているが、MAP推定値は無い。