ぎょーむ日誌 2005-03-08
2005 年 03 月 08 日 (火)
-
0700 起床．
朝飯．
コーヒー．
体温測定 36.7°C．
気分的にはまだすっきりしてないけど，
とりあえず平熱のよーで．
風邪薬の服用ヤメ．
0850 自宅発．
雪．
0905 研究室着．
今日は平常どーり 8F まで階段をあがってみたんだが，
まだ万全ではないよーですこし体温上昇が実感できた．
-
「非常階段にいたる廊下の幅は 1.6m 確保せよ，
じゃまなものは撤去」
なる責任のがれ oriented な命令に対応する雑用．
とりあえず撤去すべきロッカー上のよけーなものを撤去し，
またロッカーのうち二基は院生部屋にいれるので，
その退避先のかたづけ．
すると 8F 闇ネットの配線の一部をくみかえる必要があ，
……
と一時間ばかりを費してしまった．
-
アカマツ問題に対応するための，
parametric bootstrap 尤度比検定プログラムの改善．
まずは，
「平均は同じだけど，分散が異なる」
正規分布の最尤推定を
R
でやらせる方法を検討する．
まあ，
これは
optim()
つかって自分でやるしかないんだが．
-
などとプログラミングしてると，
数種類の雑用メイルが．
ついつい返信してしまうんだよなぁ．
-
「母分散は異なるけど平均は同じ，
と仮定する二つの正規分布の最尤推定」，
とりあえず，
関数はかけた．
こんなもんか?
まー，
関数内関数とか，
関数内スコープとか，
非負パラメーター推定のための
exp()
わざとか，
いろいろ使ってるわけですが．
estimate.gaussians.unequal.var <- function(sample1, sample2, ...)
{
logL <- function(p) # 対数尤度
{ # p == c(mean, var1, var2)
log.g <- function(sample, mean, var.log)
(
- 0.5 * log(2 * pi * exp(var.log))
- (sample - mean)**2 / (2 * exp(var.log))
)
sum(log.g(sample1, p[1], p[2])) + sum(log.g(sample2, p[1], p[3]))
}
……
(以下，うだうだと書いていたんだけど，
本日の夕方ごろにはすごくややこしくなってしまったので中略)
……
}
しかし，
きわめて標本数の少ない呪われ割りざん値
(むろん，というか残念ながらというか 連続値 / 連続値)
の「検定」のために，
「母分散の異なる正規分布で近似して，
parametric bootstrap 尤度比検定」なる
すとらてぢーは手がこんでるわりには
阿呆っぽいような気もしてきた．
ホントに標本数が少ない場合でしか使えんよーな
(今回はまさにそれなんだけど)．
本日も人口密度ひくく，
お茶部屋雑談にとらっぷされてしまう．
高校地理でならった「カナダの石油はオイルサンド (砂まじりの油)」
というような雑学を思い出したり．
あとで調べてわかったんだけど，
その埋蔵量はじつに莫大でサウジアラビア全体より大きいらしい
……
しかも精製コスト削減にも成功しつつあるらしい．
石井さん・長谷川さんに手伝ってもらって，
廊下のロッカー移動．
昼飯．
かとーさんがインフルエンザから回復してひさしぶりに登校してこられた．
なんだか竹中さんと大阪大会の次の大会に関連する
雑用に関連するメイルのやりとり．
なぜわれわれがかかるコトにまきこまれているのだろう?
生態学会といえば，
大阪大会では parasitoid のハナシをしないといけない
……
ヒトたちのための「うらかた」をしなければならないのですが．
あ，
自由集会
の準備や勉強もいろいろやりたいんだけど．
なんかいろいろ用件がくると，
私の resource allocator はいともたやすく破綻するよーです．
そういや，
低温研の加藤さんからはダケカンバ葉群推定概念図を発注されてたっけ
……
時刻は 1600．
苫小牧 parasitoid 問題はもちろん
アカマツな仕事も進捗しておらず，
雑用も処理できていない．
夜までがんばって，
parametric bootstrap 尤度比検定
(分散の異なる正規分布母集団からの標本の比較)
の改良．
いや，
改良してるつもりでいろいろとぢたばたやったんだけど，
あまり改良されたかんぢはしない．
たぶん，
単純まちがいのたぐいはすでに除去できてる，
と思うんだけど．
なにが気に入らないか，
というと
……
まさにその状況に適合した例題を作ってやっているのに，
帰無仮説をきちんと棄却してくれない場合が多々ある，
そしてうまくいく事例では Wilcoxon 検定でもうまく棄却できる，
ということだ．
# --- 例題 1 ---
> x <- rnorm(20, mean = 0, sd = 2); y <- rnorm(20, mean = 1, sd = 2.2)
> wilcox.test(x, y)
Wilcoxon rank sum test
data: x and y
W = 197, p-value = 0.9467
alternative hypothesis: true mu is not equal to 0
> r <- pb.lrt.gaussian(x, y); r$p
# gaussian(c(x, y)): logLik(heteroscedasticity) = -89.5, logLik(homoscedasticity) = -89.8
# mean.xy = 0.446592, sd.x = 2.28671, sd.y = 2.28671
[1] 0.652
# --- 例題 2 (設定は 1 と同じ; 乱数セットが異なるだけ) ---
> x <- rnorm(20, mean = 0, sd = 2); y <- rnorm(20, mean = 1, sd = 2.2)
> wilcox.test(x, y)
Wilcoxon rank sum test
data: x and y
W = 121, p-value = 0.03264
alternative hypothesis: true mu is not equal to 0
> r <- pb.lrt.gaussian(x, y); r$p
# gaussian(c(x, y)): logLik(heteroscedasticity) = -89.0, logLik(homoscedasticity) = -91.1
# mean.xy = 0.193208, sd.x = 1.728, sd.y = 2.90467
[1] 0.032
で，
かんぢんのアカマツデータセットだが
……
なにしろ 4 標本 vs 4 標本というさらにキビしい状況なんで，
分散の推定がむちゃくちゃになる，
といいますか．
# こんだけしか標本ないんですよ……
> x <- d$AmaxW[1:4]; y <- d$AmaxW[5:8]
> sort(x)
[1] 0.0342 0.0388 0.0563 0.0638
> sort(y)
[1] 0.0653 0.0717 0.0783 0.0931
# parametric bootstrap 尤度比検定
> r <- pb.lrt.gaussian(x, y); r$p
# gaussian(c(x, y)): logLik(heteroscedasticity) = 21.0, logLik(homoscedasticity) = 20.7
# mean.xy = 0.0731336, sd.x = 0.0276698, sd.y = 0.0110553
[1] 0.108
……
うーむ，
上の logLik()
(最大化対数尤度)
と「平均ちがう and (分散おなじ or ちがう)」
なモデルのそれを比較してみると
……
> logLik(lm(x ~ 1)); logLik(lm(y ~ 1))
`log Lik.' 11.956 (df=2)
`log Lik.' 12.620 (df=2)
> f <- as.factor(c(rep("T", 4), rep("C", 4)))
> logLik(lm(c(x, y) ~ f))
`log Lik.' 24.521 (df=3)
「平均は異なるけど分散はおなじ」
というモデルが最良，
か．
いやはや．
しかしながら，
この程度の対数尤度差に関する尤度比検定では上のとーりになるわけで．
まあ，
矛盾はしていないんだけどね
(検定とモデル選択は別ものだから)．
なぜ尤度比検定がここまで劣悪なのか，
よくわからん．
このような「重なりのない二標本」で
「分散の異なるけど平均の等しい母集団から得られた標本」
という仮定がナンセンスなのか?
順序統計量つかった Wilcoxon test だと「ゆーい」になっちまうんだけど，
だからといって
……
> wilcox.test(x, y)
Wilcoxon rank sum test
data: x and y
W = 0, p-value = 0.02857
これがこの標本セットにたいして妥当な検定なのかどうか，
かなり疑問だ．
なにしろ合計 8 標本とかだと，
p 値のとりうる値もすごく限定されているわけで．
たとえば，
x の最大値と y の最小値の差はそれほど大きくないんだけど，
これらを入れ換えると，
> x <- sort(x)
> y <- sort(y)
> rbind(x, y)
[,1] [,2] [,3] [,4]
x 0.0342 0.0388 0.0563 0.0638
y 0.0653 0.0717 0.0783 0.0931
> wilcox.test(c(x[1:3], y[1]), c(x[4], y[2:4]))
Wilcoxon rank sum test
data: c(x[1:3], y[1]) and c(x[4], y[2:4])
W = 1, p-value = 0.05714
p 値はちょうど 2 倍に増えるわけで．
これらの中間の値は絶対にとりえない．
(正規分布でなくてもいいけど)
何かパラメトリックな推定なら
x[4] (x の最大値)
と
y[1] (y の最小値)
なんて variance ですぐに埋められるような差にすぎない
……
しかし，
それは今回つくった parametric bootstrap 尤度比検定がへっぽこである，
という説明にはなっていない，
よねえ
……
なんかナットクできんなぁ，
といぶしみつつ撤退
……
2245 研究室発．
吹雪．
2300 帰宅．
ふう．
体温 36.2°C．
体温計はもう片づけてよさそうだ．
体重 73.8kg．
晩飯．
要約すれば，
統計学的検定なる操作をするときに，
良い帰無仮説を作ることが重要なのである
(モデル選択だと，
モデルの適切な階層性を設計すること)．
で，
何が良い帰無仮説なのか，
統計学だけでは決められない部分が存在しうる，
といったあたりがイヤなところだ．
といいますか，
統計学ってのは人間が適切なる統計学的モデル群を準備しているときに，
その状況のもとでどのような計算をなすのが妥当であるか，
を明らかにするものなんで
……
しかしながら，
計算してみないことには「適切なる」だったかどうかわからん，
って場合が多々あって
……
[今日の運動]
[今日の食卓]
- 朝 (0730):
リンゴ．
ヨーグルト．
- 昼 (1400):
研究室お茶部屋．
米麦 0.7 合．
ハクサイ・タマネギ・エノキダケ・豆腐のカレー．
- 晩 (2400):
米麦 0.7 合．
ハクサイ・タマネギ・エノキダケ・豆腐のカレー．