Stan、Python、Windows
平均値の差、および相関係数のベイズ分析
2025.11：CmdStanPy対応に改訂
StanのPython用インターフェースとしてCmdStanPyがある。CmdStanPyの解説を、このウェブサイトで行った。
なお、Stan (PyStan2)については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。
ここでは、以下において統計分析の基本的方法である対応のあるデータの平均値の差のベイズ分析、独立なデータの平均値の差のベイズ分析、および相関係数のベイズ分析のStanによるベイズ分析の例を示す。スクリプトファイルは圧縮ファイル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;
array[N] vector[2] dt;
}
transformed data {
array[N] real diff;
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スクリプト（AnalPair.py）
import matplotlib.pyplot
as plt
import numpy
as np
from cmdstanpy
import CmdStanModel
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}
model = CmdStanModel(stan_file='AnalPair.stan')
fit = model.sample(data=StanData)
print(fit.summary())
fit = fit.draws_pd() # Pandas DataFrame型へ変換
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()
リストA1のStanスクリプトファイル（AnalPair.stan）とリストA.2のPythonスクリプトファイル（AnalPair.py）を同じフォルダ内に置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境において、次のコマンドを実行する。
(stan) *****/StanFiles$ python AnalPair.py
上のコマンドを実行すると、Stanスクリプトのコンパイルが始まる。コンパイルには、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了すると、図A.1のグラフが表示される。
図A.1
図A.1には差の平均の事後分布が描かれている。平均値が9.2、95％確信区間が［7.3、11.0］であるので、健康法が血圧を下げる効果が認められると考えてよさそうである。
図A.1のWindowを閉じると、図A.2のグラフが表示される。
図A.2
差の分布の標準偏差の事後分布である。
図A.2のWindowを閉じると図A.3のグラフが表示される。
図A.3
個々の差の値は、平均値と標準偏差から確率変数の実現値として与えられるが、その分布が事後予測として示されている。事後予測分布における差が正である確率が約98.6％であるので、健康法を実行すればほぼ確実に効果が期待できると予想される。
図A.3のWindowを閉じると、スクリプトの実行終了である。
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スクリプト（AnalInd.stan）
data {
int NA;
int NB;
array[NA] real DataA;
array[NB] real DataB;
}
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スクリプト（AnalInd.py）
import matplotlib.pyplot
as plt
import numpy
as np
from cmdstanpy
import CmdStanModel
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}
model = CmdStanModel(stan_file='AnalInd.stan')
fit = model.sample(data=StanData)
print(fit.summary())
fit = fit.draws_pd() # Pandas DataFrame型への変換
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.1のStanスクリプトファイル（AnalInd.stan）とリストB.2のPythonスクリプトファイル（AnalInd.py）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされている環境で、次のコマンドを実行する。
(stan) *****/StanFiles$ python AnalInd.py
上のコマンドを実行すると、Stanスクリプトが読み込まれて、コンパイルが始まる。コンパイルには、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了すると、図B.1のグラフが表示される。
図B.1
平均の差mu_dの事後分布のグラフである。平均が−8.89、95％確信区間が［−11.63、−6.24］であり、mu_dは確信をもって負である、すなわちA市の平均値の方がB市の平均値より低いと言える。
図B.1のWindowを閉じると、図B.2のグラフが表示される。
図B.2
図B.2のフォームでは、A市とB市の平均値の事後分布が描かれている。２つの事後分布は明確に分離されているので、A市の平均値の方が低いと言えそうである。
図B.2のWindowを閉じると、図B.3のグラフが表示される。
図B.3
図B.3のグラフには、標準偏差の事後分布が描かれている。
図B.3のWindowを閉じると、図B.4のグラフが提示される。
図B.4
図B.4のグラフは、2つの市、A市とB市からランダムに一人ずつサンプルを採り、血圧の差を求めたときの事後予測分布である。A市の人の方がB市の人より血圧が高い確率が約0.059である、すなわちランダムに一人ずつ選んで血圧を比べたとき、A市の人の方が血圧が高いのは20回に１回ぐらいの頻度であると予想されている。
図B.4のWindowを閉じると、スクリプトの実行終了である。
学生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;
array[N] vector[2] Data;
}
parameters {
real muA;
real muB;
real<lower = 0.0,
upper=1.0> vsgmA;
real<lower = 0.0,
upper=1.0> vsgmB;
real<lower = 0.0, upper =
1.0> vrho;
}
transformed parameters {
real sgmA;
real sgmB;
real rho;
vector[2] mu;
cov_matrix[2]
sgm;
sgmA
= 0.001 + vsgmA*100.0;
sgmB
= 0.001 + vsgmB*100.0;
rho = 1.999*(vrho-0.5);
mu[1] = muA;
mu[2] = muB;
sgm[1,1]
= sgmA * sgmA;
sgm[1,2]
= rho * sgmA * sgmB;
sgm[2,1]
= sgm[1][2]; // rho * sgmB
* sgmA;
// sgm[1][2];
sgm[2,2]
= sgmB * sgmB;
}
model {
muA
~ normal(0.0, 1000.0);
muB
~ normal(0.0, 1000.0);
vsgmA
~ beta(1, 1);
vsgmB
~ beta(1, 1);
vrho
~ beta(1, 1);
for (i
in 1:N)
Data[i] ~ multi_normal(mu,
sgm);
}
リストC.1のStanスクリプトを用いて表Cの相関係数のベイズ分析を行うPythonスクリプトをリストC.2のように用意した。
リストC.2 相関係数のベイズ分析を行うPythonスクリプト（AnalBiNormal.py）
import matplotlib.pyplot
as plt
from cmdstanpy
import CmdStanModel
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}
model = CmdStanModel(stan_file='AnalBiNormal.stan')
fit = model.sample(data=StanData)
print(fit.summary())
fit = fit.draws_pd()
# Pandas DataFrame型への変換
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()
print(fit.keys())
tndncy = []
for ma, mb, sgm11, sgm12, sgm21, sgm22
in zip(fit['muA'], fit['muB'],
fit['sgm[1,1]'], fit['sgm[1,2]'],
fit['sgm[2,1]'], fit['sgm[2,2]']):
v1 = ss.multivariate_normal.rvs([ma,
mb],
[[sgm11,sgm12], [sgm21, sgm22]])
v2 = ss.multivariate_normal.rvs([ma,
mb],
[[sgm11, sgm12], [sgm21, sgm22]])
tndncy.append(1
if (v1[0] - v2[0])*(v1[1] - v2[1]) > 0.0 else 0)
p_tndncy = np.array(tndncy).mean()
print('\nPosterior
predictive Prob(da*db>0) = {0:.3f}\n'.format(p_tndncy))
リストC1のStanスクリプトファイル（AnalBiNormal.stan）とリストC2のPythonスクリプトファイル（AnalBiNormal.py）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境で、次のコマンドを実行する。
(stan) *****/StanFiles$ python AnalBiNormal.py
上のコマンドを実行すると、ファイルが読み込まれ、Stanスクリプトのコンパイルが始まる。コンパイルは、多少時間が掛かる。コンパイルが終わると、MCMCサンプリングが行われ、終了すると、図C.1のグラフが表示される。
図C.1
rhoの事後分布のグラフである。rhoの事後分布の平均が0.753、95％確信区間が［0.493、0.901］であるので、かなり強い相関があるといえる。
図C.1のWindowを閉じると、図C.2のグラフが表示される。
図C.2
平均muAとmuBの事後分布が描かれている。
図C.2のWindowを閉じると、図C.3のフォームが提示される。
図C.3
標準偏差sgmAとsgmBの事後分布が描かれている。
図C.3のフォームを閉じると、事後予測の計算が始まる。
想定されている母集団から２人の学生をランダムに選んだとき、２人の間で１つの科目（例えば国語）の成績の優劣が他方の成績の優劣と同じ方向である確率を予測する。
1人の学生の成績をv1=(国語１、英語１)で、もう1人の学生の成績をv2=(国語２、英語２)で表すとき、２科目の成績の優劣が同じ方向であることは次式
で表される。この式の成り立つ確率をパラメータの事後分布に基づくサンプリングで推定する。
サンプリングに基づく推定が終了すると、端末に図C.4のように出力される。
図C.4
２つの科目の優劣の方向が一致するのは、４回に３回ぐらいの頻度であるという予想である。