「統計学関連なんでもあり」の過去ログ--- 048

No.22635　正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法　　【Saito】　2018/11/29(Thu) 17:11

お世話になっております。
標記の通りなのですが，生データがなく正規分布を仮定した際の平均値と標準偏差しかない状況で，対数正規分布を仮定していた場合の信頼区間を算出する方法について，私のコーディングが間違っているだけかもしれませんが，うまく推定できないので，どなたかご教授頂けると幸いです。

下記に例を示します。

> set.seed(1)
> n <- 10000
> u <- 0.5
> sig <- 0.8
> x <- exp(rnorm(n, u, sig))
> 
> # 1. 平均値
> mean(x)
[1] 2.269916
> 
> # 2. 標準偏差
> sd(x)
[1] 2.123908
> 
> # 3. 信頼区間下限値
> mean(x)-1.96*sd(x)
[1] -1.892943
> 
> # 4. 信頼区間上限値
> mean(x)+1.96*sd(x)
[1] 6.432775

このとき，手元にあるのは，1-2の平均値，標準偏差，そして95%信頼区間（の近似値）である，3-4の値のみで，生データであるxは利用不可能な状況であるとします。

上記の信頼区間の算出方法ですと，xが負の値を取りえないのに，マイナスの値をとってしまっていることが確認できると思います。この部分を正の値の下限値を出すようにしたいと考えています。

一応，対数正規分布の信頼区間から正規分布の信頼区間を得る方法が使えるのかと思い，下記のURLにある計算式を打ち込んでみますと，実際の95%点とずれが出るように思われます。
http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf

# 推定値
 > mean(x)/exp(1.96*sqrt(log(1+(sd(x)/mean(x))^2)))
[1] 0.4797148
> mean(x)*exp(1.96*sqrt(log(1+(sd(x)/mean(x))^2)))
[1] 10.7408

# 実際の95%点
> quantile(x, c(0.025, 0.975))
     2.5%     97.5% 
0.3267348 8.0988751

たぶんこういう風に逆算でやるのはダメ，ということなのだと思うのですが・・・。
試しにlogの世界でやってみるとおおよそ合います。

> quantile(log(x), c(0.025, 0.975))
     2.5%     97.5% 
-1.523258  2.489656 
> mean(log(x)-1.96*sd(log(x)))
[1] -1.490756
> mean(log(x)+1.96*sd(log(x)))
[1] 2.477682

どこかで間違っているのだと思いますが，どなたか mean(x) と sd(x) だけを使って，quantile(x, c(0.025, 0.975)) で出力される95%点（負の値を取らない）を出すやり方をご存知の方がおられましたら，ご教授ください。

どうぞよろしくお願いいたします。

No.22636　Re: 正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法　　【青木繁伸】　2018/11/29(Thu) 19:46

いろいろあるので，順を追って
set.seed(1)
n <- 10000
u <- 0.5
sig <- 0.8
としたとき，
x0 <- rnorm(n, u, sig)
の x0 の平均値と標準偏差は当然ながら 0.5 と 0.8 ではないということ
> x0 <- rnorm(n, u, sig)
> mean(x0)
[1] 0.4947704
> sd(x0)
[1] 0.8098852
正確に指定された平均値と標準差を持つ正規乱数を生成するのは，一度正規化した後に，標準偏差を掛けて平均値を加える（標準化の逆をやる）。
> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8
この指数を取って対数正規分布にする。
> x2 <- exp(x1)
> mean(x2)
[1] 2.268486
> sd(x2)
[1] 2.130959
しかし，平均値と標準偏差を計算するのはそもそもおかしい。対数正規分布において計算される平均値と標準偏差のもつ意味は，正規分布のそれらが持つ意味とは全く異なる。たとえば，対数正規分布の平均値±1.96*標準偏差の範囲内にあるデータの割合は95%などではない。

また，mean(x)-1.96*sd(x) を「信頼区間下限値」としているが，それは上で書いたように，あるいは，あなたが最後の方で quantile といっていることからも，データの下から 2.5% 目の推定値でしょう。あなたは「母平均の信頼区間」を得たいのかな？ちがいますよね。

>> 試しにlogの世界でやってみるとおおよそ合います。
>> > quantile(log(x), c(0.025, 0.975))
とありますが，

あなたの対数正規乱数 x は，正規乱数の指数をとったもので，私が上で作った正確な例では x2 です。
> quantile(log(x2), c(0.025, 0.975))
     2.5%     97.5% 
-1.083910  2.063209 
はあたりまえですが，以下と同じですね。 log(x2) == x1, x2 ==exp(x1) ですから。
> quantile(x1, c(0.025, 0.975))
     2.5%     97.5% 
-1.083910  2.063209 
解はあるか？

その前に
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/wrong-sd.html
を読んでみてください。

そこにある図で，正規分布での平均値±標準偏差が対数正規分布では exp(平均値)*/exp(標準偏差) に対応すると書いてあります。

しかし，逆はできないのです。

あなたは図の上側の対数正規分布を対象にしているのですが，対数正規分布で計算した標準偏差は，図の下側の正規分布には変換するすべがないのが分かるかと思います。

解はあると思いますが，私は知らない。

http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf
は，正規分布のパラメータμとσと対数正規分布の平均値と分散の関係式であるが，解析的には成り立つが，実際のシミュレーションデータでは正確には成り立たない。
set.seed を外して何回かやってみると分かるが，x1 の平均値，標準偏差はいつも指定された値(0.5, 0.8)になるが，それを exp したとたん，x2 の平均値，標準偏差は毎回変わる（つまり，解析的な値とは異なる.）。そのため，解析的な関係式にしたがわない。
> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8
> 
> x2 <- exp(x1)
> mean(x2)
[1] 2.274771 ★これと
> sd(x2)
[1] 2.183789　●これと

> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8
> 
> x2 <- exp(x1)
> mean(x2)
[1] 2.266348　★これ
> sd(x2)
[1] 2.08723　●これ

No.22637　Re: 正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法　　【Saito】　2018/11/30(Fri) 11:32

コメントありがとうございます。
色々と間違っていたようで申し訳ありません。
質問のコメントについて考えているうちに，私の質問自体がよくなかったことに気づきましたので，最後のほうで新しく例を作りました。二度手間にしてしまって申し訳ありません。私の質問としては，データではなく信頼区間について聞きたかった（ことに気づきました）ということです。申し訳ありません。大変恐れ入りますが，ご教授くだされば幸いです。

先に頂いたコメントについてご返信させていただきます。

> x0 <- rnorm(n, u, sig)
> の x0 の平均値と標準偏差は当然ながら 0.5 と 0.8 ではないということ

についてですが，おっしゃる通りです。ただ，中心極限定理で，0.5と0.8の周りに正規分布してくれるような方法でもOKです。

> しかし，平均値と標準偏差を計算するのはそもそもおかしい。対数正規分布において計算される平均値と標準偏差のもつ意味は，正規分布のそれらが持つ意味とは全く異なる。

仰る通りです。ただ，そういう出力（正規分布を仮定した平均値と標準偏差）がすでにされてしまって（紙になってしまって）いるので，生データを抜きに，それをなんとかしようということでした。

> また，mean(x)-1.96*sd(x) を「信頼区間下限値」としているが，それは上で書いたように，あるいは，あなたが最後の方で quantile といっていることからも，データの下から 2.5% 目の推定値でしょう。あなたは「母平均の信頼区間」を得たいのかな？ちがいますよね。

すみません，mean(x)-1.96*sd(x) はご指摘の通り信頼区間ではありませんでした。また，この指摘を受けて，私も聞きたいことが聞けていないことに気づきましたので，下で改めてシミュレーションデータを作らせていただきました。二度手間になってしまい，申し訳ありません。

>その前に
>http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/wrong-sd.html
>を読んでみてください。

上に同じです。読んでいくうちに，私の聞きたいことが聞けない質問内容になっていることに気づきました。本当に申し訳ありません。

>http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf
>は，正規分布のパラメータμとσと対数正規分布の平均値と分散の関係式であるが，解析的には成り立つが，実際のシミュレーションデータでは正確には成り立たない。

すみません，新しい質問を下記でさせていただく前に，ここについてお聞かせください。解析的に成り立たなくとも，中心極限定理で近い値にはなるのかと思っていたのですが，うまく推定できません。
set.seed(1)
n <- 10000
u <- 0.5
sig <- 0.8

B <- 1000
x <- matrix(exp(rnorm(n*B, u, sig)), n, B)

xlo1 <- apply(x, 2, mean)/exp(1.96*sqrt(log(1+(apply(x, 2, sd)/apply(x, 2, mean))^2)))
xup1 <- apply(x, 2, mean)*exp(1.96*sqrt(log(1+(apply(x, 2, sd)/apply(x, 2, mean))^2)))

xlo2 <- apply(x, 2, quantile, 0.025)
xup2 <- apply(x, 2, quantile, 0.975)

hist(xlo1, xlim=c(0.3, 0.55))
hist(xlo2, add=T, col=2)

hist(xup1, xlim=c(7, 13.5))
hist(xup2, add=T, col=2)
xlo1とxlo2，xup1とxup2はほぼ一致するのかと思ったのですが，どこか誤解しているようです。

以下似たような話ですが，データとパラメータで知りたいことが違っていましたので，質問を修正致します。大変申し訳ありません。
set.seed(1)
n <- 10
u <- 1
sig <- 2

# 生物の密度データがあるとする（非負）。平均密度とその信頼区間を知りたい。
y <- exp(rnorm(n, u, sig))
hist(y)

# yというデータがあるもとで，正規分布を誤差分布に仮定したときの平均密度の95%信頼区間
res <- summary(lm(y~1))
lo1 <- res$coef[1]-1.96*res$coef[2]
up1 <- res$coef[1]+1.96*res$coef[2]
> lo1;up1
[1] -1.684872
[1] 22.9317
この信頼区間だと，下限がマイナスの値になってしまっていて，都合が悪いです。
ですので，誤差に対数正規分布を仮定して解析をやり直したものが下記になります。
#y というデータがあるもとで，対数正規分布を誤差分布に仮定したときの平均密度（元のスケール）の95%信頼区間
res2 <- summary(lm(log(y)~1))
lo2 <- exp(res2$coef[1]-1.96*res2$coef[2])
up2 <- exp(res2$coef[1]+1.96*res2$coef[2])
> lo2;up2
[1] 1.345521
[1] 9.318763
私がやりたいこととしては，yを使わず，res$coef[1]とres$coef[2]を使って，lo2とup2を出せないか，ということを考えています。
厳密にlo2とup2に一致しなくとも，正規分布を仮定したときの平均値と標準誤差しかないが，それらから対数正規分布を仮定したときの信頼区間を出す方法，を探しています。
乱数種を替えればyが変わるのでlo1とup1も毎回変わりますが，マイナスのlo1が出る仕組み自体がまずいのでに，それを正の値しか出ない方法に変えたいのです。

どうぞよろしくお願いいたします。