２変量正規分布によるベイズ分析
StanとPython
2025.11改訂：CmdStanPy対応
２つの変量の対でデータが得られた場合を考える。２０人の学生の国語と英語の点数のデータが図１のようであったとする。
図１ Excelデータ（TwoSubj.xlsx）
このデータに対して2変量正規分布モデルによるベイズ分析を以下のように行ってみた。
２変量
と
が次の２変量正規分布に従うとする。
上のモデルに対するStanスクリプトを本ウェブサイト後半に示すリスト１のように作成した。スクリプトファイルなどは、圧縮ファイルfilesTwoVars_nfd.zipにまとめているので、ダウンロード解凍すれば自由に利用できる。
リスト１のStanスクリプトによってベイズ分析を行うPythonスクリプトを本ウェブサイト後半に示すリスト２のように作成した。
リスト１のStanスクリプトファイル（TwoVarsNormal.stan）、リスト２のPythonスクリプトファイル（MainTwoVars.py）、および入力データファイル（TwoSubj.xlsx,図１）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境において次のコマンドを実行する。なお、CmdStanPyの簡単な説明をこのウェブサイトに用意した。
(stan) *****/filesTwoVars_nfd$ python MainTwoVars.py
上のコマンドを実行すると、入力ファイルデータ名の設定が求められる。
(stan) *****/filesTwoVars_nfd$
python MainTwoVars.py
入力データファイル（*.xlsx）= TwoSubj.xlsx
上の例では、図１のファイル名TwoSubj.xlsxが設定されている。
入力データファイルはExcel形式のファイルを図１の形式で用意する。すなわち、１行目に変数名を置き、データは２行目以降に書く。１列目はデータの通し番号であるが、文字列でもよい。２列目と３列目にそれぞれの変数のデータ値を書く。
入力データファイル名を設定して、「Enter」キーを押すと、次に画像保存用ファイル名の設定が求められる（図２）。
(stan) *****/filesTwoVars_nfd$
python MainTwoVars.py
入力データファイル（*.xlsx）= TwoSubj.xlsx
画像保存ファイル（*.png）= MyFig.png
事後分布のグラフを保存するファイルの名前であるが、拡張子が「.png」の任意のファイル名でよい。複数の画像ファイルが保存されるが、設定ファイル名の「.png」の部分を「_mu1.png」などと変更した名前のファイルで保存される。上の例では、「MyFig.png」と設定されている。
画像保存ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、処理が始まる。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.0001,
upper=10000> sgm1;
real<lower = 0.0001,
upper=10000> sgm2;
real<lower = -0.9999,
upper = 0.9999> 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.0001,
10000);
sgm2 ~ uniform(0.0001,
10000);
rho ~ uniform(-0.9999,
0.9999);
for (i
in 1:N)
X[i] ~ multi_normal(Mu, Sgm);
}
リスト２ リスト1のStanスクリプトを用いるPythonスクリプト（MainTwoVars.py）
import pandas as pd
import numpy
as np
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import japanize_matplotlib
import seaborn as sb
finnm = input('入力データファイル（*.xlsx）= ')
foutnm = input('画像保存ファイル（*.png）= ')
#
# Excelファイルの読み込み
#
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}
model = CmdStanModel(stan_file='TwoVarsNormal.stan')
fit = model.sample(data=DataStan)
print(fit.diagnose())
print(fit.summary())
fit = fit.draws_pd()
# Pandas DataFrame型に変換
SavedFiles = []
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()