Up

正規分布と対数正規分布

ベイズ分析による比較

正規分布は統計分析において広く用いられており、分布の平均値近傍の値の観測される頻度が高いと期待される。これに対して、対数正規分布では、モード（最頻値）は平均値より小さいので、平均値より小さい値の観測される頻度が高いという印象があるかもしれない。得られたデータの分布が、正規分布からのものと考えられるのか、あるいは対数正規分布からのものと考えられるのかを、ベイズ分析におけるLOO-CVにより比較した。心理学においては、正規分布は特別な扱いを受けているので、正規分布との比較を取り上げた。また、平均値の比較において対数正規分布とその対数変換された正規分布とで平均値の差の有無に違いが出ることがあるので、対数正規分布に従うデータは、対数変換せずにそのまま対数正規分布により分析することが望ましいと考えられる。

モデルの予測精度を表すベイズ的指標としてLOO-CV(leave-one-out cross-validation)があげられる（Gelman et al., 2014; Lambert, 2018）。予測精度（Expected log pointwise predictive density for a new dataset）は、以下のように与えられる（Silva and Zanella, 2024; Vehtari et al., 2016）。

LOO-CVの場合：

である。は、データからを除いたに基づく事後分布モデルによるの予測確率である。が大きい（０に近い）ことは、が全体として大きいことを意味するので、モデルのデータへの当てはまりが良いことを示していると考えられる。

確率は、データを除いたについて求める（Leave One Out）ので、にわたって計算することになり時間が掛かる。Pareto smoothed importance sampling (PSIS) では、次式により計算時間の短縮が図られている。

ここで、はvector weight of Pareto smoothed importance samplingである。

これにより、の計算値は次式で与えられる。

LOO-CVは、arvizのloo関数により得ることができる。

まず、Stanスクリプトにおいて変数log_likを用意して対数密度確率を設定する。リスト１の場合、以下のようになっている。

generated quantities {

array[n] real log_lik; // for arviz.loo

for (i in 1:n) {

log_lik[i] = normal_lpdf(x[i]|mu, sgm);

}

}

MCMCサンプルをArvizのデータ型に変換するときにパラメータlog_likelihood=Trueを指定する。Arviz用に変換したデータをarviz.looの引数とすればLOO-CVを得ることができる。リスト３では、以下のようになっている。

fit_n = sm_n.sample(data=Data)

i_data_n = az.from_cmdstanpy(fit_n, log_likelihood=True)

loo_results_n = az.loo(i_data_n)

以下では、まずモデルの簡単な説明を行い、その後、スクリプトの実行例を、正規分布が適合する場合と対数正規分布が適合する場合について示す。ファイルは、normlognormfiles.zipにまとめた。

モデル

正規分布は次式で与えられるもので、

よく知られている。

これに対して、確率変数の対数が平均、分散の正規分布に従うとき、確率変数は対数正規分布に従うという（Johnson et al., 1994）。確率密度関数は次式（１）のように書くことができる。

このとき、平均値、分散、モードは次式で与えられる（Gelman et al., 2014）。

正規分布のときは、平均値とモードは等しいが、対数正規分布のときは

である。

対数正規分布の平均値と標準偏差は、次式

で表されるので、

である。

正規分布のStanスクリプトをリスト１に、対数正規分布のStanスクリプトをリスト２に示すように作成した。リスト１のStanスクリプトファイルとリスト２のStanスクリプトファイルを用いてベイズ分析を行うPythonスクリプトをリスト３に示すように用意した。ファイルは、normlognormfiles.zipにまとめた。

実行例

実行においては、リスト１，リスト２、リスト３のファイル、および入力データファイルは、同じフォルダに置かれているものとする。

正規分布の方が適合した場合

図１に示すデータを分析した。

・

・

・

図１　仮想データファイルdata_norm.csv

入力データファイルは、CSV形式で用意する。１行目に変数名を書き、１列目はデータの識別番号を置く。分析では用いないので、任意の数値でよい。分析のデータ値は、２列目に２行目以降に書く。

リスト３のスクリプトpython ba_anal.pyを実行すると、入力データファイル名を聞いてくる。

(stan) ****/normlognormfiles$ python ba_anal.py

Input Data File (*.csv) = data_norm.csv

上の例では、図１のファイル名data_norm.csvを設定している。

入力データファイル名を設定して、Enterキーを押すと、ファイルが読み込まれて計算が始まる。

計算が終了すると、分析結果が図２のように表示される。

図２

左側のグラフが、正規分布モデルによる分析結果であり、右側のグラフが対数正規分布による分析結果である。

モデルとデータの適合度の指標であるLOO-CVの値は正規分布の方が大きい（０に近い）ので、正規分布の方がデータに適合していると言える。

図２のWindowを右上角のX印をクリックして閉じると、次のメッセージが端末に出力されて実行終了となる。

FigNormLogNorm.png was saved.

ファイルFigNormLogNorm.pngには、図２のグラフが保存されている。

対数正規分布の方が適合した場合

図３に示すデータを分析した。

・

・

・

図３　仮想データファイルdata_lognorm.csv

入力データファイルは、CSV形式で用意する。１行目に変数名を書き、１列目はデータの識別番号を置く。分析では用いないので、任意の数値でよい。分析のデータ値は、２列目に２行目以降に書く。

リスト３のスクリプトba_anal.pyを実行すると、入力データファイル名を聞いてくる。

(stan) ****/normlognormfiles$ python ba_anal.py

Input Data File (*.csv) = data_lognorm.csv

上の例では、図３のファイル名data_lognorm.csvを設定している。

入力データファイル名を設定して、Enterキーを押すと、ファイルが読み込まれて計算が始まる。

計算が終了すると、図４のグラフが表示される。

図４

左側のグラフが、正規分布モデルによる分析結果であり、右側のグラフが対数正規分布による分析結果である。

モデルとデータの適合度の指標であるLOO-CVの値は、正規分布の方より対数正規分布の方が大きい（０に近い）ので、対数正規分布の方がデータに適合していると言える。

図４の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.

Silva, L. A., & Zanella, G. (2024). Robust Leave-One-Out Cross-Validation for High-Dimensional Bayesian Models. Journal of the American Statistical Association, 119(547), 2369–2381.

Vehtari,A., Gelman,A., & Gabry,J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27,1413-1432.

スクリプト

スクリプトは、リスト１とリスト２に分けられているが、これらはデータファイルとともに同じフォルダ内と置かれているものとする。本サイトで用いたスクリプトファイルとデータファイルは、まとめて圧縮ファイルnormlognormfiles.zipとした。

リスト１　正規分布のStanスクリプト（normalmodel.stan）

//

// 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.loo

for (i in 1:n) {

log_lik[i] = normal_lpdf(x[i]|mu, sgm);

}

}

リスト２　対数正規分布のStanスクリプト（lognormalmodel.stan）

data {

int n;

array[n] real x;

}

transformed data {

real log_max_x;

real log_min_x;

log_max_x = log(max(x));

log_min_x = log(min(x));

}

parameters {

real<lower=log_min_x, upper=log_max_x> mu;

real<lower=0.0001, upper=log_max_x - log_min_x> sgm;

}

model {

mu ~ uniform(log_min_x, log_max_x);

sgm ~ uniform(0.0001, log_max_x - log_min_x);

x ~ lognormal(mu, sgm);

}

generated quantities {

real mn;

real sd;

real mode;

array[n] real log_lik; // log likelihood for WAIC and LOO

mn = exp(mu + 0.5*(sgm^2));

mode = exp(mu - (sgm^2));

sd = (exp(2*mu + (sgm^2)) * (exp(sgm^2) - 1.0))^0.5;

for (i in 1:n) {

log_lik[i] = lognormal_lpdf(x[i] | mu, sgm);

}

}

リスト３　ベイズ分析のPythonスクリプト（ba_anal.py）

#

# Yasuharu Okamoto, 2026.06

#

from cmdstanpy import CmdStanModel

import pandas as pd

import arviz as az

import matplotlib.pyplot as plt

import numpy as np

import scipy.stats as ss

pd_df = pd.read_csv(input('Input Data File (*.csv) = '))

print(pd_df)

x = pd_df[pd_df.keys()[1]].values

print(x)

Data = {'n':len(x), 'x':x}

sm_n = CmdStanModel(stan_file='normalmodel.stan')

fit_n = sm_n.sample(data=Data)

print(fit_n.summary())

i_data_n = az.from_cmdstanpy(fit_n, log_likelihood=True)

print('Summary...\n', az.summary(i_data_n, var_names=['mu', 'sgm']))

fit_dframe_n = fit_n.draws_pd()

loo_results_n = az.loo(i_data_n)

print('\nNormal model:\n', loo_results_n)

sm_l = CmdStanModel(stan_file='lognormalmodel.stan')

fit_l = sm_l.sample(data=Data)

print(fit_l.summary())

i_data_l = az.from_cmdstanpy(posterior=fit_l, log_likelihood=True)

print('Summary...\n', az.summary(i_data_l))

fit_dframe_l = fit_l.draws_pd()

loo_results_l = az.loo(i_data_l)

print('\nLognormal model:\n', loo_results_l)

"""

Using a normal distribution model

"""

print()

print('\nmu...\n', fit_dframe_n.describe()['mu'].T)

print('\nsgm...\n', fit_dframe_n.describe()['sgm'].T)

print('\nloo\n', loo_results_n)

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_n['mu'])

v_sgm = np.median(fit_dframe_n['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_n.elpd:.1f}', fontsize=16)

plt.legend(fontsize=14)

plt.tight_layout()

"""

Using a lognormal distribution model

"""

def lognormal_pdf(x, mu, sgm):

#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

print(pd_df)

x = pd_df[pd_df.keys()[1]].values

print(x)

print()

print('\nmean...\n', fit_dframe_l.describe()['mn'].T)

print('\nmode...\n', fit_dframe_l.describe()['mode'].T)

print('\nsd...\n', fit_dframe_l.describe()['sd'].T)

print('\nloo\n', loo_results_l)

plt.subplot(122) # A graph based on a lognormal distribution model

plt.hist(x, density=True, label='Data')

v_mu = np.median(fit_dframe_l['mu'])

v_sgm = np.median(fit_dframe_l['sgm'])

#v_theta = np.median(fit_dframe_l['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_l['mn'])

v_mode = np.median(fit_dframe_l['mode'])

v_sd = np.median(fit_dframe_l['sd'])

#v_theta = np.median(fit_dframe_l['theta'])

plt.plot(xcoord, ycoord, label=fr'mean={v_mean:.1f}' + '\n' +\

fr'mode={v_mode:.1f}' + '\n' +\

fr'sd={v_sd:.1f}', lw=3)

plt.xlabel(pd_df.keys()[1], fontsize=16)

plt.title(fr'Log Normal Distribution' + '\n'

fr'LOO-CV={loo_results_l.elpd:.1f}', fontsize=16)

plt.legend(fontsize=14)

plt.tight_layout()

plt.savefig('FigNormLogNorm.png')

plt.show()

print('\nFigNormLogNorm.png was saved.\n')

Home