２変量正規分布によるベイズ分析
CmdStanPy対応：2025.11改訂
２つの変量の対でデータが得られた場合を考える。２０人の学生の国語と英語の点数のデータが図１のようであったとする。
図１ Excelデータ（TwoSubj.xlsx）
このデータに対して2変量正規分布モデルによるベイズ分析を以下のように行ってみた。
２変量
と
が次の２変量正規分布に従うとする。
上のモデルに対するStanスクリプトを本ウェブサイト後半に示すリスト１のように作成した。スクリプトファイルなどは、圧縮ファイルfilesTwoVars.zipにまとめているので、ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトによってベイズ分析を行うPythonスクリプトを本ウェブサイト後半に示すリスト２のように作成した。
リスト１のStanスクリプトファイル（TwoVarsNormal.sta）、リスト２のPythonスクリプトファイル（MainTwoVars.py）、入力データファイル（図１のファイルの場合は、TwoSubj.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、、CmdStanPyのインストールされた環境において、リスト２のPythonスクリプトを実行する。
(stan) *****/filesTwoVars$ python MainTwoVars.py
上のコマンドを実行すると、まず図２の入力データファイル名設定のダイアログウィンドウが表示される。図２では、図１のデータファイル名が設定されている。
図２
入力データファイルはExcel形式のファイルを図１の形式で用意する。すなわち、１行目に変数名を置き、データは２行目以降に書く。１列目はデータの通し番号であるが、文字列でもよい。２列目と３列目にそれぞれの変数のデータ値を書く。
入力データファイル名を設定して、「Open」ボタンをクリックすると、次に画像保存用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
事後分布のグラフを保存するファイルの名前であるが、拡張子が「.png」の任意のファイル名でよい。複数の画像ファイルが保存されるが、１つは設定したファイル名のファイルであり、他のファイルは設定ファイル名の「.png」の部分を「_mu2.png」などと変更した名前のファイルで保存される。図３では、「MyFig.png」と設定されている。
出力ファイル名を設定して「Save」ボタンをクリックすると、データを読み込み、処理が始まる。Stanスクリプトのコンパイルに多少時間が掛かるが、その後、MCMCサンプリングが始まる。サンプリング終了後、まず、
の事後分布が表示される（図４）。
図４
図４のWindowを閉じると、
の事後分布が表示される（図５）。
図５
図５のWindowを閉じると、
の事後分布が表示される（図６）。
図６
図６のWindowを閉じると、
の事後分布が表示される（図７）。
図７
図７のWindowを閉じると、相関係数
の事後分布が表示される（図８）。
図８
図８のWindowを閉じると、スクリプトの実行終了である。
リスト１ ２変量正規分布モデルのStanスクリプト（TwoVarsNormal.stan）
data {
int N;
array[N] vector[2] X;
}
parameters {
real mu1;
real mu2;
real<lower = 0.01,
upper=10000> sgm1;
real<lower = 0.01,
upper=10000> sgm2;
real<lower = -0.999,
upper = 0.999> 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];
Sgm[2][2]
= sgm2 * sgm2;
}
model {
mu1 ~ normal(0.0, 10000.0);
mu2 ~ normal(0.0, 10000.0);
sgm1 ~ uniform(0.01,
10000);
sgm2 ~ uniform(0.01,
10000);
rho ~ uniform(-0.999,
0.999);
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
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import japanize_matplotlib
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}
model = CmdStanModel(stan_file='TwoVarsNormal.stan')
fit = model.sample(data=DataStan)
print(fit.summary())
fit = fit.draws_pd()
# Pandas DataFrame型に変換
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()