３パラメータ対数正規分布による１要因デザインデータのベイズ分析
岡本安晴
2026.07
分散分析は、分散（平方和）の分割に基づく分析法であるが、分析の目的は分散の分割に基づいて要因の効果の有無を調べることである。要因の効果は、要因の各水準の値が示されると分かり易い。各水準と従属変数の関係を確率モデルで表して、水準の値の事後分布を求めてみた。
分散の分割の場合、（等分散）正規分布に基づいて分析が行われるが、設定した確率モデルのパラメータにより要因の効果を調べる場合は、確率モデルは正規分布以外でもよい。確率モデルとして３パラメータ対数正規分布を設定してみた。
確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は３パラメータ対数正規分布に従うという（Johnson et al., 1994）。確率密度関数は次式（１）のように書くことができる。
このとき、平均値、分散、モードは次式で与えられる。
いま、要因Aの第
水準における
番目のデータ
が３パラメータ対数正規分布
に従うものとする。
は、式（１）の３パラメータ対数正規分布を表す。要因Aの水準数は
とする。水準におけるデータ数
は、水準ごとに異なっていてもよい。
上のモデルのStanスクリプトをリスト１のように用意した。
リスト１では、次のスクリプト
v = X[i][j] - theta[j];
v ~ lognormal(mu[j], sgm[j]);
で与えられたパラメータ値に対して、平均、標準偏差、モード（式(2)、(3)、(4)）を以下のように算出している。
mean[j] = theta[j] + exp(mu[j] + 0.5*(sgm[j]^2));
sd[j] = (exp(2*mu[j] + (sgm[j]^2))
* (exp(sgm[j]^2) - 1.0))^0.5;
mode[j] = theta[j] + exp(mu[j] - (sgm[j]^2));
リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプトをリスト２のように作成した。ファイルは、OneFacfiles.zipにまとめた。
入力データファイルは、図１のように用意する。
・
・
・
図１
1行目に水準名を書く。
各水準のデータ値を2行目から書く。水準ごとにデータ数が異なっていてもよい。欠損値があれば、空白にしておく。
リスト１とリスト２のスクリプトファイル、および入力データファイルを同じフォルダに置く。そのフォルダにカレントディレクトリを移して、リスト２のスクリプトファイルをpythonコマンドで次のように実行する。
(stan) ****/OneFacfiles$ python OneFac3PLogNorm.py
Input data
file(*.xlsx) = data.xlsx
上のPythonコマンドを実行すると、入力データファイルの設定が求められるので、設定する。上の例では、図１のファイル名が設定されているが、同じ形式のものであれば他のファイルでもよい。データ数、水準数も図１以外のものでもよい。
ファイル名を設定してEnterキーを押すと、Stanスクリプトがコンパイルされ、MCMCサンプリングが行われる。MCMCサンプリングが終了すると、トレース図が表示される（図２）。
図２
図２のWindowを閉じると、図３のグラフが表示される。
図３
各水準における平均（式（２））、標準偏差（式(3)の平方根）、モード（式(4)）、およびパラメータ
の事後分布のグラフである。
平均値の事後分布（左上のグラフ）は、V1とV3が重なり、V2とV4が重なっていて、それらは明確に分離されている。すなわち、平均値は、（V1＝V3）＜（V2＝V4）の関係にあると言える。
なお、対数変換されたデータの平均値と元のデータの平均値は、一方では同じ平均値であっても他方では異なることがある。この例は、別のウェブサイトで示した。本ウェブサイトにおける平均値meanは、変換を行わない元データの平均値である。
右上の標準偏差（分散の平方根）のグラフは、V1,V2,V4が重なっていて、それらとV3は離れている。すなわち、標準偏差は（V1＝V2＝V4）＜V3の関係にあると言える。
左下のモードのグラフでは、V1とV3のグループとV2とV4のグループに分かれていて、V3はV1より少し左に位置している。V2とV4は、ほぼ重なっている。V1とV3の平均が同じ（左上のグラフ）で、モードが違う（左下のグラフ）ことは、標準偏差の違い（右上のグラフ）に対応している。
右下のグラフは、
の事後分布である。
図３のWindowを閉じると、スクリプトの実行終了である。V1、V2、V3はほぼ重なっているが、V4は右に離れている。
図３のWindowを閉じると、スクリプトの実行終了である。
なお、現在（2026.07.15）、matplotlib 3.11で実行すると、エラーが出る。matplotlib 3.10では、エラーは出なかった。matplotlib 3.10に換えることは、次のコマンド
conda install matplotlib=3.10
で出来た。
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.
Winer, B. J., Brown, D.R., & Michels, K.M. (1991). Statistical principles in experimental design, third ed. McGraw-Hill, Inc.
リスト１ １要因データのベイズ分析Stanスクリプト（lognormal3p.stan）
data {
int N;
int J;
array[J] int nj;
array[N, J] real X;
}
transformed data {
array[J] real min_x;
for (j in 1:J) {
min_x[j] = X[1][j];
for
(i in 2:nj[j]) {
if (min_x[j] > X[i][j])
{
min_x[j] = X[i][j];
}
}
}
}
parameters {
array[J] real mu;
array[J] real<lower=0,
upper=1> vsgm;
array[J] real<lower=0,
upper=1> vtheta;
}
transformed parameters {
array[J] real sgm;
array[J] real theta;
for (j in 1:J) {
sgm[j] = 0.0001 + vsgm[j] * 1000;
theta[j] = vtheta[j] * min_x[j];
}
}
model {
for (j in 1:J) {
mu[j] ~ uniform(-1000, 1000);
vsgm[j] ~ beta(1, 1);
vtheta[j] ~ beta(1, 1);
}
real v;
for (j in 1:J) {
for
(i in 1:nj[j]) {
v = X[i][j] - theta[j];
v ~ lognormal(mu[j], sgm[j]);
}
}
}
generated quantities {
array[J] real mean;
array[J] real sd;
array[J] real mode;
for (j in 1:J) {
mean[j] = theta[j] + exp(mu[j] + 0.5*(sgm[j]^2));
sd[j] = (exp(2*mu[j] + (sgm[j]^2))
* (exp(sgm[j]^2) - 1.0))^0.5;
mode[j] = theta[j] + exp(mu[j] - (sgm[j]^2));
}
}
リスト２ 1要因デザインデータのベイズ分析Pythonスクリプト（OneFac3PLogNorm.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
df_data = pd.read_excel(input("Input
data file(*.xlsx) = "), header=0)
var_names = [v for v in df_data.keys()]
df_data_np = np.array(df_data.values)
N, J = np.shape(df_data_np)
temp_data = []
Ns = []
for v in df_data_np.T:
v = np.array(v)
v = v[np.isnan(v)
== False]
Ns.append(len(v))
if len(v)
< N:
v = np.concatenate((v, [-1]*(N - len(v))))
temp_data.append(v)
X = np.array(temp_data).T
print('Ns =', Ns)
Data = {'N':N, 'J':J, 'nj':Ns, 'X':X}
sm = CmdStanModel(stan_file='lognormal3p.stan')
fit = sm.sample(data=Data)
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())
mean_smpls = []
sd_smpls = []
mode_smpls = []
theta_smpls = []
for j in range(J):
mean_smpls.append(fit_dframe[f"mean[{j+1}]"])
sd_smpls.append(fit_dframe[f"sd[{j+1}]"])
mode_smpls.append(fit_dframe[f"mode[{j+1}]"])
theta_smpls.append(fit_dframe[f"theta[{j+1}]"])
plt.figure(figsize=(14,8))
plt.subplot(221)
plt.title('Mean', fontsize=16)
for j in range(J):
sb.kdeplot(mean_smpls[j], label=var_names[j])
plt.xlabel('mean', fontsize=14)
plt.legend()
plt.subplot(222)
for j in range(J):
sb.kdeplot(sd_smpls[j], label=var_names[j])
plt.legend()
plt.title('SD', fontsize=16)
plt.xlabel('sd', fontsize=14)
plt.subplot(223)
for j in range(J):
sb.kdeplot(mode_smpls[j], label=var_names[j])
plt.legend()
plt.title('Mode', fontsize=16)
plt.xlabel('mode', fontsize=14)
plt.subplot(224)
for j in range(J):
sb.kdeplot(theta_smpls[j], label=var_names[j])
plt.legend()
plt.title('Theta', fontsize=16)
plt.xlabel('$\\theta$', fontsize=14)
plt.tight_layout()
plt.show()