２変量正規分布によるベイズ分析
StanとPython
２つの変量の対でデータが得られた場合を考える。２０人の学生の国語と英語の点数のデータが図１のようであったとする。
図１ Excelデータ（TwoSubj.xlsx）
このデータに対して2変量正規分布モデルによるベイズ分析を以下のように行ってみた。
２変量
と
が次の２変量正規分布に従うとする。
上のモデルに対するStanスクリプトを本ウェブサイト後半に示すリスト１のように作成した。スクリプトファイルなどは、圧縮ファイルfilesTwoVars_nfd.zipにまとめているので、ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトによってベイズ分析を行うPythonスクリプトを本ウェブサイト後半に示すリスト２のように作成した。
リスト２のスクリプトを実行すると、まず、入力データファイル名の設定が求められる（図２）。
図２
図２では、図１のデータファイル名が設定されている。
入力データファイルはExcel形式のファイルを図１の形式で用意する。すなわち、１行目に変数名を置き、データは２行目以降に書く。１列目はデータの通し番号であるが、文字列でもよい。２列目と３列目にそれぞれの変数のデータ値を書く。
入力データファイル名を設定して、「Enter」キーを押すと、次に画像保存用ファイル名の設定が求められる（図２）。事後分布のグラフを保存するファイルの名前であるが、拡張子が「.png」の任意のファイル名でよい。複数の画像ファイルが保存されるが、設定ファイル名の「.png」の部分を「_mu1.png」などと変更した名前のファイルで保存される。図２では、「MyFigure.png」と設定されている。
画像保存ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、処理が始まる。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 pandas as pd
import numpy as np
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
finnm = input('入力データファイル（*.xlsx）= ')
foutnm = input('画像保存ファイル（*.png）= ')
#
# Excelファイルの読み込み
# モジュール xlrd がインストール済みであること
#
xlsx_fl = pd.ExcelFile(finnm)
# ExcelファイルからExcelFile型オブジェクトを生成
data_xlsx = pd.read_excel(xlsx_fl, header = None) # ExcelFile型オブジェクトからDataFrame型オブジェクトを生成
data = data_xlsx.values
# Dataframe型オブジェクトからnumpyの配列オブジェクトを生成
print('\n', data)
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)
SavedFiles = []
# 日本語表記の設定
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] + '_mu.png')
SavedFiles.append(foutnm[:-4] +
'_mu1.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')
SavedFiles.append(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')
SavedFiles.append(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')
SavedFiles.append(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')
SavedFiles.append(foutnm[:-4] +
'_rho.png')
plt.show()
print('\n保存ファイル：')
for nm in SavedFiles:
print(' ', nm)
print()