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

No.22639　その2 正規分布を仮定した信頼区間を。。。　　【青木繁伸】　2018/12/01(Sat) 10:59

いろいろ枝葉末節にこだわりすぎた感がありました。
サンプルサイズが小さい場合に精度が低いのはやむを得ない（当然）なので，結局は，対数正規分布のパラメータとそれを対数変換した正規分布のパラメータの関係式から解を求めるしかないようですね。
対数正規分布の標本平均，標本分散から，非線形連立方程式で正規分布のパラメータを求め，それに基づき信頼限界を計算し，それを指数変換する。

元データなしで，信頼限界を求める

getCL = function(MEAN, VAR, n) {
  func = function(p) {
    return(c(exp(p[1]+p[2]/2)-MEAN, exp(2*p[1]+p[2])*(exp(p[2])-1)-VAR))
  }
  library(nleqslv)
  ans = nleqslv(c(2,2), func) # 対数正規分布の平均値と分散から
  mean = ans$x[1] # 対数をとったときの正規分布の平均値と
  sd = sqrt(ans$x[2]) # 標準偏差を非線形連立方程式を解いて求める
  cat(sprintf("mean = %f, sd = %g\n", mean, sd))
  cl.normal = mean+c(1, -1)*qt(0.025, n-1)*sd/sqrt(n-1) # 正規分布目盛りでの信頼限界
  cl = exp(cl.normal) # 対数目盛りに変換
  cat(cl)
}

対数正規分布にしたがうデータから，標本平均・標本分散　MEAN = 3.390011 VAR = 37.060
が求められている（元データはない）

これから対数をとった正規分布の標本平均，標本標準偏差を求める
mean = 0.500339, sd = 1.20041
となる
それに基づいて，信頼区間を求めると
0.667087 4.077615

サンプルサイズが小さいと rlnorm で生成されたデータの平均値と分散は，母数とはかなりかけ離れたものになっている

No.22642　Re: その2 正規分布を仮定した信頼区間を。。。　　【Saito】　2018/12/03(Mon) 10:29

コメントありがとうございます。
しかし，一つ前のスレッドでも書きましたが，こちらの手元にあるのは平均値とその標準「誤差」であり，元データの平均値と標準「偏差」ではありませんでした。そのため，微妙に私のやりたいこととずれていることに気づきましたので，再度ご教授頂ければ幸いです。何度も何度も申し訳ありません。

まず，元データの平均値と標準偏差が使える場合に合うことを示します。

> set.seed(1)
> n <- 100000
> u <- 1
> sig <- 0.5
>
> #生物の密度データがあるとする（非負）。平均密度とその信頼区間を知りたい。
> 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]
>
> 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])
>
> MEAN = mean(y)
> VAR = var(y)
> lo2;up2
[1] 2.706802
[1] 2.723691
> getCL(MEAN, VAR, n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906>

上記のように，MEAN=mean(y), VAR=var(y)として，元データyの平均値と分散が使えるという状況であれば，getCLとlo2, up2がほぼ一致することを確認しました（なお，sigが大きいと合わないようなので小さくし，小サンプルであることの問題も省きたいため，サンプルサイズnは大きくしてあります）

ところが，私の手元にあるのは，yの平均値と標準偏差ではなく，解析結果であるyの平均値とその標準誤差となっています。つまり，res$coef[1]とres$coef[2]が手元にあるということでして，mean(y)とvar(y)の値が手元にあるわけではありませんでした・・・。私の説明が拙かったようで，本当に申し訳ありません。
確認してみましたが，当然ながら，こちらとは合いません。

> MEAN = res$coef[1]
> VAR = res$coef[2]^2
> lo2;up2
[1] 2.706802
[1] 2.723691
> getCL(MEAN, VAR, n)
mean = 1.124605, sd = 0.0016874
3.07897 3.079034>

何度もお尋ねするのは心苦しいのですが，res$coef[1]とres$coef[2]が手元にあり，mean(y)とvar(y)の値（とn）が手元にない条件で，lo2, up2に準じるものは出せるでしょうか。何度も申し訳ありません。

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

No.22643　Re: その2 正規分布を仮定した信頼区間を。。。　　【青木繁伸】　2018/12/03(Mon) 11:23

> 私の手元にあるのは，yの平均値と標準偏差ではなく，解析結果であるyの平均値とその標準誤差となっています

標準誤差=標準偏差/sqrt(サンプルサイズ）　ですけど？
標準偏差=標準誤差*sqrt(サンプルサイズ）
分散=標準偏差^2

でしょう？

> res <- summary(lm(y~1))
> lo1 <- res$coef[1]-1.96*res$coef[2]
> up1 <- res$coef[1]+1.96*res$coef[2]
lm する必要なんてなくて，
res$coef[1] = mean(y)
res$coef[2] = sd(y)/sqrt(n)
ですよね。

簡単なところを複雑にやっているので，頭が混乱します

> set.seed(1)
> n <- 100000
> u <- 1
> sig <- 0.5
> #生物の密度データがあるとする（非負）。平均密度とその信頼区間を知りたい。
> 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
[1] 3.068823
> up1
[1] 3.089189
> res$coef[1]
[1] 3.079006
> mean(y)
[1] 3.079006
> res$coef[2]
[1] 0.005195524
> sd(y)/sqrt(n)
[1] 0.005195524
> (res$coef[2]*sqrt(n))^2
[1] 2.699347
> var(y)
[1] 2.699347
>

No.22644　Re: その2 正規分布を仮定した信頼区間を。。。　　【青木繁伸】　2018/12/03(Mon) 11:39

つづき

> mean(y)
[1] 3.079006
> var(y)
[1] 2.699347
> n
[1] 1e+05

以下で mean, var を使っているけど yは必要ない。平均値と分散を引数として与えているだけ。

> getCL(mean(y), var(y), n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906

以下のように必要なのは数値だけ
ただ，ここで示したように有効数字7桁では正確ではないので，コンピュータ上で持っている正確な値を使っただけ。
あなたは，y は持っていないが，計算された平均値と分散（をサンプルサイズで割ったものの平方根である標準偏差）を持っているので，その二つの数値を指定するだけ。

> getCL(3.079006, 2.699347, 100000)
mean = 0.999331, sd = 0.500551
2.70805 2.724906

標準誤差と分散の関係をちゃんとすれば，res$coef[1], res$coef[2] を使っても，当然同じ結果になる

> res$coef[1]
[1] 3.079006
> (res$coef[2]*sqrt(n))^2
[1] 2.699347

> getCL(res$coef[1], (res$coef[2]*sqrt(n))^2, n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906　当然ながら同じになる

No.22647　Re: その2 正規分布を仮定した信頼区間を。。。　　【青木繁伸】　2018/12/03(Mon) 13:49

また枝葉末節なこと（基本的なこと）になりますが，あなたが使っている 1.96 は標準正規分布由来のものだが，本来は t 分布に基づく値でなければならない。

http://aoki2.si.gunma-u.ac.jp/lecture/Average/Mean2.html
の「母分散が不明の場合」にあたるので。

いまはテストデータでやっているから n がとても大きく，1.96 で問題ないが，実際に運用するときには問題になる。

No.22650　Re: その2 正規分布を仮定した信頼区間を。。。　　【Saito】　2018/12/04(Tue) 08:56

コメントありがとうございます。
不勉強で申し訳ありません。理解できていないのですが，標準誤差から標準偏差に戻すときに，nが未知でも戻るのでしょうか？今手元にあるのは，平均値とその標準誤差のみで，サンプルサイズnは不明なのです。getCLの中身を見る限り，何だか無理なような気がしてきています・・・。
もし戻らないとなると，当初の目的を果たすことは不可能でしょうか。

t分布については失念していました。
n=10だと，2.228ですね。ありがとうございます。

No.22652　Re: その2 正規分布を仮定した信頼区間を。。。　　【青木繁伸】　2018/12/04(Tue) 09:38

> 標準誤差から標準偏差に戻すときに，nが未知でも戻るのでしょうか？

定義式に書かれているとおり。
nが分からないなら無理でしょう。
現代統計学では，サンプルサイズこそ一番大事なものです。

妥当と思われるnを仮定するか，何通りかのnで代替してこの範囲という結果を示すしかないでしょう。

No.22653　Re: その2 正規分布を仮定した信頼区間を。。。　　【Saito】　2018/12/04(Tue) 11:25

コメントありがとうございます。
nがわからなければ不可能という旨，承知いたしました。
色々とご教授くださりありがとうございました。

以下余談です。
何故こうなるのかわからず，おそらく間違っていると思って相談しませんでしたが，この件について私ではない方（前任者）のメモで

Lower=exp(log(res$coef[1])-1.96*res$coef[2]/res$coef[1])
Upper=exp(log(res$coef[1])+1.96*res$coef[2]/res$coef[1])

として，平均値を対数を取って，なぜかCVを足す（引く）というやり方で正規分布を仮定して得られた平均値と標準誤差（res$coef[1], res$coef[2]）を使って，対数正規分布を仮定したときの95%信頼区間(Lower, Upper)になるというメモがありました。やってみても，getCLと全然違う範囲になりますし，CVという単位が消えた値を足すのも奇妙なので，間違っているのでは，と思っています。とはいえ，他の用途でも構いませんので，この式に見覚えがありましたら，ご教授くだされば幸いです。