ぎょーむ日誌 2006-11-23

2006 年 11 月 23 日 (木)

• 0800 起床． 朝飯． コーヒー． あ， 今日は祝日か， ということで洗濯． 怠業してる場合ではないんだが …… と思いつつも洗濯機がまわると怠業もーどになってしまう私．
• 怠業ついでに運動だ． 1140 自宅発北大構内走． 曇天． 路上のあちこちに雪や水たまりがのこってて走りにくい． そしていつもの往復区間である 低温研-獣医学部前の歩道は完全に雪に被覆されている． 1225 帰宅． 昼飯．
• 1310 自宅発． 曇． 1325 研究室着．
• さーて， 昨晩かえりぎわに 「あとはもうどうにでもなれ， とりあえず全種・全 stage で計算しろ」 と命じておいた MCMC 計算だが …… 小山さんから事後分布作図の R プログラムを拝領して図にしてみたところ， やっぱり収束してないとわかった． モデリングのいいかげんな状況での なげやり計算はダメですな． だめモデルの挙動の把握には役たつけど．

[事後分布の作図例]

まだぜんぜん信用できない計算結果です．

• たしかにうまくいってない計算も多い …… 一種ずつ取りだして図にしてみると， しかしここしばらく MCMC 計算の実験対象に使っている 植物コード Solidago s (アキノキリンソウ実生) はどのパラメーターもきれいに収束してやがるな． 標本数が多いせいか? 上の図でもミカヅキグサとかはひと山に見える． ふーむ ……
• お茶部屋会話．
  「久保さん， そこにおみやげのお菓子が」
  「ありがとうございます． しかしお茶部屋のお菓子を取るのは勇気がいるので， あとで誰も見てないときにこっそりといただきます」
• サロベツ ZIP 階層ベイズモデル， 無謀なモデル改造の一環としてブロック (line) ごとの random effects を多変量正規分布 dmnorm(平均, 逆数分散共分散行列) からださせるようにする． 今までは完全に独立した正規分布から 生成させるようにしていたんだけど， 収束の悪さはどうもこのあたりに何か原因があるような気がしてならない．
• さて， めんどうなのはこの新しい事前分布 (prior) たる多変量正規分布 の超事前分布 (hyper prior) だ． 単なる正規分布の場合だと， その分散 (の逆数) の事前分布はガンマ分布を「なんとなく」 というか， 共役事前分布 (conjugate prior) だからという理由で使ってきた． では， 多変量正規分布の分散共分散行列 (の逆数) の共役事前分布は何か? それは Wishart 分布 なんだよね …… ようこそ， アタマがヘンになりそうな多変量確率分布の世界へ．
• R の library(bayesm) に入ってる rwishart() や library(MCMCpack) に入ってる rwish()， すなわち Wishart 乱数はきだす関数でしばらく遊んでココロをおちつける． で， コードの改造にとりくんだわけだが …… 無情報 Wishart 分布にするためのパラメーターの与えかたがわからん． 調べてみてもナゾ．
• えーと今回の多変量正規乱数の次元は 2 …… ということで， なンか対角行列になっていればよい正定値行列 R を R の中で R <- structure(c(1, 0, 0, 1), .Dim = c(2, 2)) として data として R2WinBUGS 経由でわたす． Wishart 分布の自由度は 2． われながらでたらめだな ……
• 600 MCMC step 計算， 472 秒で終了． 収束はわるくない． DIC = 1383． 問題の多変量正規乱数由来の事後分布もそこそこ値がゆれてる．
• よくわからんので， R <- structure(c(0.01, 0, 0, 0.01), .Dim = c(2, 2)) としてみる． 482 秒． 「逆数」分散共分散行列の値はむちゃくちゃでかい， ということは分散とかすごく小さいとなるわけだが …… そのあおりをうけて， なのか DIC = 1293．
• まだよくわからんので R <- structure(c(100, 0, 0, 100), .Dim = c(2, 2))． とうぜんこれは上とは逆で 「逆数」分散共分散行列の値はむちゃくちゃ小さい． なぜか deviance 小さく DIC = 1284． すごくチェイン間差がある， つまりぜんぜん収束してない．
• なんとなく R <- structure(c(0.1, 0, 0, 0.1), .Dim = c(2, 2)) というかもっと簡潔に R <- diag(0.1, 2) でよいのかしらん， と勝手に決めて 1000 MCMC step 計算やらせてみる． n.burnin = 700． この問題に関しては 300 step のサンプリングでは足りないのかな?
• 786 秒． DIC = 1242． aP とか収束あまりよくない．

• R <- diag(0.01, 2) にしてみるか …… 810 秒， DIC = 1095． 収束はまだまだ．
• お， 伊東さん苦闘の 「Mac OS X上で、WinBUGS+R2WinBUGSを使用する」 が (とそのまわりの 経緯も) 公開されている …… CrossOver と Darwine のコンビネイションわざ， でしたか． うーむ， すごい． おつかれさまです．
• サンプル長かえてみる． n.iter = 1300, n.burnin = 700, n.thin = 2． …… 計算時間 1030 秒， DIC = 1294． たいへんよく混交している．

• 多変量正規分布の共分散の事後分布の平均値は負だった (かなり「ゆるい」けど) …… で， その解釈はどうなるわけ? 「植物の定着できる場所は限定されてる & いるところはたくさん」 VS 「どこでも定着できる & しかしどこでも低密度」 の区別がつかん， ということかしらん． そしてなぜ Gibbs sampler は無相関の正規分布ペアから 「負の相関をもつ事後分布」 を生成できないんだろうか …… いや， 生成してたのかな?
• よくわからんところもあるけど， ということで， 今晩の全種・全 stage 長時間計算は n.iter = 1600, n.burnin = 1000, n.thin = 2， 昨日とのちがいは多変量正規分布の事前分布 & Wishart 分布の超事前分布ということで． 1900 研究室発． 謎のお供をする晩飯． 2100 帰宅．
• データ解析こんさるメイル書きつづく． R FAQ． 「わくの外」に何かを描くには par(xpd = NA) が必要．
• [今日の運動]
  • 北大構内走 45 分間．
• [今日の食卓]
  • 朝 (0820): 米麦 0.7 合． 海藻スープ．
  • 昼 (1250): タラコスパゲッティー．
  • 晩 (1930): 北 14 西 1 のグロビュール． スープカレー．