対応のあるデータの平均値の差のベイズ分析
2変量正規分布
対応のあるデータ（X1, X2）が、２変量正規分布に従うとして平均値の差のベイズ分析を行った。
分析のためのStanスクリプトを以下のように用意した。
stan_code = """
data {
int
N;
vector[2] Data[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 {
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
~ exponential(0.0001);
sgm2
~ exponential(0.0001);
rho
~ uniform(-1.0, 1.0);
for
(i in 1:N)
Data[i] ~ multi_normal(mus, cov);
}
"""
変量X1の平均値に対する変量X2の平均値の変化量
変量X2の平均値 ― 変量X1の平均値
を求めている。
完全なプログラムは本ウェブサイト後ろにリスト１として示した。
リスト１では、データを
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によるベイズ分析を行うと、以下の結果を得る。リスト１の実行は、CentOS StreamにインストールしたAnacondaのJupyter Notebookにおいて行った。
MCMCサンプリングの統計量を出力するスクリプト
print(fit)
の実行結果は、以下の通りである。
Inference for Stan model:
anon_model_8d31416547bbd9a7108b8ce1b885b1ae.
4 chains, each with iter=2000;
warmup=1000; thin=1;
post-warmup draws per chain=1000,
total post-warmup draws=4000.
mean se_mean
sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu1 73.63 0.05 2.25 69.09 72.14 73.64 75.1 78.06 1727 1.0
mu2 70.32 0.05 2.29 65.77 68.81 70.34 71.82 74.84 1742 1.0
sgm1 9.72 0.04 1.68 7.02 8.54 9.5 10.69 13.6 1523 1.0
sgm2 9.92 0.04 1.73 7.16 8.72 9.69 10.88 14.01 1586 1.0
rho
0.92 1.0e-3 0.04 0.82 0.9 0.93 0.95 0.97 1537 1.0
d
-3.3 0.01 0.87 -5.03 -3.87 -3.3 -2.73 -1.58 4449 1.0
mus[1] 73.63 0.05 2.25 69.09 72.14 73.64 75.1 78.06 1727 1.0
mus[2] 70.32 0.05 2.29 65.77 68.81 70.34 71.82 74.84 1742 1.0
cov[1,1] 97.28 0.91 35.33 49.22 72.9 90.31 114.33 185.01 1497 1.0
cov[2,1] 91.62 0.93 34.96 44.12 67.96 85.03 108.62 182.94 1420 1.0
cov[1,2] 91.62 0.93 34.96 44.12 67.96 85.03 108.62 182.94 1420 1.0
cov[2,2] 101.31 0.95 37.29 51.29 76.07 93.91 118.29 196.42 1556 1.0
lp__ -88.2 0.04 1.64 -92.21 -89.06 -87.87 -86.99
-86.01 1406 1.0
Rhatの値が1.0であるので、サンプリングの収束はよいと言える。
平均値の差
の事後分布は、図１のようである。
図１
95% CIが０より小さい領域にあるので、差（平均値
―平均値
）は負である、すなわち、
平均値＜平均値
であると考えられる。
平均値と平均値の事後分布は図２のようである。かなりの重なりが見られる。
図２
変量X1とX2のそれぞれの標準偏差
と
の事後分布は図３のようである。
図３
相関係数の事後分布は、図４のようである。
図４
ベイズ分析用ライブラリStanのPythonインターフェースPyStanの説明は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
import scipy.stats as ss
import numpy as np
import matplotlib.pyplot as plt
import pystan
import seaborn as sb
stan_code = """
data {
int
N;
vector[2] Data[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 {
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
~ exponential(0.0001);
sgm2
~ exponential(0.0001);
rho
~ uniform(-1.0, 1.0);
for
(i in 1:N)
Data[i] ~ multi_normal(mus, cov);
}
"""
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]
]
DataStan = {'N': len(data), 'Data':
data}
sm = pystan.StanModel(model_code =
stan_code)
fit = sm.sampling(data = DataStan)
print(fit)
p025, p975 = np.percentile(fit['d'],
[2.5, 97.5])
sb.kdeplot(fit['d'])
plt.title('μ2 - μ1\n95%CI
= [{0:.2f}, {1:.2f}]'.format(p025, p975), fontsize = 16)
plt.yticks([])
plt.savefig('d.png')
plt.show()
sb.kdeplot(fit['mu1'], label = 'μ1')
sb.kdeplot(fit['mu2'], label = 'μ2')
plt.yticks([])
plt.legend(fontsize = 12)
plt.savefig('mus.png')
plt.show()
sb.kdeplot(fit['sgm1'], label = 'σ1')
sb.kdeplot(fit['sgm2'], label = 'σ2')
plt.yticks([])
plt.legend(fontsize = 12)
plt.savefig('sgms.png')
plt.show()
sb.kdeplot(fit['rho'])
plt.title('ρ', fontsize = 18)
plt.yticks([])
plt.savefig('rho.png')
plt.show()