２要因クロスデータの３パラメータ対数正規分布によるベイズ分析
主効果、相互作用の平均値、標準偏差、モード、シフトにおける分析
岡本安晴
2026.07
２要因クロスデザインの標準分析法は分散分析であるが、分散分析は、分散（平方和）の分割に基づく分析法である。分散の分割に基づいて要因の効果の有無が調べられる。しかし、要因の効果は、要因の各水準の値が示されると分かり易い。各水準と従属変数の関係を確率モデルで表して、水準の値の事後分布を求めてみた。
分散の分割の場合、（等分散）正規分布に基づいて分析が行われる。これに対して、設定した確率モデルのパラメータにより要因の効果を調べる場合は、確率モデルは正規分布以外でもよい。確率モデルとして３パラメータ対数正規分布を設定してみた。
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は３パラメータ対数正規分布に従うという（Johnson et al., 1994）。確率密度関数は次式（１）のように書くことができる。
このとき、平均値、分散、モードは次式で与えられる。
いま、要因Aの第
水準、要因Bの第
水準における
番目のデータ
が、３パラメータ対数正規分布
に従うものとする。
は、式（１）の３パラメータ対数正規分布を表す。要因Aの水準数は
、要因Ｂの水準数は
とする。要因Aの水準が水準で要因Bの水準が水準である条件におけるデータ数
は、条件ごとに異なっていてもよい。
各パラメータにおける主効果、交互作用を標準の分散分析における考え方に従って、以下のように求める。
平均
パラメータとしての平均
は、式(2)により与えられる。
まず、全体の平均を
とおく。
第１要因（要因A）の主効果を
とおき、第２要因（要因B）の主効果も同様に
とおく。
このとき、交互作用を
とおく。よって、
である。
標準偏差、モード、シフト
標準偏差（式（３）の平方根）、モード（式（４））、シフト（モデル（１）におけるパラメータ
）についても、主効果、交互作用を同様に与える。
上のモデルに基づくStanスクリプトをリスト１のように作成した。リスト１のStanスクリプトを用いて２要因クロスデータのベイズ分析を行うPythonスクリプトをリスト２のように作成した。ファイルは、TwoFac3PLNormfiles.zipにまとめた。
入力データファイルは、図１に示す形式で用意する。
・
・
・
図１
図１のファイルは、第１要因Factor-AがA1、A2、A3の３水準、第２要因Factor-BがB1、B2の２水準で、各条件におけるデータ数は空白セルがあり同じではない。第１列目の１行目には、要因Aの水準数（半角数字）を書き、数字の後に半角スラッシュを置いて要因名を書く。第１列目の2行目には、要因Bの水準数（半角数字）を書き、数字の後に半角スラッシュ/を置いて要因名を書く。第２列目以降は、第１行目には要因Aの水準名を書き、第２行目に要因Bの水準名を書く。要因Bの水準は要因Aの水準の下に書き、水準Bの順序は、要因Aの各水準の下で同じ順序とする。第１列目の３行目以降は、データの識別値を書く。
データ値は第３行目以降に、第１列目のデータの識別値に続けて、第２列目以降に該当する要因Aの水準と要因Bの水準の組み合わせの下に書く。なお、データの識別値は、図２では文字列であるが、これは通し番号でもよい。また、同じ識別値が複数回用いられてもよい。なお、データ値のセルが空白であれば欠損値として扱われる。
入力データファイルとリスト１およびリスト２のスクリプトファイルを同じフォルダ内に置く。そのフォルダにカレントディレクトリを移して、次のpythonコマンドを実行する。
(stan) ****/TwoFac3PLNormfiles$ python
ba2Fctrs3PLNormal.py
Input data file(*.xlsx) = data.xlsx
入力データファイル名の設定が求められる。上の例では、図１のファイル名data.xlsxが設定されている。
なお、現在（2026.07.20）、matplotlib 3.11ではエラーが出る。matplotlib 3.10では、大丈夫であった。matplotlib 3.10への交換は、次のコマンド
conda install
matplotlib=3.10
で、出来た。
図１の形式のデータに欠損値があれば、次のスクリプトで処理され、欠損値のセルは詰められる。
ipos = 0
for j in range(J):
for k in range(K):
v =
X[:,ipos]
print(v[:5])
v =
np.array(v)
v =
v[np.isnan(v) == False]
Ns[j][k] = len(v)
for
i in range(len(v)):
X3dim[j][k][i] = v[i]
if
len(v) < N:
for i in range(len(v), N):
X3dim[j][k][i] = -1
ipos
+= 1
ファイルが読み込まれると、Stanスクリプトがコンパイルされ、MCMCサンプリングが始まる。
MCMCサンプリングが終了すると、トレース図が表示される（図２）。
図２
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータmean（式（2））、標準偏差sd（式（3）の平方根）、モードmode（式（4））、theta（グラフではと表記）の各条件における事後分布のグラフである。
図３のWindowを閉じると、図４のグラフが表示される。
図４
パラメータmeanの要因Aおよび要因Bにおける各条件の主効果の事後分布である。要因Aにおいては、水準a1と水準a3の事後分布が重なっており、水準a2の事後分布がそれらとは離れて右側に位置している。すなわち、（水準a1≒水準A3）＜水準a2である。要因Bにおいては、２つの水準b1とb2の主効果の事後分布は離れている。すなわち、水準b1＜水準b2である。
図４のWindowを閉じると、図５のグラフが表示される。
図５
パラメータmeanにおける交互作用のグラフである。要因Aの各水準
と要因Bの各水準
との組み合わせ
における交互作用の事後分布である。破線を０の位置に引いて、交互作用の事後分布を箱図で表している。下の段には、水準水準におけるデータの平均値と主効果のみで表される値
のグラフが描かれている。データの平均値と主効果のみで表される値の差が上段の交互作用のグラフに反映されている。
なお、対数変換したデータの平均値と元のデータの平均値は、一方では同じ平均値であっても、他方では等しくない平均値であることがある。この例を示したウェブサイトも用意した。本ウェブサイトの平均値meanは、変換を行わない元データの平均値である。
図５のWindowを閉じると、図６のグラフが表示される。
図６
標準偏差パラメータ
の要因Aおよび要因Bの各水準に対する主効果の事後分布のグラフである。要因Aにおいては（左のグラフ）、水準a1の右に水準a3のグラフがあり、その間に水準a2が水準a1寄りに位置している。水準a1＝＜水準a2＜水準a3と言えそうである。
図６のWindowを閉じると、図７のグラフが表示される。
図７
標準偏差パラメータsdにおける交互作用のグラフである。上段に交互作用の事後分布の箱図が描かれている。下段には、水準値ajとbkの組み合わせに他するデータの標準偏差（破線）と主効果のみから算出される値（標準偏差の予測値）（実線）のグラフが描かれている。データから算出された標準偏差の値と主効果のみからの予測値の差が、上段の交互作用の箱図に対応している。
図７のWindowを閉じると、図８のグラフが表示される。
図８
パラメータモードmodeにおける主効果の事後分布である。要因Aにおいては、水準a1とa3に対して水準a2の主効果が右側に離れている。水準a1とa3はともにa2から離れてまとまっているが、事後分布はお互いに異なっている。水準３＜水準１＜＜水準２となっている。
右側のグラフは要因Bの主効果のグラフである。水準２の主効果のグラフは水準１より右に離れて位置しており、水準１＜水準２となっている。
図８のWindowを閉じると、図９のグラフが表示される。
図９
パラメータmodeにおける交互作用のグラフである。上段のグラフは、要因Aの水準j（aj）と要因Bの水準k（bk）における交互作用の事後分布を箱図で表している。
下段のグラフは、ajとbkの組み合わせにおけるmodeの推定値（事後分布の中央値、モードはデータから直接算出できない）（破線）と主効果のみによる予測値（実線）のグラフである。パラメータmodeの推定値と主効果のみによる予測値の差が、上段の交互作用のグラフに反映されている。
図９のWindowを閉じると、図10のグラフが表示される。
図10
シフト値theta（）における主効果のグラフである。左のグラフは、要因Aの水準１（a1）、水準２（a2）、水準３（a3）の主効果の事後分布である。水準１と水準３の事後分布はほぼ重なっているが、水準２の事後分布はそれらに比べると右に寄っている。右側のグラフは、要因Bにおける主効果のグラフである。２つの水準のグラフは離れていて、水準２のグラフが水準１のグラフの右側にある。
図10のWindowを閉じると、図11のグラフが表示される。
図11
シフトパラメータtheta（）における交互作用のグラフである。上段は、各交互作用の事後分布を箱図で表したものである。上段の箱図の分布を見ると、箱（中央の50％区間）が0（緑の破線）から離れており、交互作用は統計学的にすべて認められると言える。しかし、下段のthetaの推定値（事後分布の中央値）のグラフ（破線）と、主効果のみによる推定値（実線）との関係を見ると、主効果による変動量（b1＜b2）と比べると交互作用の大きさは小さく、交互作用は実質的には重要でない可能性がある。
図11のWindowを閉じると、スクリプトの実行終了である。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B.,
Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed.
CRC Press.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
Johnson, N. L., Kotz, S., & Barakrishnan, N. (1994) Continuous univariate distributions, Vol.1, 2nd Ed. John Wiley & Sons, Inc.
Kirk, R.E. (1995). Experimental design: Procedures for the behavioral sciences, third ed. Brooks/Cole Publishing Company.
Myers, J. L. & Well, A. D. (2003). Research design and statistical analysis, second ed. Lawrence Erlbaum Associates, Publishers.
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Winer, B. J., Brown, D.R., & Michels, K.M. (1991). Statistical principles in experimental design, third ed. McGraw-Hill, Inc.
リスト１ ２要因デザインデータ分析のStanスクリプト（twofctrLNorn.stan）
data {
int J;
int K;
int N;
array[J, K] int Ns;
array[J, K, N] real X;
}
transformed data {
array[J, K] real min_x;
for (j in 1:J) {
for (k
in 1:K) {
min_x[j, k] = X[j][k][1];
for (i in 2:Ns[j, k]) {
if (min_x[j,k] > X[j][k][i]) {
min_x[j,k] = X[j][k][i];
}
}
}
}
}
parameters {
array[J,K] real<lower=0,
upper=1> vmu;
array[J,K] real<lower=0,
upper=1> vsgm;
array[J,K] real<lower=0,
upper=1> vtheta;
}
transformed parameters {
array[J,K] real mu;
array[J,K] real sgm;
array[J,K] real theta;
for (j in 1:J) {
for (k in 1:K) {
mu[j,k] = -2000 + vmu[j,k]*4000;
sgm[j,k] = 0.0001 + vsgm[j,k] * 1000;
theta[j,k] = vtheta[j,k] * min_x[j,k];
}
}
}
model {
for (j in 1:J) {
vmu[j] ~ beta(1, 1);
vsgm[j] ~ beta(1, 1);
vtheta[j] ~ beta(1, 1);
}
for (j in 1:J) {
for
(k in 1:K) {
for (i in 1:Ns[j][k]) {
X[j][k][i] - theta[j][k] ~ lognormal(mu[j][k], sgm[j][k]);
}
}
}
}
generated quantities {
array[J, K] real mean;
array[J, K] real sd;
array[J, K] real mode;
for (j in 1:J) {
for
(k in 1:K) {
mean[j,k] = theta[j,k] + exp(mu[j,k] + 0.5*(sgm[j,k]^2));
sd[j,k] = (exp(2*mu[j,k] + (sgm[j,k]^2)) *
(exp(sgm[j,k]^2) - 1.0))^0.5;
mode[j,k] = theta[j,k] + exp(mu[j,k] - (sgm[j,k]^2));
}
}
/************************************************
mean
*************************************************/
real mean_pp;
// The global mean
array[J] real mean_jp; // Means at level aj
array[K] real mean_pk; // Means at level bk
array[J] real a_mean;
// Main effects
array[K] real b_mean;
// Main effects
array[J, K] real
ab_mean; // Interactions
real v;
mean_pp = 0.0;
for (j in 1:J) {
for
(k in 1:K) {
mean_pp += mean[j][k];
}
}
mean_pp /= J * K;
for (j in 1:J) {
v =
0.0;
for
(k in 1:K) {
v += mean[j][k];
}
mean_jp[j] = v / K;
a_mean[j] = mean_jp[j] - mean_pp;
}
for (k in 1:K) {
v =
0.0;
for
(j in 1:J) {
v += mean[j][k];
}
mean_pk[k] = v / J;
b_mean[k] = mean_pk[k] - mean_pp;
}
for (j in 1:J) {
for
(k in 1:K) {
ab_mean[j][k] = mean[j][k] - mean_jp[j] - mean_pk[k] + mean_pp;
}
}
/****************************************************
SD
*****************************************************/
real sd_pp;
// The global mean
array[J] real sd_jp; // Means at level aj
array[K] real sd_pk; // Means at level bk
array[J] real a_sd;
// Main effects
array[K] real b_sd;
// Main effects
array[J, K] real ab_sd; // Interactions
//real v;
sd_pp = 0.0;
for (j in 1:J) {
for (k
in 1:K) {
sd_pp += sd[j][k];
}
}
sd_pp /= J * K;
for (j in 1:J) {
v =
0.0;
for
(k in 1:K) {
v += sd[j][k];
}
sd_jp[j] = v / K;
a_sd[j] = sd_jp[j] - sd_pp;
}
for (k in 1:K) {
v =
0.0;
for
(j in 1:J) {
v += sd[j][k];
}
sd_pk[k] = v / J;
b_sd[k] = sd_pk[k] - sd_pp;
}
for (j in 1:J) {
for
(k in 1:K) {
ab_sd[j][k] = sd[j][k] - sd_jp[j] - sd_pk[k] + sd_pp;
}
}
/********************************************
Mode
**********************************************/
real mode_pp;
// The global mean
array[J] real mode_jp; // Means at level aj
array[K] real mode_pk; // Means at level bk
array[J] real a_mode;
// Main effects
array[K] real b_mode;
// Main effects
array[J, K] real
ab_mode; // Interactions
mode_pp = 0.0;
for (j in 1:J) {
for
(k in 1:K) {
mode_pp += mode[j][k];
}
}
mode_pp /= J * K;
for (j in 1:J) {
v =
0.0;
for
(k in 1:K) {
v += mode[j][k];
}
mode_jp[j] = v / K;
a_mode[j] = mode_jp[j] - mode_pp;
}
for (k in 1:K) {
v =
0.0;
for
(j in 1:J) {
v += mode[j][k];
}
mode_pk[k] = v / J;
b_mode[k] = mode_pk[k] - mode_pp;
}
for (j in 1:J) {
for (k
in 1:K) {
ab_mode[j][k] = mode[j][k] - mode_jp[j] - mode_pk[k] + mode_pp;
}
}
/***********************************************
Theta
************************************************/
real theta_pp;
// The global mean
array[J] real theta_jp; // Means at level aj
array[K] real theta_pk; // Means at level bk
array[J] real a_theta;
// Main effects
array[K] real b_theta;
// Main effects
array[J, K] real
ab_theta; // Interactions
theta_pp = 0.0;
for (j in 1:J) {
for
(k in 1:K) {
theta_pp += theta[j][k];
}
}
theta_pp /= J * K;
for (j in 1:J) {
v =
0.0;
for
(k in 1:K) {
v += theta[j][k];
}
theta_jp[j] = v / K;
a_theta[j] = theta_jp[j] - theta_pp;
}
for (k in 1:K) {
v =
0.0;
for
(j in 1:J) {
v += theta[j][k];
}
theta_pk[k] = v / J;
b_theta[k] = theta_pk[k] - theta_pp;
}
for (j in 1:J) {
for
(k in 1:K) {
ab_theta[j][k] = theta[j][k] - theta_jp[j] - theta_pk[k]
+ theta_pp;
}
}
}
リスト２ ２要因デザインデータ分析のPythonスクリプト（ba2Fctrs3PLNormal.py）
from cmdstanpy import CmdStanModel
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import arviz as az
data_pd =
pd.read_excel(input("Input data file(*.xlsx) = "), header=[0,1])
print(data_pd)
print(data_pd.keys())
print(data_pd.keys()[0])
J =
int(data_pd.keys()[0][0].split('/')[0])
print(J)
K =
int(data_pd.keys()[0][1].split('/')[0])
print(K)
X = np.array(data_pd.values[:,1:],
dtype=float)
print(X[:5], '\n', X[-5:])
N = len(X)
X3dim = np.empty((J,K,N))
Ns = np.empty((J,K))
ipos = 0
for j in range(J):
for k in range(K):
v =
X[:,ipos]
print(v[:5])
v =
np.array(v)
v =
v[np.isnan(v) == False]
Ns[j][k] = len(v)
for i
in range(len(v)):
X3dim[j][k][i] = v[i]
if
len(v) < N:
for i in range(len(v), N):
X3dim[j][k][i] = -1
ipos
+= 1
print(X3dim[:,:,:5], '\n',
X3dim[:,:,-5:])
Ns = np.array(Ns, dtype=int)
print('Ns_1 =\n', Ns)
Data = {'N':N, 'J':J, 'K':K, 'Ns':Ns, 'X':X3dim}
sm =
CmdStanModel(stan_file='twofctrLNorm.stan')
fit = sm.sample(data=Data,
max_treedepth=20)
print(fit.diagnose())
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# inf_dataを用いてのトレースプロット
az.plot_trace_dist(inf_data,
var_names=['mu','sgm','theta'])
plt.tight_layout()
plt.show()
fit_dframe = fit.draws_pd()
# PandasのDataFrame型に変換
print(fit_dframe.keys())
n_smpls = len(fit_dframe['mean[1,1]'])
means = np.empty((J,K,n_smpls))
sds = np.empty((J,K,n_smpls))
modes = np.empty((J,K,n_smpls))
thetas = np.empty((J,K,n_smpls))
for j in range(J):
for k in range(K):
means[j][k] = fit_dframe[f"mean[{j+1},{k+1}]"]
sds[j][k] = fit_dframe[f"sd[{j+1},{k+1}]"]
modes[j][k] = fit_dframe[f"mode[{j+1},{k+1}]"]
thetas[j][k] = fit_dframe[f"theta[{j+1},{k+1}]"]
plt.figure(figsize=(10,6))
plt.subplot(221)
plt.title('Posterior distributions of
mean', fontsize=16)
for j in range(J):
for k in range(K):
sb.kdeplot(means[j][k], label=f"a{j+1}b{k+1}")
plt.legend()
plt.xlabel('mean', fontsize=14)
plt.subplot(222)
plt.title('Posterior distributions of
sd', fontsize=16)
for j in range(J):
for k in range(K):
sb.kdeplot(sds[j][k], label=f"a{j+1}b{k+1}")
plt.legend()
plt.xlabel('sd', fontsize=14)
plt.subplot(223)
plt.title('Posterior distributions of
mode', fontsize=16)
for j in range(J):
for k in range(K):
sb.kdeplot(modes[j][k], label=f"a{j+1}b{k+1}")
plt.legend()
plt.xlabel('mode', fontsize=14)
plt.subplot(224)
plt.title('Posterior distributions of
$\\theta$', fontsize=16)
for j in range(J):
for k in range(K):
sb.kdeplot(thetas[j][k], label=f"a{j+1}b{k+1}")
plt.legend()
plt.xlabel('$\\theta$', fontsize=14)
plt.tight_layout()
plt.show()
##########################################
#
#
mean
#
#################################333333333
mean_pp =
fit_dframe["mean_pp"]
a_mean = np.empty((J, n_smpls))
b_mean = np.empty((K, n_smpls))
ab_mean = np.empty((J,K, n_smpls))
for j in range(J):
a_mean[j] =
fit_dframe[f"a_mean[{j+1}]"]
for k in range(K):
b_mean[k] =
fit_dframe[f"b_mean[{k+1}]"]
for j in range(J):
for k in range(K):
ab_mean[j][k] = fit_dframe[f"ab_mean[{j+1},{k+1}]"]
plt.figure(figsize=(10,3))
plt.subplot(121)
for j in range(J):
sb.kdeplot(a_mean[j], label=
f"a{j+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior distributions of
a/mean', fontsize=16)
plt.subplot(122)
for k in range(K):
sb.kdeplot(b_mean[k],
label=f"b{k+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior distributions of
b/mean', fontsize=16)
plt.tight_layout()
plt.show()
plt.figure(figsize=(10,6))
plt.subplot(2,2,(1,2))
ab_smpls = {}
for j in range(J):
for k in range(K):
ab_smpls[f"a{j+1}b{k+1}"] = ab_mean[j,k]
plt.plot([-1, J*K+1], [0, 0],
color='g', ls='--')
sb.boxplot(ab_smpls)
plt.xticks(fontsize=14)
plt.title('Posterior distributions of
interactions/mean', fontsize=16)
Xmeans = np.empty((J,K))
for j in range(J):
for k in range(K):
Xmeans[j][k] = np.mean(X3dim[j,k,:Ns[j,k]])
print(Xmeans)
print(Xmeans.shape)
mean_pp_hat = np.median(mean_pp)
a_hat = []
for j in range(J):
a_hat.append(np.median(a_mean[j]))
a_hat = np.array(a_hat)
print('a_hat =\n', a_hat)
b_hat = []
for k in range(K):
b_hat.append(np.median(b_mean[k]))
b_hat = np.array(b_hat)
plt.subplot(223)
for k in range(K):
plt.plot(range(J),
Xmeans[:,k], label=f"data in b{k+1}", ls='--')
for k in range(K):
plt.plot(range(J), a_hat +
b_hat[k] + mean_pp_hat,
label=f"main effects aj+b{k+1}")
plt.xticks(np.arange(0,J),
[f"a{j+1}" for j in range(J)], fontsize=14)
plt.legend()
plt.title('Main effects/mean',
fontsize=16)
plt.subplot(224)
for j in range(J):
plt.plot(range(K),
Xmeans[j], label=f"data in a{j+1}", ls='--')
for j in range(J):
plt.plot(range(K), b_hat +
a_hat[j] + mean_pp_hat,
label=f"main effects a{j+1}+bk")
plt.xticks(np.arange(0,K),
[f"b{k+1}" for k in range(K)], fontsize=14)
plt.legend()
plt.title('Main effects/mean',
fontsize=16)
plt.tight_layout()
plt.show()
##############################################
#
#
SD
#
##############################################
Xsds = np.empty((J,K))
for j in range(J):
for k in range(K):
Xsds[j][k] = np.std(X3dim[j,k,:Ns[j,k]])
print(Xsds)
print(Xsds.shape)
sd_pp = fit_dframe["sd_pp"]
a_sd = np.empty((J, n_smpls))
b_sd = np.empty((K, n_smpls))
ab_sd = np.empty((J,K, n_smpls))
for j in range(J):
a_sd[j] =
fit_dframe[f"a_sd[{j+1}]"]
for k in range(K):
b_sd[k] =
fit_dframe[f"b_sd[{k+1}]"]
for j in range(J):
for k in range(K):
ab_sd[j][k]
= fit_dframe[f"ab_sd[{j+1},{k+1}]"]
plt.figure(figsize=(10,3))
plt.subplot(121)
for j in range(J):
sb.kdeplot(a_sd[j], label=
f"a{j+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior
distributions/sd', fontsize=16)
plt.subplot(122)
for k in range(K):
sb.kdeplot(b_sd[k],
label=f"b{k+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior
distributions/sd', fontsize=16)
plt.show()
plt.figure(figsize=(10,6))
plt.subplot(2,2,(1,2))
ab_smpls = {}
for j in range(J):
for k in range(K):
ab_smpls[f"a{j+1}b{k+1}"] = ab_sd[j,k]
plt.plot([-1, J*K+1], [0, 0],
color='g', ls='--')
sb.boxplot(ab_smpls)
plt.xticks(fontsize=14)
plt.title('Posterior distributions of
interactions/sd', fontsize=16)
plt.tight_layout()
sd_pp_hat = np.median(sd_pp)
a_hat = []
for j in range(J):
a_hat.append(np.median(a_sd[j]))
a_hat = np.array(a_hat)
print('a_hat =\n', a_hat)
b_hat = []
for k in range(K):
b_hat.append(np.median(b_sd[k]))
b_hat = np.array(b_hat)
plt.subplot(223)
for k in range(K):
plt.plot(range(J),
Xsds[:,k], label=f"data in b{k+1}", ls='--')
for k in range(K):
plt.plot(range(J), a_hat +
b_hat[k] + sd_pp_hat,
label=f"main effects aj+b{k+1}")
plt.xticks(np.arange(0,J),
[f"a{j+1}" for j in range(J)], fontsize=14)
plt.legend()
plt.title('Main effects/sd',
fontsize=16)
plt.tight_layout()
plt.subplot(224)
for j in range(J):
plt.plot(range(K), Xsds[j],
label=f"data in a{j+1}", ls='--')
for j in range(J):
plt.plot(range(K), b_hat +
a_hat[j] + sd_pp_hat,
label=f"main effects a{j+1}+bk")
plt.xticks(np.arange(0,K),
[f"b{k+1}" for k in range(K)], fontsize=14)
plt.legend()
plt.title('Main effects/sd',
fontsize=16)
plt.tight_layout()
plt.show()
#################################################
#
#
Mode
#
#################################################
Xmodes = np.empty((J,K))
for j in range(J):
for k in range(K):
Xmodes[j][k] = np.median(fit_dframe[f"mode[{j+1},{k+1}]"])
print(Xmodes)
print(Xmodes.shape)
mode_pp =
fit_dframe["mode_pp"]
a_mode = np.empty((J, n_smpls))
b_mode = np.empty((K, n_smpls))
ab_mode = np.empty((J,K, n_smpls))
for j in range(J):
a_mode[j] =
fit_dframe[f"a_mode[{j+1}]"]
for k in range(K):
b_mode[k] =
fit_dframe[f"b_mode[{k+1}]"]
for j in range(J):
for k in range(K):
ab_mode[j][k] = fit_dframe[f"ab_mode[{j+1},{k+1}]"]
plt.figure(figsize=(10, 3))
plt.subplot(121)
for j in range(J):
sb.kdeplot(a_mode[j], label=
f"a{j+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior
distributions/mode', fontsize=16)
plt.subplot(122)
for k in range(K):
sb.kdeplot(b_mode[k],
label=f"b{k+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior
distributions/mode', fontsize=16)
plt.show()
plt.figure(figsize=(10, 6))
plt.subplot(2,2,(1,2))
ab_smpls = {}
for j in range(J):
for k in range(K):
ab_smpls[f"a{j+1}b{k+1}"] = ab_mode[j,k]
plt.plot([-1, J*K+1], [0, 0],
color='g', ls='--')
sb.boxplot(ab_smpls)
plt.xticks(fontsize=14)
plt.title('Posterior distributions of
interactions/mode', fontsize=16)
plt.tight_layout()
mode_pp_hat = np.median(mode_pp)
a_hat = []
for j in range(J):
a_hat.append(np.median(a_mode[j]))
a_hat = np.array(a_hat)
print('a_hat =\n', a_hat)
b_hat = []
for k in range(K):
b_hat.append(np.median(b_mode[k]))
b_hat = np.array(b_hat)
plt.subplot(223)
for k in range(K):
plt.plot(range(J),
Xmodes[:,k], label=f"est.in b{k+1}", ls='--')
for k in range(K):
plt.plot(range(J), a_hat +
b_hat[k] + mode_pp_hat,
label=f"main effects aj+b{k+1}")
plt.xticks(np.arange(0,J),
[f"a{j+1}" for j in range(J)], fontsize=14)
plt.legend()
plt.title('Main effects/mode',
fontsize=16)
plt.tight_layout()
plt.subplot(224)
for j in range(J):
plt.plot(range(K),
Xmodes[j], label=f"est.in a{j+1}", ls='--')
for j in range(J):
plt.plot(range(K), b_hat +
a_hat[j] + mode_pp_hat,
label=f"main effects a{j+1}+bk")
plt.xticks(np.arange(0,K),
[f"b{k+1}" for k in range(K)], fontsize=14)
plt.legend()
plt.title('Main effects/mode',
fontsize=16)
plt.tight_layout()
plt.show()
################################################
#
#
Theta
#
################################################
Xthetas = np.empty((J,K))
for j in range(J):
for k in range(K):
Xthetas[j][k] = np.median(fit_dframe[f"theta[{j+1},{k+1}]"])
print(Xthetas)
print(Xthetas.shape)
theta_pp =
fit_dframe["theta_pp"]
a_theta = np.empty((J, n_smpls))
b_theta = np.empty((K, n_smpls))
ab_theta = np.empty((J,K, n_smpls))
for j in range(J):
a_theta[j] =
fit_dframe[f"a_theta[{j+1}]"]
for k in range(K):
b_theta[k] =
fit_dframe[f"b_theta[{k+1}]"]
for j in range(J):
for k in range(K):
ab_theta[j][k] = fit_dframe[f"ab_theta[{j+1},{k+1}]"]
plt.figure(figsize=(10, 3))
plt.subplot(121)
for j in range(J):
sb.kdeplot(a_theta[j],
label= f"a{j+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior
distributions/theta', fontsize=16)
plt.subplot(122)
for k in range(K):
sb.kdeplot(b_theta[k],
label=f"b{k+1}")
plt.xlabel('')
plt.legend()
plt.title('Posterior
distributions/theta', fontsize=16)
plt.show()
plt.figure(figsize=(10, 6))
plt.subplot(2,2,(1,2))
ab_smpls = {}
for j in range(J):
for k in range(K):
ab_smpls[f"a{j+1}b{k+1}"] = ab_theta[j,k]
plt.plot([-1, J*K+1], [0, 0],
color='g', ls='--')
sb.boxplot(ab_smpls)
plt.xticks(fontsize=14)
plt.title('Posterior distributions of
interactions/theta', fontsize=16)
plt.tight_layout()
theta_pp_hat = np.median(theta_pp)
a_hat = []
for j in range(J):
a_hat.append(np.median(a_theta[j]))
a_hat = np.array(a_hat)
print('a_hat =\n', a_hat)
b_hat = []
for k in range(K):
b_hat.append(np.median(b_theta[k]))
b_hat = np.array(b_hat)
plt.subplot(223)
for k in range(K):
plt.plot(range(J),
Xthetas[:,k], label=f"est.in b{k+1}", ls='--')
for k in range(K):
plt.plot(range(J), a_hat +
b_hat[k] + theta_pp_hat,
label=f"main effects aj+b{k+1}")
plt.xticks(np.arange(0,J),
[f"a{j+1}" for j in range(J)], fontsize=14)
plt.legend()
plt.title('Main effects/theta',
fontsize=16)
plt.tight_layout()
plt.subplot(224)
for j in range(J):
plt.plot(range(K),
Xthetas[j], label=f"est.in a{j+1}", ls='--')
for j in range(J):
plt.plot(range(K), b_hat +
a_hat[j] + theta_pp_hat,
label=f"main effects a{j+1}+bk")
plt.xticks(np.arange(0,K),
[f"b{k+1}" for k in range(K)], fontsize=14)
plt.legend()
plt.title('Main effects/theta',
fontsize=16)
plt.tight_layout()
plt.show()