２変量正規分布によるベイズ分析
StanとPython
２つの変量の対でデータが得られた場合を考える。２０人の学生の国語と英語の点数のデータが図１のようであったとする。
図１ Excelデータ（TwoSubj.xlsx）
このデータに対して2変量正規分布モデルによるベイズ分析を以下のように行ってみた。
２変量
と
が次の２変量正規分布に従うとする。
上のモデルに対するStanスクリプトを本ウェブサイト後半に示すリスト１のように作成した。スクリプトファイルなどは、圧縮ファイルfilesTwoVars.zipにまとめているので、ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトによってベイズ分析を行うPythonスクリプトを本ウェブサイト後半に示すリスト２のように作成した。
リスト２のスクリプトを実行すると、まず図２の入力データファイル名設定のダイアログウィンドウが表示される。図２では、図１のデータファイル名が設定されている。
図２
入力データファイルはExcel形式のファイルを図１の形式で用意する。すなわち、１行目に変数名を置き、データは２行目以降に書く。１列目はデータの通し番号であるが、文字列でもよい。２列目と３列目にそれぞれの変数のデータ値を書く。
入力データファイル名を設定して、「開く」ボタンをクリックすると、次に画像保存用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
事後分布のグラフを保存するファイルの名前であるが、拡張子が「.png」の任意のファイル名でよい。複数の画像ファイルが保存されるが、１つは設定したファイル名のファイルであり、他のファイルは設定ファイル名の「.png」の部分を「_mu2.png」などと変更した名前のファイルで保存される。図３では、「MyFigure.png」と設定されている。
出力ファイル名を設定して「保存」ボタンをクリックすると、データを読み込み、処理が始まる。Stanスクリプトのコンパイルとビルドに多少時間が掛かるが、その後、MCMCサンプリングが始まる。サンプリング終了後、まず、
の事後分布が表示される（図４）。
図４
図４のウィンドウを閉じると、
の事後分布が表示される（図５）。
図５
図５のウィンドウを閉じると、
の事後分布が表示される（図６）。
図６
図６のウィンドウを閉じると、
の事後分布が表示される（図７）。
図７
図７のウィンドウを閉じると、相関係数
の事後分布が表示される（図８）。
図８
図８のウィンドウを閉じると、プログラムの実行終了である。
リスト１ ２変量正規分布モデルのStanスクリプト（TwoVarsNormal.stan）
data {
int
N;
vector[2] X[N];
}
parameters {
real mu1;
real mu2;
real<lower = 0.0> sgm1;
real<lower = 0.0>
sgm2;
real<lower = -1.0, upper
= 1.0> rho;
}
transformed parameters {
vector[2] Mu;
cov_matrix[2]
Sgm;
Mu[1] = mu1;
Mu[2] = mu2;
Sgm[1][1]
= sgm1 * sgm1;
Sgm[1][2]
= rho * sgm1 * sgm2;
Sgm[2][1]
= Sgm[1][2]; // rho * sgm2 * sgm1;
Sgm[2][2]
= sgm2 * sgm2;
}
model {
mu1 ~ normal(0.0, 10000.0);
mu2 ~ normal(0.0, 10000.0);
sgm1 ~ exponential(0.0001);
sgm2 ~ exponential(0.0001);
rho ~ uniform(-1.0, 1.0);
for (i
in 1:N)
X[i] ~ multi_normal(Mu, Sgm);
}
リスト２ リスト1のStanスクリプトを用いるPythonスクリプト（MainTwoVars.py）
import tkinter
as tk
import tkinter.filedialog
as fd
import pandas as pd
import numpy
as np
import pystan
import matplotlib.pyplot
as plt
import seaborn
as sb
def askinfilename(msg = None, types = [('', '*.*')]):
"""
入力用ファイル名の設定
"""
rt
= tk.Tk()
rt.withdraw()
flnm
= fd.askopenfilename(title = msg,
filetypes = types)
rt.destroy()
return flnm
def askoutfilename(msg = None, types = [('', '*.*')], deftype
= '.txt'):
"""
出力用ファイル名の設定
"""
rt
= tk.Tk()
rt.withdraw()
flnm
= fd.asksaveasfilename(title = msg,
filetypes = types,
defaultextension = deftype)
rt.destroy()
return flnm
finnm = askinfilename('入力データExcelファイル', types = [('Excel(*.xlsx)', '*.xlsx')])
foutnm = askoutfilename(msg = '画像保存ファイル', types = [('画像 (*.png)',
'*.png')],
deftype = '.png')
xlsx_fl = pd.ExcelFile(finnm)
data_xlsx = pd.read_excel(xlsx_fl, header
= None)
data = data_xlsx.values
VarName1 = data[0, 1]
VarName2 = data[0, 2]
print(VarName1, VarName2)
X = data[1:, 1:3]
X = [[v[0], v[1]] for v in X]
print(X)
DataStan = {'N': len(X),
'X': X}
sm = pystan.StanModel(file
= 'TwoVarsNormal.stan')
fit = sm.sampling(data
= DataStan, n_jobs = 1)
print(fit)
# 日本語表記の設定
plt.rcParams['font.family']
= ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
mu1_025, mu1_50, mu1_975 = np.percentile(fit['mu1'], [2.5, 50, 97.5])
plt.title(('{0:}の平均値の事後分布'.format(VarName1)) + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
mu1_50, mu1_025, mu1_975), fontsize = 20)
sb.kdeplot(fit['mu1'])
plt.xlabel(VarName1 + 'の平均値', fontsize
= 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '.png')
plt.show()
mu2_025, mu2_50, mu2_975 = np.percentile(fit['mu2'], [2.5, 50, 97.5])
plt.title(('{0:}の平均値の事後分布'.format(VarName2)) + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
mu2_50,
mu2_025, mu2_975), fontsize = 20)
sb.kdeplot(fit['mu2'])
plt.xlabel(VarName2 + 'の平均値', fontsize
= 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_mu2.png')
plt.show()
sgm1_025, sgm1_50, sgm1_975 = np.percentile(fit['sgm1'], [2.5, 50, 97.5])
plt.title(('{0:}の標準偏差の事後分布'.format(VarName1)) + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
sgm1_50, sgm1_025, sgm1_975), fontsize = 20)
sb.kdeplot(fit['sgm1'])
plt.xlabel(VarName1 + 'の標準偏差', fontsize
= 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_sgm1.png')
plt.show()
sgm2_025, sgm2_50, sgm2_975 = np.percentile(fit['sgm2'], [2.5, 50, 97.5])
plt.title(('{0:}の標準偏差の事後分布'.format(VarName2)) + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
sgm2_50, sgm2_025, sgm2_975), fontsize = 20)
sb.kdeplot(fit['sgm2'])
plt.xlabel(VarName2 + 'の標準偏差', fontsize
= 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_sgm2.png')
plt.show()
rho_025, rho_50, rho_975 = np.percentile(fit['rho'], [2.5, 50, 97.5])
plt.title('相関係数の事後分布' + \
'\nMed = {0:.2f}, 95%CI = [{1:.2f}, {2:.2f}]'.format(
rho_50, rho_025, rho_975), fontsize = 20)
sb.kdeplot(fit['rho'])
plt.xlabel('相関係数', fontsize = 16)
plt.yticks([])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_rho.png')
plt.show()