対数正規分布と正規分布での平均値の比較

対数正規分布に従う確率変数は、対数変換すると正規分布に従う確率変数になる。一方での平均値に統計学的差が認められても、他方では差が認められないことがある。対数正規分布において平均値に差があるときに、その対数変換したデータでは差が認められない例と、対数正規分布では平均値に差がないがその対数変換した正規分布では差が認められる例を用意した。ファイルは、normallognormalfiles.zipにまとめた。

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

よく知られている。

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

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

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

である。

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

で表されるので、

である。

２つの独立なサンプルの対数正規分布による分析は、正規分布による分析と同じように、Stanを用いるとベイズ分析を簡単に行うことができる。正規分布による分析と対数正規分布による分析のStanスクリプト例をリスト１およびリスト２のように用意した。

対数正規分布において平均値に差が認められるが、対数変換した正規分布では差が認められない例

図１のファイルの変数x1とx2は、平均値に統計学的差が認められるものである。

図１　サンプルデータdata_diffMean.xlsx

図１のデータの変数x1とx2のヒストグラム、およびlog(x1)とlog(x2)のヒストグラムを描くと、図２のようになる。

図２　データとその対数変換値のヒストグラム

左側のグラフがx1とx2のグラフであり、それぞれ平均値が168.16と242.08であることがグラフの右上に表示されている。

対数変換された変数log(x1)とlog(x2)のヒストグラムが右側に示されている。log(x1)の平均値が5.00、log(x2)の平均値が5.00と表示されている。

左のグラフにおける平均値の差と、右側のグラフにおける平均値の差をベイズ分析するためのスクリプトをリスト３のように用意した。

変数x1とx2に対数正規分布を当てはめてベイズ分析を行い、それぞれの平均値を表すパラメータmn1とmn2の事後分布を図３の左側のグラフに示した。

図３　データとその対数変換値のベイズ分析

mn1とmn2の事後分布は明瞭に分離されている。パラメータmn1の95％CIは[166.46, 169.98]であり、パラメータmn2の95%CIは[236.17, 248.19]であり、２つの区間は重なりがない。事後分布においてmn1がmn2より小さい比率はであり、事後分布においてmn1の方がmn2より小さいと言える。

対数変換されたデータlog(x1)とlog(x2)に適用された正規分布モデルの平均値パラメータmu1とmu2の事後分布が図３の右側に示されている。事後分布の重なりは明瞭である。パラメータmu1の95%CIは[4.99, 5.01]であり、mu2の95%CIは[4.98, 5.02]であり、mu1の95%CIを含んでいる。事後分布において、mu1がmu2より小さい比率はであるので、mu1とmu2は同じぐらいの値であると考えられる。

対数正規分布において平均値に差が認められないが、対数変換した正規分布では差が認められる例

図４のデータdata_sameMean.xlsxについて考える。

図４　サンプルデータdata_sameMean.xlsx

ヒストグラムを描くと図５のようになる。

図５　データとその対数変換後のヒストグラム

変量x1とx2の平均はほぼ同じであるが、対数変換した変量nx1=log(x1)とnx2=log(x2)の平均値は、3.56と2.52であり異なる。この平均値の比較を統計学的に検討するために、リスト４のスクリプトによりベイズ分析を行った。

対数正規分布の変量x1とx2の平均値パラメータmn1とmn2、および対数変換した変量log(x1)とlog(x2)の正規分布の平均値パラメータmu1とmu2の事後分布は、図６に示すものが得られた。

図６　データとその対数変換値のベイズ分析

右側のグラフに示されている対数変換された変量log(x1)とlog(x2)の正規分布の平均値パラメータmu1とmu2の事後分布は明瞭に分離されている。しかし、左側のグラフに示されている変量x1とx2の対数正規分布の平均値パラメータmn1とmn2の事後分布には、かなりの重なりが認められる。パラメータmn1の95%CIは[49.20, 51.12]であり、mn2の95%CIは[48.20, 53.47]であり、mn1の95%CIはMN2の95%CIに含まれている。事後分布において、パラメータmn1がmn2より小さい比率はであり、mn1とmn2に差はほぼないと言える。

以上より、変量x1とx2の平均値に差はないと考えられる。

なお、２つの独立なサンプルの、正規分布、対数正規分布、３パラメータ対数正規分布による分析のウェブサイトも用意している。

参考文献

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian Data Analysis, third edition. 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.

岡本安晴（2019）　いまさら聞けないPythonでデータ分析．丸善出版

付録

リスト１　２つの独立なサンプルの正規分布による分析スクリプト（ba_normal.stan）

data {

int n1;

array[n1] real x;

int n2;

array[n2] real y;

}

transformed data {

real min1;

real min2;

real max1;

real max2;

min1 = min(x);

min2 = min(y);

max1 = max(x);

max2 = max(y)

;}

parameters {

real<lower=min1, upper=max1> mu1;

real<lower=0.0001, upper=max1-min1> sgm1;

real<lower=min2, upper=max2> mu2;

real<lower=0.0001, upper=max2-min2> sgm2;

}

transformed parameters {

}

model {

mu1 ~ uniform(min1, max2);

sgm1 ~ uniform(0.0001, max1-min1);

mu2 ~ uniform(min2, max2);

sgm2 ~ uniform(0.0001, max2-min2);

x ~ normal(mu1, sgm1);

y ~ normal(mu2, sgm2);

}

リスト２　２つの独立なサンプルの対数正規分布による分析スクリプト（ba_lognormal.stan）

data {

int n1;

array[n1] real x1;

int n2;

array[n2] real x2;

}

transformed data {

real log_max_x1;

real log_min_x1;

real log_max_x2;

real log_min_x2;

log_max_x1 = log(max(x1));

log_min_x1 = log(min(x1));

log_max_x2 = log(max(x2));

log_min_x2 = log(min(x2));

}

parameters {

real<lower=log_min_x1, upper=log_max_x2> mu1;

real<lower=log_min_x2, upper=log_max_x2> mu2;

real<lower=0.0001, upper=log_max_x1 - log_min_x1> sgm1;

real<lower=0.0001, upper=log_max_x2 - log_min_x2> sgm2;

}

model {

mu1 ~ uniform(log_min_x1, log_max_x1);

sgm1 ~ uniform(0.0001, log_max_x1 - log_min_x1);

mu2 ~ uniform(log_min_x2, log_max_x2);

sgm2 ~ uniform(0.0001, log_max_x2 - log_min_x2);

x1 ~ lognormal(mu1, sgm1);

x2 ~ lognormal(mu2, sgm2);

}

generated quantities {

real mn1;

real mn2;

real sd1;

real sd2;

real diff_mns;

real diff_sds;

real diff_mus;

real diff_sgms;

real mode1;

real mode2;

real diff_modes;

mn1 = exp(mu1 + 0.5*(sgm1^2));

mode1 = exp(mu1 - (sgm1^2));

sd1 = (exp(2*mu1 + (sgm1^2)) * (exp(sgm1^2) - 1.0))^0.5;

mn2 = exp(mu2 + 0.5*(sgm2^2));

mode2 = exp(mu2 - (sgm2^2));

sd2 = (exp(2*mu2 + (sgm2^2)) * (exp(sgm2^2) - 1.0))^0.5;

diff_mns = mn1 - mn2;

diff_sds = sd1 - sd2;

diff_modes = mode1 - mode2;

diff_sds = sd1 - sd2;

diff_mus = mu1 - mu2;

diff_sgms = sgm1 - sgm2;

diff_modes = mode1 - mode2;

}

リスト３　データ図１（data_diffMean.xlsx）の分析スクリプト

import numpy as np

import matplotlib.pyplot as plt

import pandas as pd

from cmdstanpy import CmdStanModel

import seaborn as sb

pdframe = pd.read_excel('data_diffMean.xlsx')

print(pdframe)

print(pdframe.keys())

print(pdframe['x1'])

print(pdframe['x1'].values)

x1 = pdframe['x1'].values

x2 = pdframe['x2'].values

model = CmdStanModel(stan_file="ba_lognormal.stan") # コンパイル

fit = model.sample(data={'n1':len(x1), 'x1':x1,'n2':len(x2), 'x2':x2})

print(fit.summary()) # MCMCサンプリングの統計

fit_dframe = fit.draws_pd() # PandasのDataFrame型に変換

# 変数の対数変換

x = np.log(x1)

y = np.log(x2)

model_n = CmdStanModel(stan_file="ba_normal.stan") # コンパイル

fit_n = model_n.sample(data={'n1':len(x), 'x':x, 'n2':len(y), 'y':y})

print(fit_n.summary()) # MCMCサンプリングの統計

fit_dframe_n = fit_n.draws_pd() # PandasのDataFrame型に変換

plt.figure(figsize=(12,5))

plt.subplot(121)

p025, p975 = np.percentile(fit_dframe['mn1'], [2.5, 97.5])

sb.kdeplot(fit_dframe['mn1'], label=f'mn1: \n95$CI=[{p025:.2f}, {p975:.2f}]')

p025, p975 = np.percentile(fit_dframe['mn2'], [2.5, 97.5])

sb.kdeplot(fit_dframe['mn2'], label=f'mn2: \n95%CI=[{p025:.2f}, {p975:.2f}]')

pmn1Lmn2 = np.mean(fit_dframe['mn1'] - fit_dframe['mn2'] < 0)

plt.title(f'P(mn1 < mn2)={pmn1Lmn2:.3f}', fontsize= 16)

plt.xlabel('mn1, mn2', fontsize= 14)

plt.legend(fontsize=14)

plt.subplot(122)

p025, p975 = np.percentile(fit_dframe['mu1'], [2.5, 97.5])

sb.kdeplot(fit_dframe_n['mu1'], label=f'mu1: \n95%CI=[{p025:.2f}, {p975:.2f}]')

p025, p975 = np.percentile(fit_dframe['mu2'], [2.5, 97.5])

sb.kdeplot(fit_dframe_n['mu2'], label=f'mu2\n95%CI=[{p025:.2f}, {p975:.2f}]')

pmu1Lmu2 = np.mean(fit_dframe['mu1'] - fit_dframe['mu2'] < 0)

plt.title(f'P(mu1 < mu2)={pmu1Lmu2:.3f}', fontsize= 16)

plt.xlabel('mu1, mu2', fontsize= 12)

plt.legend(fontsize=12, loc='lower center')

plt.show()

リスト４　データ図４（data_sameMean.xlsx）の分析スクリプト