対数正規分布に従う乱数の作成

目次

●対数正規分布
●乱数の特徴をみてみる
●平均の等しい2標本の差の検定の比較
●中央値の等しい2標本の差の検定の比較
●結論

対数正規分布

　対数正規分布を，対数変換すると正規分布になります．

　昔は，よく「データを対数変換してt検定を適用する」などといいましたが，取得したデータが対数正規分布に従うから，対数変換をして正規分布に従うデータに変える，という意味があるのです¹．

　ノンパラメトリックな手法のシミュレーションのために，正規分布に従わない分布の代表として対数正規分布を使おうかと思いますが，対数正規分布の確率変数$X$の平均$E(X)$と分散$V(X)$は，対応する正規分布の平均$\mu$と分散$\sigma^2$と，

\begin{align}
E(X)&=e^{\mu+\frac{\sigma^2}{2}}  \\　
V(X)&= e^{2\mu+\sigma^2}(e^{\sigma^2}-1) \\ 
\end{align}

の関係があります．
　さらに中央値$E_{Med}(X)$は，

\begin{align}
E_{Med}(X)&=e^{\mu+\frac{\sigma^2}{2}}  \\　
\end{align}

の関係があります．

正規分布と対数正規分布のパラメータ

　平均$\mu$と分散$\sigma^2$を指定した対数正規分布に従う$E(X)$，$V(X)$の乱数$X$は，以上で述べた関係式を利用して発生できます．

対数正規分布の平均と標準偏差

　対数正規分布に従う平均$E(X)$，分散$V(X)$は，

\begin{align}
E(X)&=e^{\mu+\frac{\sigma^2}{2}} \tag{1} \\　
V(X)&= e^{2\mu+\sigma^2}(e^{\sigma^2}-1) \tag{2} \\ 
\end{align}

の関係があると述べました．以降は，$E(X)=E，V(X)=V^2$として表します．
　$\mu$と$\sigma^2$から$E$と$V^2$を求める方法を述べます．

平均

　$(1)$式を変形して，

\begin{align}
E&=e^{\mu+\frac{\sigma^2}{2}}  \tag{(1)式再掲} \\　
log(E)&={\mu+\frac{\sigma^2}{2}}  \\　
\therefore　\mu&={log(E)-\frac{\sigma^2}{2}}  \tag{3} \\　
\end{align}

標準偏差

　$(1)$を2乗して，

\begin{align}
E^2&= \Bigl(e^{\mu+\frac{\sigma^2}{2}} \Bigr)^2 \\
E^2&= e^{2\mu+\sigma^2} \tag{4} \\ 
\end{align}

$(4)$式を，$(2)$式に代入します．

\begin{align}
V^2&= E^2(e^{\sigma^2}-1) \\
\frac{V^2}{E^2}&= e^{\sigma^2}-1 \\ 
e^{\sigma^2}&=\frac{V^2}{E^2}+1 \\
\sigma^2&=log\biggl(\frac{V^2}{E^2}+1 \biggr) \\
\therefore　\sigma&=\sqrt{log\biggl(\frac{V^2}{E^2}+1 \biggr)} \tag{5}  \\
\end{align}

対数正規分布の中央値と標準偏差

中央値$E_{Med}(X)$，それに対する分散$V(X)$として，

\begin{align}
E_{Med}(X)&=e^{\mu} \tag{6} \\　
V(X)&= e^{2\mu+\sigma^2}(e^{\sigma^2}-1) \tag{7} \\ 
\end{align}

です．改めて以降は，$Med(X)=Med，V(X)=V^2$として表します．
　$Med$と$V^2$から，対数正規分布の特性値$\mu$と$\sigma$を求めます．

中央値

　$(6)$式を変形して，

\begin{align}
Med&=e^{\mu}  \tag{(6)式再掲} \\　
log(Med)&={\mu}  \\　
\therefore　\mu&=log(Med)  \tag{8} \\　
\end{align}

標準偏差

　公式$(a^m)^n=a^{nm}，a^n・a^m=a^{n+m}$から，$(7)$式に$(6)$式を代入，変形します．

\begin{align}
V^2&= (e^\mu)^2・e^{\sigma^2}(e^{\sigma^2}-1)  \\
V^2&= Med^2・e^{\sigma^2}(e^{\sigma^2}-1)  \\
\frac{V^2}{Med^2}&= e^{\sigma^2}(e^{\sigma^2}-1)  \\
(e^{\sigma^2})^2-e^{\sigma^2}&=\frac{V^2}{Med^2}  \\
\end{align}

ここで，全体に4をかけます．

\begin{align}
4(e^{\sigma^2})^2-4e^{\sigma^2}&= 4\frac{V^2}{Med^2}  \\
\end{align}

まわりくどいですが，左辺を変形します$(4a^2-4a+1=(2a^2-1)^2$より）．

\begin{align}
\Big(2(e^{\sigma^2})^2-1 \Bigr)^2&= 4\frac{V^2}{Med^2}+1  \\
2(e^{\sigma^2})^2-1&= \sqrt{4\frac{V^2}{Med^2}+1}  \\
e^{\sigma^2}&= \frac{\sqrt{4\frac{V^2}{Med^2}+1}+1}{2}  \\
\sigma^2&= log\Bigg( \frac{\sqrt{4\frac{V^2}{Med^2}+1}+1}{2} \Bigg)  \\
\therefore　\sigma&=\sqrt{ log\Bigg( \frac{\sqrt{4\frac{V^2}{Med^2}+1}+1}{2} \Bigg) } \tag{9}\\
\end{align}

$(3)，(5)，(8)，(9)$式を使って，乱数を作成します．

Rを用いて作成

　Rでは対数正規分布の乱数をrlnorm関数で発生できます．
　

平均と標準偏差を基に発生

$(3)，(5)$式に基づき，設定したい平均をmean=$E$，標準偏差をsd=$V$として，variableに代入します．

mean.log<- log(mean) - (sdlog^2) / 2
sd.log <- sqrt(log((sd/mean)^2 + 1))
variable<-rlnorm(n, meanlog = mean.log, sdlog = sd.log)

　

中央値と標準偏差を基に発生

$(8)，(9)$式に基づき，設定したい平均をmedian=$Med$，標準偏差をsd=$V$として，variableに代入

mean.log <- log(median)
sd.log <- sqrt(log((1 + sqrt((2*sd/median)^2 + 1))/2))
variable<-rlnorm(n, meanlog = mean.log, sdlog = sd.log)

早速，実験開始．

乱数の特徴をみてみる

平均を指定した場合

　上記の$(3),(5)$式で，$E=10$，$V＝10$の対数正規乱数をn=100,000発生させてヒストグラムを描くと，

　ほぼ対数正規分布そのまま．いずれにせよ，平均とsd，もしくは組み合わせを使って計算するだけなので，2標本が同一の分布であれば，何れの検定も理論に近似するはず．
　ということは…sdを変えると分布が変わる可能性がある（図は，それぞれn=10,000で，うえからsd=10，20，30）．

　分布が変わってます．

E＝10，V＝10，20，30のそれぞれn=10,000の乱数を発生させその平均を回収します．その結果（うえからsd=10，20，30），

そんなに悪くはない．
　それでは，n=10でやってみます．

　標本の大きさが小さくなると，標本平均の分布は対数正規分布に近くなります．平均は10と計算されますが，低い値の方に偏り，ときには大きな値の出現する確率が高まります．2標本$X_m$,$Y_n$の$n：m$または$sd_x：sd_y$の比が異なるほど期待値10から外れた値が出やすくなり，さらに$n+m$の大きさが大きいほど，$\alpha$は大きくなるでしょう．

中央値を指定した場合

上記の$(8),(9)$式で，$E_{Med}=10$，$V＝10$の対数正規乱数をn=100,000発生させてヒストグラムを描くと，

基本的に平均と同じです．$sd$を変えても同様です．
　Median＝10，V＝10，20，30のそれぞれn=10,000の乱数を発生させその平均を回収します．その結果（うえからsd=10，20，30），

　平均よりは，ばらつきが少ない．
　それでは，n=10でやってみます．

　圧倒的に平均よりも安定しています．まあ，正規分布ではないので中央値の方が安定して当然です．

平均の等しい2標本の差の検定の比較

　まずは，正規分布に従うか否か，等分散性が成り立つかどうかは無視して，対数正規分布に従う2標本のデータに対して，

• 2標本$t$検定（正規分布．等分散）
• $Welch$の検定（正規分布．等分散に寄らない）
• $Mann-Whitney$の検定（分布に寄らない．等分散）
• $Brunner-Munzel$の検定（分布に寄らない．等分散に寄らない）

を適用します．$Brunner-Munzel$の検定は，特に頑健な手といわれます．
以降のシミュレーションの有意水準は$p=0.05$とします．

シミュレーションの手順

　シミュレーションは，$1-\alpha$，検出力（$1-\beta$）の2通り行います．$1-\alpha$は，相変わらず個人の好みで$\alpha$を算出します．

①乱数の発生

　R4.2.3を使用します．
　$X_1,…,X_n～LN(\mu_x，\sigma_x^2)，Y_1,…,Y_m～LN(\mu_y，\sigma_y^2)$の対数正規乱数を発生させます．以下の各条件とします．
　$n：m=10：10，20，30，50$（4通り）
　$V_x：V_y＝10：10，20，30，40，50，60 $（6通り）
に変化させます．検出力では，$E_y=20$とします．これらを，$\mu_x，\sigma_x，\mu_y，\sigma_y$に変換して代入します．
$XとY$をvariable列，$X＝0，Y＝1$としたgroup列（as.factor）のdataを作成します．
ややこしいので，以降，$V_x$を$sd(X)$，$V_y$を$sd(Y)$と記載します．

②検定の手順

dataに対して，2標本の差の検定を行います．
（１）$X，Y$に対して，2標本$t$検定
t.test(variable~groups, alternative='two.sided', conf.level=.95, var.equal=TRUE, data=data)
を行います．→p値を回収
（２）$X，Y$に対して，$Welch$の修正による2標本t検定
t.test(variable~groups, alternative='two.sided', conf.level=.95, var.equal=FALSE, data=data)
を行います．→p値を回収
（３）$X，Y$に対して，$Mann-Whitney$の検定
wilcox.exact(variable~groups,data,paired=F)を行います．→p値を回収
（４）$X，Y$に対して，$Brunner-Munzel$の検定
brunner.munzel.test(variable,group)を行います．→p値を回収
　これらを，各10,000回繰り返し計算します．$p<0.05$の数を数え，割合を求めます．10,000回検定して500回分，5％（＝0.05）の$\alpha$エラーが理想です．検出力では，1に近づくほど良いと判断します．

注意
　Rで$Brunner-Munzel$の検定を行うためには，brunnermunzelパッケージを使用します．しかしバージョンによっては，インストールされないときもあります．代わりにlawstatパッケージでも使えます．

αエラーの結果

　　以上の結果です．0.05（5％）が理論値です．これに近いほど良い（理論通り）といえます．
　もともと正規分布に従う標本ではないので，2標本$t$検定，$Welch$の検定は，だめですね．

結果（表1）の要約

• $sd(X)：sd(Y)=10：10$のところでは，$Welch$の検定以外は，まあまあといったところ．
• $sd(X)：sd(Y)=10：20$を超えると，全て理論値を超えてしまう．差が開くに従って，どんどん0.05を上回っていく．有意差が出やすくなっていくわけです．
  • そもそも， 2標本$t$検定と$Welch$の検定は正規分布を前提にしているし，$Mann-Whitney$の検定と$Brunner-Munzel$の検定は平均の差の検定ではありませんから，当たり前といえば当たり前ですが．

検出力（1-β）の結果

　$mean(X)=10，mean(Y)=20$でシミュレーションしています．

結果（表2）の要約

• $sd(X)：sd(Y)=10：10～20$のところでは，全般に良好ではある．
• $sd(X)：sd(Y)=10：30$を超えると，検出力は下がる．しかし$sd$の差が開くに従って，検出力が上がる．
  • 今回発生させている，分布の特徴が出ています．

中央値の等しい2標本の差の検定の比較

　上記同様に，対数正規分布に従う2標本のデータに対して，

• 2標本$t$検定（正規分布．等分散）
• $Welch$の検定（正規分布．等分散に寄らない）
• $Mann-Whitney$の検定（分布に寄らない．等分散）
• $Brunner-Munzel$の検定（分布に寄らない．等分散に寄らない）

を適用します．
以降のシミュレーションの有意水準は$p=0.05$とします．

シミュレーションの手順

　シミュレーションは，$1-\alpha$，検出力（$1-\beta$）の2通り行います．$1-\alpha$は，$\alpha$を算出します．

①乱数の発生

　R4.2.3を使用します．
　中央値の等しい対数正規乱数を発生させます．以下の各条件とします．
　$n：m=10：10，20，30，50$（4通り）
　$V_x：V_y＝10：10，20，30，40，50，60 $（6通り）
に変化させます．検出力では，$E_{Med}_y=20$とします．これらを，$\mu_x，\sigma_x，\mu_y，\sigma_y$に変換して代入します．
$XとY$をvariable列，$X＝0，Y＝1$としたgroup列（as.factor）のdataを作成します．
ややこしいので，以降，$V_x$を$sd(X)$，$V_y$を$sd(Y)$と記載します．

②検定の手順

dataに対して，2標本の差の検定を行います．
（１）$X，Y$に対して，2標本$t$検定
t.test(variable~groups, alternative='two.sided', conf.level=.95, var.equal=TRUE, data=data)
を行います．→p値を回収
（２）$X，Y$に対して，$Welch$の修正による2標本t検定
t.test(variable~groups, alternative='two.sided', conf.level=.95, var.equal=FALSE, data=data)
を行います．→p値を回収
（３）$X，Y$に対して，$Mann-Whitney$の検定
wilcox.exact(variable~groups,data,paired=F)を行います．→p値を回収
（４）$X，Y$に対して，$Brunner-Munzel$の検定
brunner.munzel.test(variable,group)を行います．→p値を回収
　これらを，各10,000回繰り返し計算します．$p<0.05$の数を数え，割合を求めます．10,000回検定して500回分，5％（＝0.05）の$\alpha$エラーが理想です．検出力では，1に近づくほど良いと判断します．

αエラーの結果

　　以上の結果です．0.05（5％）が理論値です．これに近いほど良い（理論通り）といえます．
　もともと正規分布に従う標本ではないので，2標本$t$検定，$Welch$の検定は，だめですね．

結果（表1）の要約

• $n=m$のときは，$Mann-Whitney$の検定はほぼ完ぺきです．
• $Brunner-Munzel$の検定は，さすがに優秀です．

検出力（1-β）の結果

　$median(X)=10，median(Y)=20$でシミュレーションしています．

結果（表2）の要約

• どの手法も，似たり寄ったりです．

結論

　対数正規分布に従う標本に対して，シミュレーションを行いました．

• 対数正規分布の標本では，平均を基にした考え方が間違っている．
  • 今更ですが，正規分布に従うか従わないかの確認は，重要と思います．
• 対数正規分布の標本では，中央値を基にする．
  • 対応する検定も，$Mann-Whitney$の検定，$Brunner-Munzel$の検定が望ましい．
  • $Brunner-Munzel$の検定は，より頑健である．積極的に使用する価値はある．
  • ただし，上記両者とも”中央値の差の検定”ではないことも留意しておく．

　やはり事前に正規分布に従うか否かを確認した方が，良いと思います．どうやって確認するかですが，検定に頼らざるを得ないのではないかと思います．

1. 取得したデータが対数正規分布に従うならば，なので，何でもかんでも対数変換というわけにはいきません． ↩