対応のあるデータの平均値の差のベイズ分析
2変量正規分布を用いた場合
2025.11：CmdStanPy対応に改訂
Stanはベイズ分析のためのソフトであり、CmdStanPyはそのPythonインターフェースである。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
対応のあるデータ（X1, X2）が、２変量正規分布に従うとして平均値の差のベイズ分析を行った。
分析のためのStanスクリプトはリスト１のように用意した。
リスト１ ２変量正規分布のStanスクリプト（DiffMeans.stan）
data {
int N;
array[N] vector[2] Data;
}
parameters {
real mu1;
real mu2;
real<lower = 0.0001,
upper=1000> sgm1;
real<lower = 0.0001,
upper=1000> sgm2;
real<lower = -0.9999,
upper = 0.9999> rho;
}
transformed parameters {
real d;
vector[2] mus;
cov_matrix[2]
cov;
d = mu2 - mu1;
mus[1]
= mu1;
mus[2]
= mu2;
cov[1][1]
= sgm1*sgm1;
cov[1][2]
= rho*sgm1*sgm2;
cov[2][1]
= rho*sgm2*sgm1;
cov[2][2]
= sgm2*sgm2;
}
model {
mu1 ~ normal(0.0, 1000.0);
mu2 ~ normal(0.0, 1000.0);
sgm1 ~ uniform(0.0001,
1000);
sgm2 ~ uniform(0.0001,
1000);
rho ~ uniform(-0.9999,
0.9999);
for (i
in 1:N)
Data[i] ~ multi_normal(mus, cov);
}
変量X1の平均値（mu1）に対する変量X2の平均値（mu2）の変化量
変量X2の平均値（mu2） ― 変量X1の平均値（mu2）
を求めている。
リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプト例を、本ウェブサイト後半のリスト２に示した。ファイルは、Diff2VarNormal.zipにまとめた。
リスト2では、データを
data = [
[73,
70], [68, 61], [71, 68], [68, 62], [61, 60],
[73,
70], [80, 81], [71, 72], [75, 65], [59, 52],
[89,
85], [68, 71], [96, 91], [88, 83], [79, 78],
[68,
69], [71, 65], [66, 64], [72, 69], [78, 72]
]
と置いている。データ対[X1, X2]のリストである。
リスト１のStanスクリプトファイル（DiffMeans.stan）とリスト２のPythonスクリプトファイル（MainDiffMeans.py）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、次のコマンドをCmdStanPyのインストールされた環境で実行する。CmdStanPyの簡単な説明は、このウェブサイトに用意した。
(stan)
*****/Diff2VarNormal$ python MainDiffMeans.py
上のコマンドを実行すると、まず、データの散布図が表示される（図１）。
図１
図１のWindowを閉じると、Stanスクリプトのコンパイルが始まる。コンパイルは、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、サンプリング後、トレース図が表示される（図２）。
図２
左側に事後分布のグラフ、右側にMCMCサンプリングのトレース図が表示されている。
図２のWindowを閉じると、図３のグラフが表示される。
図３
左上に平均値の差mu2―mu1の事後分布のグラフが表示されている。95％CIは、負の領域にある。事後確率P(mu2＞mu1)は0.00であり、mu1の方がmu2より大きいと言える。差mu2―m1の事後分布の中央値はー3.28である。
右上のグラフは、平均値mu1とmu2の事後分布のグラフである。mu1の事後分布の方が、mu2の事後分布より右にシフトしている。
左下のグラフは、標準偏差sgm1とsgm2の事後分布のグラフである。ほとんど重なっている。
右下のグラフは、相関係数rhoの事後分布のグラフである。事後分布が、0.9を超える領域にかなりの分布があるので、強い相関があると言える。
図３のWindowを閉じると、図４のグラフが表示される。
図４
データ（X1、X2）の事後予測分布である。パラメータの事後分布に対して、データ（X1、X2）を次のスクリプトで生成している。
xy_smpls = []
for mu1, mu2, sgm1, sgm2, r in
zip(fit['mu1'], fit['mu2'],
fit['sgm1'], fit['sgm2'],fit['rho']):
mean = [mu1, mu2]
cov
= [[sgm1**2, r*sgm1*sgm2], [r*sgm2*sgm1, sgm2**2]]
xy
= ss.multivariate_normal.rvs(mean, cov)
plt.plot(xy[0], xy[1], 'o', lw=0)
xy_smpls.append(xy)
xy_smpls = np.array(xy_smpls)
ppos_xy = np.mean(xy_smpls[:,1] > xy_smpls[:,0])
赤色の直線は、Y＝Xを表しているので、この直線より上にあるデータは、Y＞Xが成り立っている。この予測事後確率は、タイトル欄に
と表示されている。５対に１対の割合でYの方がXより大きいという事後予測である。
ベイズ分析用ライブラリソフトStanのPythonインターフェースPyStan（PyStan2）の説明は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
リスト２ ２変量正規分布による対応のあるデータのベイズ分析（MainDiffMeans.py）
import scipy.stats
as ss
import numpy
as np
import matplotlib.pyplot
as plt
from cmdstanpy
import CmdStanModel
import seaborn as sb
import arviz
as az
data = [
[73,
70], [68, 61], [71, 68], [68, 62], [61, 60],
[73,
70], [80, 81], [71, 72], [75, 65], [59, 52],
[89,
85], [68, 71], [96, 91], [88, 83], [79, 78],
[68,
69], [71, 65], [66, 64], [72, 69], [78, 72]
]
plt.scatter(np.array(data)[:,0],
np.array(data)[:,1])
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Scattergram of (X,Y)', fontsize=16)
plt.tight_layout()
plt.savefig('FigScatter.png')
plt.show()
saved_files = ['FigScatter.png']
DataStan = {'N': len(data),
'Data': data}
model = CmdStanModel(stan_file='DiffMeans.stan')
fit = model.sample(data=DataStan)
print(fit.summary())
infdata = az.from_cmdstanpy(fit)
# Arviz InfereceData型へ変換
az.plot_trace(infdata, figsize=(12,8))
plt.tight_layout()
plt.savefig('FigTrace.png')
plt.show()
saved_files.append('FigTrace.png')
fit = fit.draws_pd()
# Pndas
DataFrame型への変換
plt.figure(figsize=(12,
8))
plt.subplot(221)
p025, med, p975 = np.percentile(fit['d'],
[2.5, 50, 97.5])
pposd = np.mean(fit['d']>0)
sb.kdeplot(fit['d'])
plt.title(fr'$\mu2-\mu1:$ 95%CI=[{p025:.2f}, {p975:.2f}]' + '\n' +
fr'$P(\mu2>\mu1)=${pposd:.2f}, Med.={med:.2f}', fontsize=14)
plt.yticks([])
plt.xlabel(r'$\mu2-\mu1$', fontsize=14)
plt.subplot(222)
sb.kdeplot(fit['mu1'], label = 'μ1')
sb.kdeplot(fit['mu2'], label = 'μ2')
plt.title('Posterior Distributions', fontsize=16)
plt.yticks([])
plt.xlabel('')
plt.legend(fontsize =
12)
plt.subplot(223)
sb.kdeplot(fit['sgm1'], label = 'σ1')
sb.kdeplot(fit['sgm2'], label = 'σ2')
plt.title('Posterior Distributions', fontsize=16)
plt.xlabel('')
plt.yticks([])
plt.legend(fontsize =
12)
plt.subplot(224)
sb.kdeplot(fit['rho'])
plt.title('Posterior Distribution of ρ', fontsize
= 18)
plt.yticks([])
plt.xlabel(r'$\rho$', fontsize=14)
plt.tight_layout()
plt.savefig('Figs.png')
plt.show()
saved_files.append('Figs.png')
xy_smpls = []
for mu1, mu2, sgm1, sgm2, r in
zip(fit['mu1'], fit['mu2'],
fit['sgm1'], fit['sgm2'],fit['rho']):
mean = [mu1, mu2]
cov
= [[sgm1**2, r*sgm1*sgm2], [r*sgm2*sgm1, sgm2**2]]
xy
= ss.multivariate_normal.rvs(mean, cov)
plt.plot(xy[0], xy[1], 'o', lw=0)
xy_smpls.append(xy)
xy_smpls = np.array(xy_smpls)
ppos_xy = np.mean(xy_smpls[:,1] > xy_smpls[:,0])
xy_min = np.min(xy_smpls)
xy_max = np.max(xy_smpls)
plt.plot([xy_min, xy_max], [xy_min, xy_max], c='r', lw=2,
label='Y=X')
plt.title('Predictive Distribution of (X, Y)' +
'\n' +
f'P(Y>X)={ppos_xy:.5f}', fontsize=16)
plt.xlabel('X', fontsize=14)
plt.ylabel('Y', fontsize=14)
plt.legend(fontsize=14)
plt.tight_layout()
plt.savefig('FigPredXY.png')
plt.show()
saved_files.append('FigPredXY.png')
print('\nSaved
files...')
for nm in saved_files:
print(nm)