正規分布と対数正規分布
LOO-CVとWAICによる比較
正規分布は統計分析において広く用いられており、分布の平均値近傍の値の観測される頻度が高いと期待される。これに対して、対数正規分布では、モード（最頻値）は平均値より小さいので、平均値より小さい値の観測される頻度が高いという印象があるかもしれない。得られたデータの分布が、正規分布からのものと考えられるのか、あるいは対数正規分布からのものと考えられるのかを、ベイズ分析におけるLOO-CVとWAICにより比較した。心理学においては、正規分布は特別な扱いを受けているので、正規分布との比較を取り上げた。
まず、モデルの簡単な説明を行い、その後、スクリプトの実行例を、正規分布が適合する場合と対数正規分布が適合する場合について示す。スクリプトは、リスト１とリスト２として示す。
モデル
正規分布は次式で与えられるもので、
よく知られている。
これに対して、確率変数
の対数
が平均
、分散
の正規分布に従うとき、確率変数は３パラメータの対数正規分布に従うという（Johnson et al., 1994）。
のときは、単に対数正規分布という。確率密度関数は次式（１）のように書くことができる。
このとき、平均値、分散、モードは次式で与えられる。
正規分布のときは、平均値とモードは等しいが、対数正規分布のときは
である。
対数正規分布のパラメータをと
ではなく、対数正規分布の平均値
と標準偏差
をパラメータとすると、パラメータと分布との関係が分かり易くなる。
式（２）、式（３）より
とおけば、
であり、
である。
上のモデルに従ってベイズ分析を行うスクリプトを、リスト１およびリスト２として作成した。リスト１とリスト２、およびデータファイルは同じフォルダ内に置かれているものとする。リスト１のスクリプトをPython3.12/MiniConda/Ubuntu22/WSL/Windowsで実行した後、リスト２のスクリプトをPython3.13/Windowsで実行した。
図１に示すデータを分析した。
・
・
・
図１ 仮想データファイルdata_norm.csv
入力データファイルは、CSV形式で用意する。１行目に変数名を書き、１列目はデータの識別番号を置く。分析では用いないので、任意の数値でよい。分析のデータ値は、２列目に２行目以降に書く。
リスト１のスクリプトmcmcsampling.pyを実行すると、入力データファイル名を聞いてくる。
(py312) *****/normlognormfiles$ python mcmcsampling.py
Input Data File (*.csv) = data_norm.csv
上の例では、図１のファイル名data_norm.csvを設定している。
入力データファイル名を設定して、Enterキーを押すと、ファイルが読み込まれて計算が始まる。
計算が終了すると、以下のメッセージが端末に出力されて、スクリプトの実行終了となる。
fit_dframe_N.pkl
was saved.
fit_dframe_LN.pkl
was saved.
(py312)
*****/normlognormfiles$
MCMCサンプリングのサンプルは、fit_dframe_N.pklとfit_dframe_LN.pklに保存されている。これらは、次のリスト２の実行において読み込まれ、MCMCサンプルの分析が行われる。
リスト１の実行後、リスト２をPython3.13/Windowsで実行した。リスト１の実行で保存されたMCMCサンプルが読み込まれ、分析結果が図２のように表示される。
図２
左側のグラフが、正規分布モデルによる分析結果であり、右側のグラフが３パラメータ対数正規分布による分析結果である。
モデルとデータの適合度の指標であるLOO-CVとWAICの値は、いずれも正規分布の方が小さいので、正規分布の方がデータに適合していると言える。
図２のWindowを右上角のX印をクリックして閉じると、次のメッセージが端末に出力されて実行終了となる。
FigNormLogNorm.png was saved.
ファイルFigNormLogNorm.pngには、図２のグラフが保存されている。
図３に示すデータを分析した。
・
・
・
図３ 仮想データファイルdata_lognorm.csv
入力データファイルは、CSV形式で用意する。１行目に変数名を書き、１列目はデータの識別番号を置く。分析では用いないので、任意の数値でよい。分析のデータ値は、２列目に２行目以降に書く。
リスト１のスクリプトmcmcsampling.pyを実行すると、入力データファイル名を聞いてくる。
(py312) *****/normlognormfiles$
python mcmcsampling.py
Input Data
File (*.csv) = data_lognorm.csv
上の例では、図３のファイル名data_lognorm.csvを設定している。
入力データファイル名を設定して、Enterキーを押すと、ファイルが読み込まれて計算が始まる。
計算が終了すると、以下のメッセージが端末に出力されて、スクリプトの実行終了となる。
fit_dframe_N.pkl was saved.
fit_dframe_LN.pkl was saved.
(py312) *****/normlognormfiles$
MCMCサンプリングのサンプルは、fit_dframe_N.pklとfit_dframe_LN.pklに保存されている。これらは、次のリスト２の実行において読み込まれ、MCMCサンプルの分析が行われる。
リスト１の実行後、リスト２をPython3.13/Windowsで実行した。リスト１の実行で保存されたMCMCサンプルが読み込まれ、分析結果が図４のように表示される。
図４
左側のグラフが、正規分布モデルによる分析結果であり、右側のグラフが３パラメータ対数正規分布による分析結果である。
モデルとデータの適合度の指標であるLOO-CVとWAICの値は、いずれも正規分布の方より対数正規分布の方が小さいので、対数正規分布の方がデータに適合していると言える。
図４のWindowを右上角のX印をクリックして閉じると、次のメッセージが端末に出力されて実行終了となる。
FigNormLogNorm.png
was saved.
ファイルFigNormLogNorm.pngには、図４のグラフが保存されている。
参考文献
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.
スクリプトは、リスト１とリスト２に分けられているが、これらはデータファイルとともに同じフォルダ内と置かれているものとする。本サイトで用いたスクリプトファイルとデータファイルは、まとめて圧縮ファイルnormlognormfiles.zipとした。
リスト１のスクリプトを実行したのち、リスト２のスクリプトを実行する。リスト１のスクリプトは、Python3.12/MiniConda/Ubuntu22/WSL/Windowsで実行した。リスト１の実行後、リスト２のスクリプトを、Python3.13/Windowsで実行した。
#
# Yasuharu
Okamoto, 2025.08
#
import stan
import pandas as pd
import arviz as az
import pickle
stan_model = """
//
// Normal
Distribution
//
data {
int
n;
array[n] real x;
}
transformed data {
real
xmin;
real
xmax;
xmin
= min(x);
xmax
= max(x);
}
parameters {
real<lower=0, upper=1> vmu;
real<lower=0, upper=1> vsgm;
}
transformed parameters {
real
mu;
real
sgm;
mu =
xmin + (xmax - xmin) * vmu;
sgm
= 0.0001 + (xmax - xmin) * vsgm;
}
model {
vmu
~ beta(1, 1);
vsgm
~ beta(1, 1);
x ~
normal(mu, sgm);
}
generated quantities {
array[n] real log_lik;
// for arviz.waic and
loo
for (i
in 1:n) {
log_lik[i] = normal_lpdf(x[i]|mu, sgm);
}
}
"""
pd_df = pd.read_csv(input('Input Data
File (*.csv) = ')) # pd.read_csv('data_norm.csv')
print(pd_df)
x = pd_df[pd_df.keys()[1]].values
print(x)
Data = {'n':len(x), 'x':x}
sm = stan.build(stan_model, data=Data)
fit = sm.sample(num_chains=4,
init=[{'vmu':0.5, 'vsgm':0.5}]*4)
print(fit)
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print('Summary...\n',
az.summary(i_data, var_names=['mu', 'sgm']))
fit_dframe = fit.to_frame()
loo_results = az.loo(i_data,
scale='deviance')
print('\n', loo_results)
waic_results = az.waic(i_data,
scale='deviance')
print('\n', waic_results)
with open('fit_dframe_N.pkl', 'wb') as
f:
pickle.dump((pd_df, fit_dframe,
loo_results, waic_results), f) # MCMCサンプルの保存
stan_lognormal_model =
"""
//
// LogNormal Distribution
//
data {
int
n;
array[n] real x;
}
transformed data {
real
xmin;
real
xmax;
xmin
= min(x);
xmax
= max(x);
}
parameters {
real<lower=0.0, upper=1.0> vmn;
real<lower=0.0, upper=1.0> vsd;
real<lower=0.0, upper=1.0> vtheta;
}
transformed parameters {
real
mn;
real
sd;
real
mu;
real
sgm;
real
theta;
/*
パラメータのスケーリング
model宣言部で変域の広いパラメータに
出される警告を、スケーリングで避ける
*/
mn =
0.0001 + vmn * xmax;
sd =
0.0001 + vsd * xmax;
sgm
= (log((sd/mn)^2 + 1))^0.5;
mu =
log(mn) - 0.5 * (sgm^2);
theta = vtheta * xmin;
}
model {
vmn
~ beta(1, 1);
vsd
~ beta(1, 1);
vtheta ~ beta(1, 1);
for
(i in 1:n) {
x[i] - theta ~ lognormal(mu, sgm);
}
}
generated quantities {
real
mean;
real
mode;
array[n] real log_lik;
// for arviz.waic and
loo
mean
= theta + mn;
mode
= theta + exp(mu - sgm^2);
for
(i in 1:n) {
log_lik[i] = lognormal_lpdf(x[i]-theta|mu, sgm);
}
}
"""
Data = {'n':len(x), 'x':x}
sm = stan.build(stan_lognormal_model,
data=Data)
fit = sm.sample(num_chains=4, init=[{'vmu':0.5,
'vsgm':0.5, 'vtheta':0.5}]*4)
print(fit)
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print('Summary...\n',
az.summary(i_data))
fit_dframe = fit.to_frame()
loo_results = az.loo(i_data,
scale='deviance')
print('\n', loo_results)
waic_results = az.waic(i_data,
scale='deviance')
print('\n', waic_results)
with open('fit_dframe_LN.pkl', 'wb')
as f:
pickle.dump((pd_df,
fit_dframe, loo_results, waic_results), f)
# MCMCサンプルの保存
print('\nfit_dframe_N.pkl was saved.')
print('\nfit_dframe_LN.pkl was
saved.\n')
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import scipy.stats as ss
import pickle
"""
Using a normal distribution model
"""
with open('fit_dframe_N.pkl', 'rb') as
f:
(pd_df, fit_dframe,
loo_results, waic_results) = pickle.load(f)
# MCMCサンプルの保存
print(pd_df)
x = pd_df[pd_df.keys()[1]].values
print(x)
print()
print('\nmu...\n',
fit_dframe.describe()['mu'].T)
print('\nsgm...\n',
fit_dframe.describe()['sgm'].T)
print('\nloo\n', loo_results)
print('\nwaic\n', waic_results)
x = pd_df[pd_df.keys()[1]].values
print(x)
plt.figure(figsize=(10,4))
plt.subplot(121) # A graph based on a normal distribution
model
plt.hist(x, label='Data',
density=True)
v_mu = np.median(fit_dframe['mu'])
v_sgm = np.median(fit_dframe['sgm'])
xcoord = np.linspace(np.min(x),
np.max(x), 1000)
ycoord = ss.norm.pdf(xcoord, v_mu,
v_sgm)
plt.plot(xcoord, ycoord,
label=fr'$\mu$={v_mu:.1f}' + '\n' +
fr'$\sigma$={v_sgm:.1f}', lw=3)
plt.xlabel(pd_df.keys()[1],
fontsize=14)
plt.title('Normal Distribution\n' + \
fr'LOO-CV={loo_results.elpd_loo:.1f},' + \
fr'
WAIC={waic_results.elpd_waic:.1f}',
fontsize=16)
plt.legend(fontsize=14)
plt.tight_layout()
"""
Using a lognormal distribution model
"""
def lognormal_pdf(x, mu, sgm, theta):
x -= theta
if x <= 0.0:
return 0.0
else:
v1 =
sgm * x * ((2*np.pi)**0.5)
v2 =
np.exp(-((np.log(x) - mu)**2) / (2 * (sgm**2)))
return v2 / v1
with open('fit_dframe_LN.pkl', 'rb')
as f:
(pd_df, fit_dframe, loo_results,
waic_results) = pickle.load(f)
# MCMCサンプルの保存
print(pd_df)
x = pd_df[pd_df.keys()[1]].values
print(x)
print()
print('\nmean...\n',
fit_dframe.describe()['mean'].T)
print('\nmode...\n',
fit_dframe.describe()['mode'].T)
print('\nsd...\n',
fit_dframe.describe()['sd'].T)
print('\nloo\n', loo_results)
print('\nwaic\n', waic_results)
x = pd_df[pd_df.keys()[1]].values
print(x)
plt.subplot(122) # A graph based on a lognormal
distribution model
plt.hist(x, density=True,
label='Data')
v_mu = np.median(fit_dframe['mu'])
v_sgm = np.median(fit_dframe['sgm'])
v_theta =
np.median(fit_dframe['theta'])
xcoord = np.linspace(np.min(x),
np.max(x), 1000)
ycoord = []
for v in xcoord:
ycoord.append(lognormal_pdf(v,
v_mu, v_sgm, v_theta))
v_mean = np.median(fit_dframe['mean'])
v_mode = np.median(fit_dframe['mode'])
v_sd = np.median(fit_dframe['sd'])
v_theta =
np.median(fit_dframe['theta'])
plt.plot(xcoord, ycoord,
label=fr'mean={v_mean:.1f}' + '\n' +\
fr'mode={v_mode:.1f}' + '\n' +\
fr'sd={v_sd:.1f}' + '\n' +\
fr'$\theta$={v_theta:.1f}', lw=3)
plt.xlabel(pd_df.keys()[1],
fontsize=16)
plt.title(fr'Log Normal Distribution'
+ '\n'
fr'LOO-CV={loo_results.elpd_loo:.1f},' + \
fr'
WAIC={waic_results.elpd_waic:.1f}',
fontsize=16)
plt.legend(fontsize=14)
plt.tight_layout()
plt.savefig('FigNormLogNorm.png')
plt.show()
print('\nFigNormLogNorm.png was
saved.\n')