繰り返し１要因分散分析（固定効果）
階層的ベイズ分析
2025.12：CmdStanPy対応に改訂
ベイズ分析のソフトStanのPython用インターフェースCmdStanPy対応に改訂した。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
繰り返し１要因固定効果のstrictly additive modelは次式で与えられる（Winer et al., 1991, p. 261）。
ここで、
である。
式（１）においては、パラメータ値は不定であるので、以下のようにパラメータを再設定する。
このとき、式（１）は以下のようになる。
式（２）、（４）より、
である。
式（５）、（６）のモデルをリスト１のStanスクリプトで表す。
リスト１ 繰り返し１要因固定効果モデルのStanスクリプト（OneFctrFixedRep.stan）
//
Yasuharu Okamoto, 2020.08, 2025.12
data {
int a;
int N;
array[N, a] real X;
}
parameters {
real<lower=1.0e-9>
sgm_s; // Hyper parameter
real mu0;
// Hyper parameter
real<lower=1.0e-9>
sgm0;
// Hyper parameter
real mu_t;
// Hyper parameter
real<lower=1.0e-9>
sgm_t; // Hyper parameter
real<lower = 1.0e-9>
sgm;
array[a - 1] real tau;
array[N] real mu;
}
transformed parameters {
array[a] real f_a;
real v = 0.0;
for (j in 2:a) {
v +=
tau[j - 1];
}
f_a[1] = -v / a;
for (j in 2:a) {
f_a[j] = tau[j - 1] + f_a[1];
}
}
model {
sgm_s - 1.0e-9 ~
exponential(0.001); // Hyper distribution
mu0 ~ normal(0.0,
1000.0);
// Hyper distribution
sgm0 - 1.0e-9 ~
exponential(0.001);
// Hyper distribution
mu_t ~ normal(0.0,
1000.0);
// Hyper distribution
sgm_t - 1.0e-9 ~
exponential(0.001); // Hyper distribution
sgm - 1.0e-9 ~
exponential(sgm_s);
for (i in 1:N) {
mu[i] ~ normal(mu0, sgm0);
}
for (j in 1:(a-1)){
tau[j] ~ normal(mu_t, sgm_t);
}
for (i in 1:N)
for
(j in 1:a) {
if (j == 1) {
X[i][j] ~ normal(mu[i], sgm);
}
else {
X[i][j] ~ normal(mu[i] + tau[j-1], sgm);
}
}
}
リスト１では、水準の効果tau[j]、および個人要因mu[i]は、それぞれデータ全体の状況を表す共通の分布からのサンプルであると考えて、ハイパー分布を設定した。スクリプトは以下のようである。
sgm_s - 1.0e-9 ~
exponential(0.001); // Hyper distribution
mu0 ~ normal(0.0,
1000.0);
// Hyper distribution
sgm0 - 1.0e-9 ~
exponential(0.001);
// Hyper distribution
mu_t ~ normal(0.0,
1000.0);
// Hyper distribution
sgm_t - 1.0e-9 ~
exponential(0.001); // Hyper distribution
sgm - 1.0e-9 ~
exponential(sgm_s);
for (i in 1:N) {
mu[i] ~ normal(mu0,
sgm0);
}
for (j in 1:(a-1)){
tau[j] ~ normal(mu_t, sgm_t);
}
事前分布
mu[i] ~
normal(mu0, sgm0);
および
tau[j] ~ normal(mu_t, sgm_t);
は、mu[i]およびtau[j]がrandom effectであると設定しているとも考えられる。すなわち、fixed effectモデルに階層的モデルを設定すると、random effectモデルと確率モデルとしては同じとなる。ベイズ分析においては、事前分布の設定によりrandom effectモデルとfixed effectモデルが統合されると考えられる。
リスト１のStanスクリプトを用いて繰り返し１要因固定効果モデルのベイズ分析を行うスクリプトを、本ウェブサイト後半にリスト２として示した。これらのスクリプトファイルなどは、圧縮ファイルfilesOneFctrFixedRep.zipとしてまとめた。ダウンロード展開すれば自由に利用できる。PythonにおけるStan（PyStan2）の使い方については＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞第10章～第12章で説明している。
リスト１のStanスクリプトファイル（OneFctrFixedRep.stan）、リスト２のPythonスクリプトファイル（OneFctrFixedRep.py）、および入力データファイル（図１のファイルの場合は、datarep.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移して、CmdStanPyのインストールされた環境において、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan)
*****/filesOneFctrFixedRep$ python OneFctrFixedRep.py
上のコマンドを実行すると、入力データファイル名、出力テキストファイル名、画像の保存ファイル名の設定が求められる。
(stan) *****/filesOneFctrFixedRep$
python OneFctrFixedRep.py
入力データファイル（*.xlsx）= datarep.xlsx
出力ファイル（*.txt）= Results.txt
画像保存ファイル（*.png）= MyFigure.png
入力データファイルは、図１のように用意する。
図１ 入力データ例（datarep.xlsx）
入力データファイルは、Excelファイルとして用意する。１行目に、変数名と要因の水準名を書き、データ値は２行目から設定する。1列目は、実験参加者など各データ組の識別値であるが、同じ値を複数回用いてもよい。文字列でも通し番号でもよい。２列目からそれぞれの水準のデータ値を書く。
入力データファイル名を設定して「Enter」キーを押すと、次に計算結果の出力ファイル名の設定が求められる。ファイル拡張子が「.txt」である任意のテキストファイル名を設定する。
出力ファイル名の設定後、「Enter」キーを押すと、グラフの画像保存ファイル名の設定が求められる。ファイル拡張子が「.png」である任意のファイル名を設定する。設定したファイル名の「.png」を画像に合わせて編集した名前でファイルが保存される。
保存画像ファイル名を設定して「Enter」キーを押すと、データが読み込まれ、Stanスクリプトのコンパイルが始まる。このコンパイルは、多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、終了するとトレース図が表示される（図２）。
図２
左側に事後分宇のグラフ、右側にMCMCサンプリングのトレース図が表示されている。
図2のWindowを閉じると、図3のグラフが表示される。
図３
ハイパーパラメータに対応するハイパー事前分布のグラフである。
左側のグラフは、mu[i]に対して設定される分布normal(mu0, sgm0)である。
mu[i] ~ normal(mu0, sgm0);
右側のグラフは、tau[j]に対して設定される分布normal(mu_t, sgm_t)である。
tau[j] ~ normal(mu_t, sgm_t);
それぞれパラメータ(mu0とsgm0、およびmu_tとsgm_tをMCMCサンプルからランダムに100対を選んで描画している。赤の太い曲線はそれらの平均値である。スクリプトを以下に示す。
rng = np.random.default_rng()
idx_slct = rng.choice(len(fit['mu0']),
size=100, replace=False)
plt.figure(figsize=(12,5))
plt.subplot(121)
xcoord = np.linspace(min_X, max_X,
1000)
curve_list = []
for idx in idx_slct:
ycoord = ss.norm.pdf(xcoord,
loc=fit['mu0'][idx], scale=fit['sgm0'][idx])
plt.plot(xcoord,ycoord,
alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r',
label='Mean values')
plt.legend()
plt.title('正規分布（mu0、sgm0）', fontsize=16)
plt.subplot(122)
mu_t_coord = np.linspace(min_X -
max_X, max_X - min_X, 1000)
curve_list = []
for idx in idx_slct:
ycoord =
ss.norm.pdf(mu_t_coord,
loc=fit['mu_t'][idx], scale=fit['sgm_t'][idx])
plt.plot(mu_t_coord,ycoord,
alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(mu_t_coord, mn_curve, c='r',
label='Mean values')
plt.legend()
plt.title('正規分布（mu_t、sgm_t）', fontsize=16)
plt.legend()
plt.tight_layout()
plt.savefig(foutnm[:-4] +
'_h_distri.png')
savedFiles.append(foutnm[:-4] +
'_h_distri.png')
plt.show()
図３のWindowを閉じると、図４のグラフが表示される。
図４
各水準の値
の事後分布のグラフである。ラベルに示されている値は、事後分布の中央値である。Gelman et al. (2021)は、事後分布の代表値として、平均値より安定した中央値を勧めている。
A2、A3、A４の分布は互いに離れているが、A４の分布は、原点０を中央に含んでいるので、A1＝０と重なっていると解釈できる。
図4のWindowを閉じると、図５のグラフが表示される。
図５
各水準の値
の事後分布のグラフである。はf_aで表されている。
との関係は、式（４）で表されている。ただし、Stanスクリプト中では、添字は１から始まるという制約があることに注意。A３、A１とA４、A2の３グループに事後分布が分かれていると解釈できる。
図５のWindowを閉じると、図６のグラフが表示される。
図６
効果量の事後分布のグラフである。分布の代表値として中央値が表示されている。
効果量は次式で与えられている（cf. Kirk, 1995, p. 181）。
スクリプトは、以下の通りである。
effcts = []
for a_s, sd
in zip(fit_f_a.T, fit['sgm']):
v = np.std(a_s)
effcts.append(v / sd)
図６のWindowを閉じると、図７のグラフが表示される。
図７
データの事後予測値のグラフである。スクリプトを以下に示す。
for j in range(a):
v_smpls = []
for v_mus, v_taus, v_sgm
in \
zip(fit_mu.T, fit_tau.T, fit['sgm']):
v_list = []
for
i in range(N):
if j == 0:
v = ss.norm.rvs(v_mus[i], v_sgm)
else:
v = ss.norm.rvs(v_mus[i] + v_taus[j-1], v_sgm)
v_list.append(v)
v_smpls.append(np.mean(v_list))
sb.kdeplot(v_smpls,
label=A_name[j])
図７のWindowを閉じると、スクリプトの実行終了である。
出力ファイルの内容は以下のようになっている。
入力データファイル: datarep.xlsx
A = ['A1' 'A2' 'A3' 'A4']
a = 4
N = 10
[[182. 199. 185. 183.]
[186. 212. 184. 180.]
[222. 233. 207. 217.]
[207. 209. 186. 207.]
[187. 203. 181. 192.]
[235. 245. 227. 227.]
[197. 206. 189. 190.]
[205. 228. 197. 202.]
[195. 198. 181. 196.]
[204. 213. 180. 214.]]
Mean
MCSE
StdDev
MAD ... ESS_bulk ESS_tail ESS_bulk/s R_hat
lp__ -133.415000 0.113106 3.868120 3.732200 ...
1233.610 1960.970 2390.72 1.002060
sgm_s 0.332303 0.003368 0.243500 0.203627 ...
4126.020 1930.980 7996.17 0.999697
mu0 202.246000 0.123114 5.907820 5.344910 ...
2419.390 2060.300 4688.73 1.001780
sgm0 16.986800 0.101236 5.156860 4.371600 ...
3192.440 2346.430 6186.91 1.003940
mu_t -1.177490 1.097950 31.684300 12.795700 ...
1617.880 981.689 3135.42 1.001350
sgm_t 39.256100 1.919560 59.694000 15.547100 ...
1463.690 1672.600 2836.61 1.001950
sgm
6.216030 0.018794 0.859566 0.824404 ...
2205.900 2716.390 4275.00 1.000150
tau[1] 12.351400 0.076368 2.793390 2.674530 ...
1365.500 2291.330 2646.32 0.999924
tau[2] -10.177700 0.077105 2.770560 2.661170 ...
1324.600 2083.940 2567.06 1.000270
tau[3] -1.255060 0.073326 2.757750 2.697540 ...
1439.080 2072.520 2788.91 0.999909
mu[1] 187.599000 0.070008 3.527060 3.410420 ...
2581.670 2647.930 5003.24 1.000040
mu[2] 190.690000 0.076783 3.537180 3.494460 ...
2153.070 2530.950 4172.63 1.000190
mu[3] 218.830000 0.079997 3.502970 3.513950 ...
1973.880 2497.200 3825.34 1.000370
mu[4] 202.012000 0.075066 3.581150 3.344790 ...
2353.560 2349.330 4561.16 1.000040
mu[5] 190.981000 0.075189 3.463790 3.307510 ...
2154.840 2472.380 4176.04 1.000600
mu[6] 231.984000 0.079996 3.625130 3.522170 ...
2147.600 2209.670 4162.01 1.000070
mu[7] 195.503000 0.074129 3.481360 3.411680 ...
2273.750 2243.440 4406.49 1.000560
mu[8] 207.510000 0.079818 3.606390 3.502920 ...
2072.290 2562.620 4016.06 1.001100
mu[9] 192.665000 0.074323 3.482140 3.431040 ...
2237.350 2252.650 4335.95 1.000060
mu[10] 202.470000 0.078172 3.482650 3.378570 ...
2025.580 2737.000 3925.54 0.999797
f_a[1] -0.229651 0.055014 1.709910 1.665250 ... 993.141 1555.820 1924.69 1.000040
f_a[2] 12.121700 0.026550 1.711200 1.634020 ...
4223.680 3044.020 8185.43 1.000730
f_a[3] -10.407400 0.026968 1.679480 1.646620 ...
3969.940 2664.260 7693.68 1.000470
f_a[4] -1.484710 0.024706 1.648160 1.613000 ...
4538.860 2871.480 8796.24 1.000340
v
0.918605 0.220056 6.839640 6.661000 ... 993.141 1555.820 1924.69 1.000040
[25 rows x 11 columns]
tau-水準A1:
0に固定
tau-水準A2:
平均値 = 12.351
中央値 = 12.352
Q1 = 10.534 Q3 = 14.148
95% CI = [6.870, 17.854]
tau-水準A3:
平均値 = -10.178
中央値 = -10.189
Q1 = -11.990 Q3 = -8.399
95% CI = [-15.496, -4.586]
tau-水準A4:
平均値 = -1.255
中央値 = -1.288
Q1 = -3.046 Q3 = 0.572
95% CI = [-6.740, 4.261]
f_a水準A1:
平均値 = -0.230
中央値 = -0.190
Q1 = -1.352 Q3 = 0.893
95% CI = [-3.676, 3.109]
f_a水準A2:
平均値 = 12.122
中央値 = 12.127
Q1 = 11.029 Q3 = 13.233
95% CI = [8.714, 15.492]
f_a水準A3:
平均値 = -10.407
中央値 = -10.409
Q1 = -11.520 Q3 = -9.299
95% CI = [-13.654, -7.091]
f_a水準A4:
平均値 = -1.485
中央値 = -1.469
Q1 = -2.578 Q3 = -0.410
95% CI = [-4.824, 1.737]
効果量:
平均値 = 1.333
中央値 = 1.330
Q1 = 1.170 Q3 = 1.491
95% CI = [0.875, 1.800]
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
Kirk, R.E. (1995). Experimental design: Procedures for the behavioral sciences, third ed. Brooks/Cole Publishing Company.
Winer, B. J., Brown, D.R., & Michels, K.M. (1991). Statistical principles in experimental design, third ed. McGraw-Hill, Inc.
リスト２ 繰り返し１要因固定効果モデルのPythonスクリプト例（OneFctrFixedRep.py）
import numpy as np
import scipy.stats as ss
import scipy.optimize as so
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
from cmdstanpy import CmdStanModel
import seaborn as sb
import arviz as az
finnm = input('入力データファイル（*.xlsx）= ')
txtoutnm = input('出力ファイル（*.txt）= ')
foutnm = input('画像保存ファイル（*.png）= ')
txtout = open(txtoutnm, 'w')
#
# Excelファイルの読み込み
# モジュール xlrd がインストール済みであること
#
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)
txtout.write('\n入力データファイル: {}\n'.format(finnm))
A_name = data[0, 1:]
a = len(A_name)
N = len(data) - 1
X = np.empty((N, a), dtype = float)
print('A = ', A_name)
print('a = ', a)
print('N = ', N)
txtout.write('\nA =
{}\n'.format(A_name))
txtout.write('a = {}\n'.format(a))
txtout.write('N = {}\n'.format(N))
for i in range(N):
for j in range(a):
X[i][j] = data[i + 1][j + 1]
print(X)
txtout.write('\n{}\n'.format(X))
min_X = np.min(X)
max_X = np.max(X)
StanData = {'a': a, 'N': N, 'X': X}
model =
CmdStanModel(stan_file='OneFctrFixedRep.stan')
fit = model.sample(data=StanData,
adapt_delta=0.99)
print(fit.diagnose())
print(fit.summary())
txtout.write(f'{fit.summary()}')
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
az.plot_trace(infdata, figsize=(14,
9))
plt.tight_layout()
plt.savefig(foutnm[:-4] +
'_trace.png')
savedFiles = [foutnm[:-4] +
'_trace.png']
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型への変換
print(fit.keys())
rng = np.random.default_rng()
idx_slct = rng.choice(len(fit['mu0']),
size=100, replace=False)
plt.figure(figsize=(12,5))
plt.subplot(121)
xcoord = np.linspace(min_X, max_X,
1000)
curve_list = []
for idx in idx_slct:
ycoord = ss.norm.pdf(xcoord,
loc=fit['mu0'][idx], scale=fit['sgm0'][idx])
plt.plot(xcoord,ycoord,
alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r',
label='Mean values')
plt.legend()
plt.title('正規分布（mu0、sgm0）', fontsize=16)
plt.subplot(122)
mu_t_coord = np.linspace(min_X -
max_X, max_X - min_X, 1000)
curve_list = []
for idx in idx_slct:
ycoord =
ss.norm.pdf(mu_t_coord,
loc=fit['mu_t'][idx], scale=fit['sgm_t'][idx])
plt.plot(mu_t_coord,ycoord,
alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(mu_t_coord, mn_curve, c='r',
label='Mean values')
plt.legend()
plt.title('正規分布（mu_t、sgm_t）', fontsize=16)
plt.legend()
plt.tight_layout()
plt.savefig(foutnm[:-4] +
'_h_distri.png')
savedFiles.append(foutnm[:-4] +
'_h_distri.png')
plt.show()
fit_tau = []
for j in range(a-1):
fit_tau.append(fit[f'tau[{j+1}]'])
fit_tau = np.array(fit_tau)
print('fit_tau: ', np.shape(fit_tau))
v_meds = []
for j in range(0, a):
if j == 0:
print('\ntau-{}:'.format(A_name[j]))
print(' 0に固定')
txtout.write('\n\ntau-水準{}:\n'.format(A_name[j]))
txtout.write(' 0に固定\n')
else:
print('\ntau-{}:'.format(A_name[j]))
txtout.write('\n\ntau-水準{}:\n'.format(A_name[j]))
v_mean = fit_tau[j-1].mean()
v_025p, v_Q1, v_med, v_Q3, v_975p = \
np.percentile(fit_tau[j-1], [2.5, 25, 50, 75, 97.5])
v_meds.append(v_med)
print(' 平均値 = {0:.3f}'.format(v_mean))
print(' 中央値 = {0:.3f}'.format(v_med))
print(' Q1 =
{0:.3f} Q3 =
{1:.3f}'.format(v_Q1, v_Q3))
print(' 95% CI =
[{0:.3f}, {1:.3f}]'.format(v_025p, v_975p))
txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))
txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))
txtout.write('
Q1 = {0:.3f} Q3
= {1:.3f}\n'.format(v_Q1, v_Q3))
txtout.write('
95% CI = [{0:.3f}, {1:.3f}]\n'.format(v_025p, v_975p))
plt.figure(figsize = (8,5))
plt.title('水準の値 $\\tau$' + \
'\n$\\tau$({}) = 0, 固定'.format(A_name[0]), fontsize = 20)
for j in range(1, a):
sb.kdeplot(fit_tau[j-1],
label=A_name[j] + f': {v_meds[j-1]:.3f}')
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.savefig(foutnm[:-4] + '_tau.png')
savedFiles = [foutnm[:-4] +
'_tau.png']
plt.show()
fit_f_a = []
for j in range(a):
fit_f_a.append(fit[f'f_a[{j+1}]'])
fit_f_a = np.array(fit_f_a)
print('fit_f_a: ', np.shape(fit_f_a))
v_meds = []
for j in range(a):
print('\nf_a水準{}:'.format(A_name[j]))
txtout.write('\n\nf_a水準{}:\n'.format(A_name[j]))
v_mean =
fit_f_a[j].mean()
v_025p, v_Q1, v_med, v_Q3,
v_975p = \
np.percentile(fit_f_a[j], [2.5, 25, 50, 75, 97.5])
v_meds.append(v_med)
print(' 平均値 = {0:.3f}'.format(v_mean))
print(' 中央値 = {0:.3f}'.format(v_med))
print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_Q1,
v_Q3))
print(' 95% CI = [{0:.3f},
{1:.3f}]'.format(v_025p, v_975p))
txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))
txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))
txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_Q1,
v_Q3))
txtout.write(' 95% CI = [{0:.3f},
{1:.3f}]\n'.format(v_025p, v_975p))
plt.figure(figsize = (8,5))
plt.title('水準の値 f_a', fontsize = 20)
for j in range(a):
sb.kdeplot(fit_f_a[j],
label=A_name[j]+f': {v_meds[j]:.3f}')
plt.yticks([])
plt.legend()
plt.savefig(foutnm[:-4] + '_f_a.png')
savedFiles.append(foutnm[:-4] +
'_f_a.png')
plt.show()
effcts = []
for a_s, sd in zip(fit_f_a.T,
fit['sgm']):
v = np.std(a_s)
effcts.append(v / sd)
effcts = np.array(effcts)
print('\n効果量:')
txtout.write('\n\n効果量:\n')
v_mean = effcts.mean()
v_025p, v_Q1, v_med, v_Q3, v_975p =
np.percentile(effcts, [2.5, 25, 50, 75, 97.5])
print(' 平均値 = {0:.3f}'.format(v_mean))
print(' 中央値 = {0:.3f}'.format(v_med))
print(' Q1 = {0:.3f} Q3 = {1:.3f}'.format(v_Q1,
v_Q3))
print(' 95% CI = [{0:.3f},
{1:.3f}]'.format(v_025p, v_975p))
txtout.write(' 平均値 = {0:.3f}\n'.format(v_mean))
txtout.write(' 中央値 = {0:.3f}\n'.format(v_med))
txtout.write(' Q1 = {0:.3f} Q3 = {1:.3f}\n'.format(v_Q1,
v_Q3))
txtout.write(' 95% CI = [{0:.3f},
{1:.3f}]\n'.format(v_025p, v_975p))
plt.title('効果量' + f'Med.={v_med:.3f}' +
'\n' +
f'95%CI=[{v_025p:.3f}, {v_975p:.3f}]', fontsize = 20)
sb.kdeplot(effcts)
plt.ylabel([])
plt.tight_layout()
plt.savefig(foutnm[:-4] + '_es.png')
savedFiles.append(foutnm[:-4] +
'_es.png')
plt.show()
fit_mu = []
for i in range(N):
fit_mu.append(fit[f'mu[{i+1}]'])
fit_mu = np.array(fit_mu)
fit_tau = []
for j in range(a-1):
fit_tau.append(fit[f'tau[{j+1}]'])
fit_tau = np.array(fit_tau)
for j in range(a):
v_smpls = []
for v_mus, v_taus, v_sgm
in \
zip(fit_mu.T, fit_tau.T, fit['sgm']):
v_list = []
for
i in range(N):
if j == 0:
v = ss.norm.rvs(v_mus[i], v_sgm)
else:
v = ss.norm.rvs(v_mus[i] + v_taus[j-1], v_sgm)
v_list.append(v)
v_smpls.append(np.mean(v_list))
sb.kdeplot(v_smpls,
label=A_name[j])
plt.title('事後予測値', fontsize=16)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig(foutnm[:-4] + '_pred.png')
savedFiles.append(foutnm[:-4] +
'_pred.png')
plt.show()
txtout.close()
savedFiles.append(txtoutnm)
print('\n保存ファイル:')
for nm in savedFiles:
print(' ', nm)
print()