PyStan 3によるポアッソン回帰分析
拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第12章で扱ったPyStan 2によるポアッソン回帰分析StanスクリプトをPyStan 3で実行した。Stanスクリプトは変更なしで用いることが出来たが、PythonスクリプトはPyStan 3用に書き換えた。書き換えたスクリプトを、リスト１に、Stanスクリプト（拙著、リスト12.1）をリスト２に示す。スクリプトファイルは、圧縮ファイルfiles_pssn.zipにまとめた。
リスト１では、PyStan 3によるビルドにおいてStanスクリプトがファイルから読み込んだ文字列で与えられている。すなわち
with open('Poisson.stan', 'r') as f:
sm = stan.build(f.read(),
data = Data)
である。PyStan 3では、Stanスクリプトは文字列として関数stan.build()の第１パラメータとして設定する。また、データも関数stan.build()のパラメータとして設定する。
なお、PyStan 3のインストールは、AnacondaのPython 3.10仮想環境では
pip install
pystan
であるが、Pythonスクリプトにおけるimportでは
import stan
である。
リスト１ PyStan 3によるMCMCサンプリングを行うスクリプト（PoissonReg_PS3.py）
import stan
import pickle
import numpy as np
import arviz as az
RawData = [[1, 1, 0, 0, 1, 0, 0],
[1, 1, 0, 0, 0, 1, 0],
[10, 1, 0, 0, 0, 0, 1],
[3, 0, 1, 0, 1, 0, 0],
[2, 0, 1, 0, 0, 1, 0],
[10, 0, 1, 0, 0, 0, 1],
[6, 0, 0, 1, 1, 0, 0],
[2, 0, 0, 1, 0, 1, 0],
[4, 0, 0, 1, 0, 0, 1]]
F = []
Xr1 = []
Xr2 = []
Xr3 = []
Xc1 = []
Xc2 = []
Xc3 = []
for v in RawData:
F.append(v[0])
Xr1.append(v[1])
Xr2.append(v[2])
Xr3.append(v[3])
Xc1.append(v[4])
Xc2.append(v[5])
Xc3.append(v[6])
for i in range(9):
print(' ', F[i], ' ', Xr1[i], ' ', Xr2[i], ' ', Xr3[i],
' ', Xc1[i], ' ', Xc2[i], ' ', Xc3[i])
N = 9
Data = {'N': N, 'F': F, 'Xr1': Xr1,
'Xr2': Xr2, 'Xr3': Xr3,
'Xc1': Xc1, 'Xc2': Xc2, 'Xc3': Xc3}
with open('Poisson.stan', 'r') as f:
sm = stan.build(f.read(),
data = Data)
fit = sm.sample()
print(fit)
infrnc_d = az.from_pystan(posterior =
fit, posterior_model = sm)
print('\nSummary...\n',
az.summary(infrnc_d))
with open('infrnc_d.pkl', 'wb') as f:
pickle.dump(infrnc_d, f)
d_frame = fit.to_frame()
print(d_frame.keys())
with open('d_frame.pkl', 'wb') as f:
pickle.dump(d_frame, f)
リスト１では、次のコード
with open('Poisson.stan', 'r') as f:
sm = stan.build(f.read(),
data = Data)
によりStanスクリプトのビルドを行い、MCMCサンプリングは次のコードで行っている。
fit = sm.sample()
MCMCサンプリングの結果は、コード
infrnc_d = az.from_pystan(posterior =
fit, posterior_model = sm)
print('\nSummary...\n',
az.summary(infrnc_d))
with open('infrnc_d.pkl', 'wb') as f:
pickle.dump(infrnc_d, f)
により、FitオブジェクトをInferenceDataオブジェクトinfrnc_dにまとめ、pickleでファイルに保存して、次のリスト３で用いる。
また、InferenceDataオブジェクトに纏められたものから関数az.summary()によってMCMCサンプリングの統計量を求めている。
次のコードでは、MCMCサンプリングの結果をPandas DataFrameオブジェクトにまとめて、pickleで保存している。Stanスクリプトにおけるパラメータ名とDataFrameでのキーとの対応は、print文の出力で確認できる。
d_frame = fit.to_frame()
print(d_frame.keys())
with open('d_frame.pkl', 'wb') as f:
pickle.dump(d_frame, f)
data{
int N;
int F[N];
int Xr1[N];
int Xr2[N];
int Xr3[N];
int Xc1[N];
int Xc2[N];
int Xc3[N];
}
parameters{
real mu;
real a2;
real a3;
real b2;
real b3;
real g22;
real g23;
real g32;
real g33;
}
transformed parameters{
real Lmbd[N];
real LogLmbd[N];
for (i in 1:N)
LogLmbd[i] = mu + a2 * Xr2[i] + a3 * Xr3[i]
+ b2 * Xc2[i] + b3 * Xc3[i]
+ g22 * Xr2[i] * Xc2[i] + g23 * Xr2[i] * Xc3[i]
+ g32 * Xr3[i] * Xc2[i] + g33 * Xr3[i] * Xc3[i];
for (i in 1:N)
Lmbd[i] = exp(LogLmbd[i]);
}
model{
for (i in 1:N)
F[i]
~ poisson(Lmbd[i]);
}
リスト１の実行後、リスト３を実行する。ただし、Windowsの場合は、リスト１は直接Windows上では実行できない。WindowsにWSL（Windows
Subsystem for Linux）を用意して、その上にインストールされたUbuntu上でPyStan3を用いる。詳しくはこのウェブサイトを参照されたい。PyStan 3は、直接Windows上で実行できないので、リスト１とリスト３に分けたが、Windows以外では分ける必要はない。
リスト３ リスト１のスクリプトに続いて実行するスクリプト（AnalPoissonReg_PS3.py）
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
import numpy as np
import arviz as az
import pandas as pd
pd.set_option('display.max_rows',
10000)
# 統計量表示サイズ行数の制限
pd.set_option('display.max_columns',
10000) # 統計量表示サイズ列数の制限
with open('infrnc_d.pkl', 'rb') as f:
infrnc_d = pickle.load(f)
with open('results.txt', 'w') as f:
f.write(az.summary(infrnc_d).__str__())
print('Summary...\n',
az.summary(infrnc_d, fmt = 'wide'))
for s in az.summary(infrnc_d):
print(s)
with open('d_frame.pkl', 'rb') as f:
d_frame = pickle.load(f)
print(d_frame.describe(percentiles =
[0.25, 0.5, 0.75]).T)
sb.kdeplot(d_frame['a2'], label =
'a2')
sb.kdeplot(d_frame['a3'], label =
'a3')
plt.xlabel('a2/a3', fontsize = 14)
plt.legend()
plt.savefig('Figa2a3.png')
plt.show()
sb.kdeplot(d_frame['b2'], label =
'b2')
sb.kdeplot(d_frame['b3'], label =
'b3')
plt.xlabel('b2/b3', fontsize = 14)
plt.legend()
plt.savefig('Figb2b3.png')
plt.show()
sb.kdeplot(d_frame['g22'], label =
'g22')
sb.kdeplot(d_frame['g23'], label =
'g23')
sb.kdeplot(d_frame['g32'], label =
'g32')
sb.kdeplot(d_frame['g33'], label =
'g33')
plt.xlabel('g**', fontsize = 14)
plt.legend()
plt.savefig('Figg__.png')
plt.show()
リスト３は、拙著のリスト12.2の後半に対応する処理を行うものである。リスト３を実行すると、まず、コード
with open('d_frame.pkl', 'rb') as f:
d_frame = pickle.load(f)
print(d_frame.describe(percentiles =
[0.25, 0.5, 0.75]).T)
sb.kdeplot(d_frame['a2'], label =
'a2')
sb.kdeplot(d_frame['a3'], label =
'a3')
plt.xlabel('a2/a3', fontsize = 14)
plt.legend()
plt.savefig('Figa2a3.png')
plt.show()
の実行により図１のグラフが表示される。
図１
主効果のパラメータa2とa3の事後分布が正の領域に寄っているので、大学受験時に「心理学≒臨床心理学」と考えていた学生において、心理学をはじめて知った時期が、小学生より中学生、高校生のときである学生の方が多いという傾向が反映されている。また、「心理学≒臨床心理学」と考えていた学生においては、中学生より高校生のときに心理学を知ったという学生の方が多いということが、a2の分布よりa3の分布の方が正の方向に寄っていることに反映されている。
図１のフォームを閉じると、次のコード
sb.kdeplot(d_frame['b2'], label =
'b2')
sb.kdeplot(d_frame['b3'], label =
'b3')
plt.xlabel('b2/b3', fontsize = 14)
plt.legend()
plt.savefig('Figb2b3.png')
plt.show()
により、図２のグラフが表示される。
図２
パラメータb2の分布は０を中心とするものであり、b3の分布は正の領域にあることは、「心理学を知った時期が小学生のとき」という学生においては、大学受験のとき「心理学とは臨床心理学だと思っていた」学生と、「心理学は知らなかった」という学生の数は同じであり、これらの学生に対して「臨床心理学は心理学の一部である」と理解していた学生の方が多いという集計を反映している。
図２のフォームを閉じると、次のコード
sb.kdeplot(d_frame['g22'], label =
'g22')
sb.kdeplot(d_frame['g23'], label =
'g23')
sb.kdeplot(d_frame['g32'], label =
'g32')
sb.kdeplot(d_frame['g33'], label =
'g33')
plt.xlabel('g**', fontsize = 14)
plt.legend()
plt.savefig('Figg__.png')
plt.show()
により、図３のグラフが表示される。
図３
パラメータg33の分布が負の領域にあることは、「心理学を知った時期が高校生のとき」であり、大学受験時に「臨床心理学が心理学の一部である」と理解していた学生が、他の学生に比べて少数派であるという集計結果を反映するものである。
図３のフォームを閉じると、スクリプトの実行終了である。