Stan、Python、Windows、C++コンパイラ
平均値の差、および相関係数のベイズ分析
StanのPython用インターフェースとしてPyStanがあるが、Windowsなどにインストールするときは、PythonのバージョンとPyStanのバージョン、Stanスクリプトの実行時に使用されるC++コンパイラとの相性などを考える必要がある。しかし、これらの問題は、Anacondaの仮想環境内で利用するときは、インストール時にAnacondaによって自動的に調整されてC++コンパリラも含めて適切な組み合わせで仮想環境内にインストールされる。Anacondaそのものも簡単にインストールできるが、簡単なインストール法をこのウェブサイトで説明している。
なお、Windowsの場合、Anacondaの仮想環境内であっても、StanによるMCMCサンプリングにおいては
n_jobs = 1
の制約があることに注意。
Windowsの場合、Pythonをウェブサイト
https://www.python.org/downloads/
からダウンロード・インストールできるが、このときIDLE（統合開発環境）もインストールされる。Anacondaの仮想環境内においても、このIDLEのエディタが利用できる。
図１では、エディタとしてPython IDLEのエディタを用い、Anacondaの仮想環境においてスクリプトファイルmain.pyを実行している様子が示されている。
図１
なお、Stan (PyStan)については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。
ここでは、以下において統計分析の基本的方法である対応のあるデータの平均値の差のベイズ分析、独立なデータの平均値の差のベイズ分析、および相関係数のベイズ分析の例を示す。スクリプトファイルは圧縮ファイルStanFiles.zipとしてまとめたので、ダウンロード解凍すれば自由に利用できる。
ある健康法を20人の人に実行してもらい、実行前と一カ月間の実行後の血圧を測ったところ表Aの結果を得たとする。
|
表A 健康法実行前と実行後の血圧mmHg。 （実行前血圧、実行後血圧）；仮想データ |
|
(133, 123),
(133, 127), (136, 129), (136, 132), (134, 127),
(136, 127), (131, 118), (139, 134), (135, 129),
(140, 130), (132, 124), (130, 122), (138, 127),
(138, 124), (130, 118), (137, 119), (131, 123),
(135, 128), (136, 120), (134, 129) |
健康法の効果を調べるために、人
の実行前の血圧
と実行後の血圧
の平均値を比べるが、次式の関係
に注目して、差
の平均値のベイズ分析を行う。個人ごとの差を見ることにより、各個人における変化に注目することができる。
StanスクリプトをリストA.1のように用意した。
リストA.1 対応のあるデータの平均値のベイズ分析Stanスクリプト（AnalPair.stan）
data {
int N;
vector[2] dt[N];
}
transformed data {
real diff[N];
for (i in 1:N)
diff[i] = dt[i][1] - dt[i][2];
}
parameters {
real mu;
real<lower = 0.0> sgm;
}
model {
mu ~ normal(0.0,
1000.0);
sgm ~ exponential(0.001);
for (i in 1:N)
diff[i] ~ normal(mu, sgm);
}
リストA.1のStanスクリプトを用いて表Aのデータを分析するPythonスクリプトをリストA.2のように作成した。
リストA.2 表Aのベイズ分析のPythonスクリプト
import matplotlib.pyplot as plt
import numpy as np
import pystan
import seaborn as sb
import scipy.stats as ss
Data = [(133, 123), (133, 127), (136,
129), (136, 132),
(134, 127), (136, 127), (131, 118), (139, 134),
(135, 129), (140, 130), (132, 124), (130, 122),
(138,
127), (138, 124), (130, 118), (137, 119),
(131, 123), (135, 128), (136, 120), (134, 129)
]
StanData = {'N': len(Data), 'dt':
Data}
sm = pystan.StanModel(file =
'AnalPair.stan')
fit = sm.sampling(data = StanData,
n_jobs = 1)
print(fit)
mu025P, mu975P =
np.percentile(fit['mu'], [2.5, 97.5])
sb.kdeplot(fit['mu'])
plt.yticks([])
plt.xlabel('mu', fontsize = 12)
plt.title('Posterior distribution of
mu' +
'\nMean = {0:.1f}, 95%CI = [{1:.1f}, {2:.1f}]'.format(
fit['mu'].mean(), mu025P, mu975P) , fontsize = 14)
plt.show()
sgm025P, sgm975P =
np.percentile(fit['sgm'], [2.5, 97.5])
sb.kdeplot(fit['sgm'])
plt.yticks([])
plt.xlabel('sigma', fontsize = 12)
plt.title('Posterior distribution of
sigma' +
'\nMean = {0:.1f}, 95%CI = [{1:.1f}, {2:.1f}]'.format(
fit['sgm'].mean(), sgm025P, sgm975P) , fontsize = 14)
plt.show()
#
# Posterior prediction of difference
#
d_post = []
for m, s in zip(fit['mu'],
fit['sgm']):
d_post.append(ss.norm.rvs(m,
s))
d_post = np.array(d_post)
p_pos = (d_post > 0.0).sum() /
len(d_post)
sb.kdeplot(d_post)
plt.xlabel('diff.', fontsize = 14)
plt.title('Posterior prediction of
difference' +
'\nProb(diff. > 0) = {0:.3f}'.format(p_pos), fontsize = 14)
plt.show()
リストA.2のスクリプトを実行すると、Stanスクリプトのコンパイルと続くMCMCサンプリングの後、図A.1のグラフが表示される。
図A.1
図A.1には差の平均の事後分布が描かれている。平均値が9.2、95％確信区間が［7.3、11.0］であるので、健康法が血圧を下げる効果が認められると考えてよさそうである。
図A.1のフォームを閉じると、図A.2のフォームが表示される。
図A.2
差の分布の標準偏差の事後分布である。図A.2のフォームを閉じると図A.3のグラフが表示される。
図A.3
個々の差の値は、平均値と標準偏差から確率変数の実現値として与えられるが、その分布が事後予測として示されている。事後予測分布における差が正である確率が約98％であるので、健康法を実行すればほぼ確実に効果が期待できると予想される。
A市から15人、B市から20人をランダムに選び、血圧を測ると表Bの結果を得たとする。
|
表B 血圧mmHg（仮想データ） |
|
|
A市 |
126, 122,
120, 132, 129, 124, 124, 125, 124, 124, 121, 121,
130, 125, 129 |
|
B市 |
126, 139,
137, 130, 139, 132, 128, 137, 134, 137, 131, 133,
136, 131, 130, 134, 141, 135, 136, 133 |
A市、B市それぞれの市民の血圧が平均muA、muB、標準偏差がsgmの正規分布に従うとする。すなわち
である。このとき、平均の差
を調べるStanスクリプトをリストB.1のように用意した。
リストB.1 独立なデータの平均の差のベイズ分析Stanスクリプト
data {
int NA;
int NB;
real DataA[NA];
real DataB[NB];
}
parameters {
real muA;
real muB;
real<lower =
0.0> sgm;
}
transformed parameters {
real mu_d;
mu_d = muA - muB;
}
model {
muA ~ normal(0.0, 1000.0);
muB ~ normal(0.0, 1000.0);
sgm ~ exponential(0.001);
for (i in 1:NA)
DataA[i] ~ normal(muA, sgm);
for (i in 1:NB)
DataB[i] ~ normal(muB, sgm);
}
表BのデータをリストB.1のStanスクリプトで分析するPythonスクリプトをリストB.2のように用意した。
リストB.2 表B.2のベイズ分析用Pythonスクリプト
import matplotlib.pyplot as plt
import numpy as np
import pystan
import seaborn as sb
import scipy.stats as ss
dataA = [126, 122, 120, 132, 129,
124, 124, 125, 124, 124,
121, 121, 130, 125, 129
]
dataB = [126, 139, 137, 130, 139,
132, 128, 137, 134, 137,
131, 133, 136, 131, 130,
134, 141, 135, 136, 133
]
StanData = {'NA': len(dataA), 'DataA':
dataA,
'NB': len(dataB), 'DataB': dataB}
sm = pystan.StanModel(file =
'AnalInd.stan')
fit = sm.sampling(data = StanData,
n_jobs = 1)
print(fit)
mean_mu_d = fit['mu_d'].mean()
mu_d025P, mu_d975P = np.percentile(fit['mu_d'],
[2.5, 97.5])
sb.kdeplot(fit['mu_d'])
plt.title('Posterior distribution of
mu_d' +
'\nmean = {0:.2f}
95%CI = [{1:.2f}, {2:.2f}]'.format(
mean_mu_d, mu_d025P, mu_d975P), fontsize = 16)
plt.yticks([])
plt.xlabel('mu_d', fontsize = 14)
plt.show()
plt.title('Posterior distributions',
fontsize = 16)
sb.kdeplot(fit['muA'], label = 'muA')
sb.kdeplot(fit['muB'], label = 'muB')
plt.yticks([])
plt.legend(fontsize = 12)
plt.show()
sb.kdeplot(fit['sgm'])
plt.title('Posterior distribution of
sgm', fontsize = 16)
plt.yticks([])
plt.show()
#
# Posterior
prediction of difference
#
diff = []
for ma, mb, s in zip(fit['muA'],
fit['muB'], fit['sgm']):
va = ss.norm.rvs(ma, s)
vb = ss.norm.rvs(mb, s)
diff.append(va - vb)
diff = np.array(diff)
p_pos_diff = (diff > 0.0).sum() /
len(diff)
sb.kdeplot(diff)
plt.title('Posterior prediction of
Differences A - B' +
'\nProb(A-B > 0) = {0:.3f}'.format(p_pos_diff), fontsize = 16)
plt.xlabel('A - B', fontsize = 14)
plt.yticks([])
plt.show()
リストB.2のスクリプトを実行すると、リストB.1のStanスクリプトのコンパイル後、MCMCサンプリングが行われ、その後、まずmu_dの事後分布が表示される（図B.1）。
図B.1
平均が−8.89、95％確信区間が［−11.49、−6.271］であり、mu_dは確信をもって負である、すなわちA市の平均値の方がB市の平均値より低いと言える。
図B.1のフォームを閉じると、図B.2のフォームが提示される。
図B.2
図B.2のフォームでは、A市B市の平均値の事後分布が描かれている。２つの事後分布は明確に分離されているので、A市の平均値の方が低いと言えそうである。
図B.2のフォームを閉じると、図B.3のフォームが提示される。
図B.3
図B.3のフォームには、標準偏差の事後分布が描かれている。
図B.3のフォームを閉じると、図B.4のフォームが提示される。
図B.4
図B.4のフォームには、2つの市、A市とB市からランダムに一人ずつサンプルを採り、血圧の差を求めたときの事後予測分布である。A市の人の方がB市の人より血圧が高い確率が約0.06である、すなわちランダムに一人ずつ選んで血圧を比べたとき、A市の人の方が血圧が高いのは20回に１回ぐらいの頻度であることが予想されている。
学生20人の国語と英語の成績が表Cのようであったとする。
|
表C 学生20人の成績（国語、英語）＜仮想データ＞ |
|
(74, 48),
(67, 55), (70, 76), (68, 81), (81, 80), (70, 75),
(65, 44), (63, 45), (64, 52), (78, 82), (71, 73),
(96, 96), (72, 70), (53, 43), (66, 62), (59, 44),
(78, 77), (57, 54), (72, 62), (71, 65) |
表Cのデータに対して国語と英語の相関係数を調べるStanスクリプトをリストC.1のように用意した。
リストC.1 相関係数を調べるStanスクリプト（AnalBiNormal.stan）
data {
int N;
vector[2] Data[N];
}
parameters {
real muA;
real muB;
real<lower = 0.0>
sgmA;
real<lower = 0.0>
sgmB;
real<lower = -1.0, upper
= 1.0> rho;
}
transformed parameters {
vector[2] mu;
cov_matrix[2] sgm;
mu[1] = muA;
mu[2] = muB;
sgm[1][1] = sgmA * sgmA;
sgm[1][2] = rho * sgmA *
sgmB;
sgm[2][1] = sgm[1][2];
sgm[2][2] = sgmB * sgmB;
}
model {
muA ~ normal(0.0, 1000.0);
muB ~ normal(0.0, 1000.0);
sgmA ~ exponential(0.001);
sgmB ~ exponential(0.001);
rho ~ uniform(-1.0, 1.0);
for (i in 1:N)
Data[i] ~ multi_normal(mu, sgm);
}
リストC.1のStanスクリプトを用いて表Cの相関係数のベイズ分析を行うPythonスクリプトをリストC.2のように用意した。
リストC.2 相関係数のベイズ分析を行うPythonスクリプト
import matplotlib.pyplot as plt
import pystan
import seaborn as sb
import numpy as np
import scipy.stats as ss
Data = [(74, 48), (67, 55), (70, 76),
(68, 81), (81, 80),
(70,
75), (65, 44), (63, 45), (64, 52), (78, 82),
(71,
73), (96, 96), (72, 70), (53, 43), (66, 62),
(59,
44), (78, 77), (57, 54), (72, 62), (71, 65)
]
StanData = {'N': len(Data), 'Data':
Data}
sm = pystan.StanModel(file =
'AnalBiNormal.stan')
fit = sm.sampling(data = StanData,
n_jobs = 1)
print(fit)
rho_mean = fit['rho'].mean()
rho_025P, rho_975P =
np.percentile(fit['rho'], [2.5, 97.5])
sb.kdeplot(fit['rho'])
plt.title('Posterior distribution of
rho' +
'\nMean = {0:.3f} 95%CI
= [{1:.3f}, {2:.3f}]'.format(
rho_mean, rho_025P, rho_975P), fontsize = 16)
plt.xlabel('rho', fontsize = 14)
plt.yticks([])
plt.show()
muA_mean = fit['muA'].mean()
muA_025P, muA_975P =
np.percentile(fit['muA'], [2.5, 97.5])
muB_mean = fit['muB'].mean()
muB_025P, muB_975P =
np.percentile(fit['muB'], [2.5, 97.5])
plt.figure(figsize = (10, 4))
plt.subplot(121)
plt.title('Posterior distribution of
muA' +
'\nmean = {0:.1f}
95%CI = [{1:.1f}, {2:.1f}]'.format(
muA_mean,
muA_025P, muA_975P))
sb.kdeplot(fit['muA'])
plt.yticks([])
plt.xlabel('muA')
plt.subplot(122)
plt.title('Posterior distribution of
muB' +
'\nmean = {0:.1f}
95%CI = [{1:.1f}, {2:.1f}]'.format(
muB_mean, muB_025P, muB_975P))
sb.kdeplot(fit['muB'])
plt.yticks([])
plt.xlabel('muB')
plt.show()
sgmA_mean = fit['sgmA'].mean()
sgmA_025P, sgmA_975P =
np.percentile(fit['sgmA'], [2.5, 97.5])
sgmB_mean = fit['sgmB'].mean()
sgmB_025P, sgmB_975P =
np.percentile(fit['sgmB'], [2.5, 97.5])
plt.figure(figsize = (10, 4))
plt.subplot(121)
plt.title('Posterior distribution of
sgmA' +
'\nmean = {0:.1f}
95%CI = [{1:.1f}, {2:.1f}]'.format(
sgmA_mean, sgmA_025P, sgmA_975P))
sb.kdeplot(fit['sgmA'])
plt.yticks([])
plt.xlabel('sgmA')
plt.subplot(122)
plt.title('Posterior distribution of
sgmB' +
'\nmean = {0:.1f}
95%CI = [{1:.1f}, {2:.1f}]'.format(
sgmB_mean, sgmB_025P, sgmB_975P))
sb.kdeplot(fit['sgmB'])
plt.yticks([])
plt.xlabel('sgmB')
plt.show()
tndncy = []
for ma, mb, sgm in zip(fit['muA'],
fit['muB'], fit['sgm']):
v1 =
ss.multivariate_normal.rvs([ma, mb], sgm)
v2 =
ss.multivariate_normal.rvs([ma, mb], sgm)
#print(v1, v2)
tndncy.append(1 if (v1[0] -
v2[0])*(v1[1] - v2[1]) > 0.0 else 0)
#print(tndncy[-1])
p_tndncy = np.array(tndncy).mean()
print('\nPosterior predictive
Prob(da*db>0) = {0:.3f}\n'.format(p_tndncy))
リストC.2のスクリプトを実行すると、Stanスクリプトのコンパイル後、MCMCサンプリングが始まる。MCMCサンプリングが終了すると、まずrhoの事後分布が表示される（図C.1）。
図C.1
rhoの事後分布の平均が0.757、95％確信区間が［0.513、0.907］であるので、かなり強い相関があるといえる。
図C.1のフォームを閉じると、図C.2のフォームが提示される。
図C.2
平均muAとmuBの事後分布が描かれている。
図C.2のフォームを閉じると、図C.3のフォームが提示される。
図C.3
標準偏差sgmAとsgmBの事後分布が描かれている。
図C.3のフォームを閉じると、事後予測の計算が始まる。
想定されている母集団から２人の学生をランダムに選んだとき、２人の間で１つの科目（例えば国語）の成績の優劣が他方の成績の優劣と同じ方向である確率を予測する。
1人の学生の成績をv1=(国語１、英語１)で、もう1人の学生の成績をv2=(国語２、英語２)で表すとき、２科目の成績の優劣が同じ方向であることは次式
で表される。この式の成り立つ確率をパラメータの事後分布に基づくサンプリングで推定する。
サンプリングに基づく推定が終了すると、端末に図C.4のように出力される。
図C.4
２つの科目の優劣の方向が一致するのは、４回に３回ぐらいの頻度であるという予想である。