ぎょーむ日誌 2009-06-19

2009 年 06 月 19 日 (金)

• 0700 起床． 朝飯． コーヒー． 0800 自宅発． 曇． 0810 研究室着．
• Paul に「これおもしろいよ」とススめられた論文を読んでみる．

  Lele SR, Dennis B, and Lutscher F. 2007. Data cloning: easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods. Ecology Letters 10:551-563 (url)

  これ， 以前にちらっとながめたときには， 「EM アルゴリズム の再発明にすぎないんでは?」 などと思っていたのだが …… くわしく読んでみたらちょっとちがった．
• EM アルゴリズムとの最大のちがいは， 数値的な最尤推定法がまったく必要ないこと． ただ単に k 回にわけて事後分布をサンプリングして， それぞれの平均をとったり， そのまた平均 (標本数は同じなので単なる全体の平均) をとれば OK というもの …… Appendix にこの証明らしきモノがあるんだけど， ちょっと意味不明ぎみなんだよね (私には)．

  Because the mean of a degenerate distribution is the point at which it is degenerate, the mean of the posterior distribution for large enough k approaches MLE of θ.

  小わけしたサンプルの平均値の分布は正規分布に漸近しそうだし， また事後確率密度関数を何重にも重ねていけば最尤推定値 (MAP) 周辺でツブれていきそうだけど …… 両者は別モノなのでは?
• ついでながら， このような推定値 (平均値だの平均値の分散だの) を得たいのであれば， k 個の「データクローン」なるものを生成する必要はなく， ふつーに MCMC 計算していって長いサンプル列をえて， それを k 個に等分割すればよいだけ， という気がする．
• といったことを Paul にメイルしてみる．
• で， その後になって上の理解はあまり正しくない， と気づいた． まず最初に，
  • 一度の k 通りのデータ (コピー) にもとづいてサンプリングする
  • データコピーしないで長くサンプルする
  の事後分布からのサンプル平均などが異なるのであれば (事後分布の定義は形式的には異なっているわけだが， 実際のサンプリング結果も異なるのか?)， 上記引用の退化分布うんぬんは 「そんなものかもなぁ」 といった気がしないでもない …… といったあたりからアタマの中での再検討がはじまり ……
• あ， この論文の Figure 1 ってのが退化分布 (degenerate distribution) ぢゃん． われながらうかつ． データを 100 とか 200 回コピーすると， 事後分布は最尤推定値のまわりに「退化」するみたいだね．
• まあ， この方式で データのコピーとかとると MCMC 計算の収束とかすごく遅くなるかもね． コピー数は 100 とか 200 ぐらい必要なかんじなので， 単純に計算時間は 100 倍・200 倍になる． さらに「パラメーターがうごきにくくなる」 のであれば …… で， そこまで時間かけて最尤推定した場合のご利益としては， AIC とかが使えるといったことぐらいかな．
• 工藤さんからの連絡で， 工学部疎開生活の電力問題はかなりキビしいものになりそう， とのこと． 院生部屋にまわす電流は 170A になるとか， 冬季にこの部屋を加熱するためにガスストーブが必要になりそうだから (ちなみにガスストーブは疎開時にすべて廃棄される予定だった)， ガス配管を今からやってみるとか …… いったいどうなるんだろうね．

• データ解析こんさる． また (呪われぎみな) library(lme4) の mcmcsamp() を調べてみる． まず， どうでもいいことだが， 図示は xyplot() だの densityplot() だのでできるのね (cf. 昨日)．
• で， mcmcsamp() は GLMM なときにうまく動かないかも， とのことなので (ほとんど使ったことないのに) library(MCMCglmm) はどうでしょうかね， などとてきとーな返信をしてみる．
• なかまさんから R 作図における日本語文字のあつかいに関して， いろいろ教えていただいたの返信メイルかき．
• 疎開先での院生部屋のネットワークがたいへんなことになりそうなので， いっそうの無線 LAN 化を推進してみることにする． 雪野さんに予算その他いろいろ教えていただいてから， 北大生協に連絡メイル．
• 面談式のデータ解析こんさる．
• A802 院生部屋の CANON LBP-5400 がよーやくにして修理される …… といっても Y トナーカートリッジを交換するだけなのだが． 2 週間ほど放置されてた．
• Paul (もうアメリカにかえるそうで) とちょっと議論． くだんの data cloning 論文の appendix， 「k ブロックごとの平均値の平均値は最尤推定値」 まわりはそんなにきちんと証明されてない， そして 「事前分布にあまり影響されない」 うんぬんはそんなに厳密には示されてないんでは， といったハナシなんかを．
• めいる式のデータ解析こんさる． CRAN ではなく R-forge に入ってる MuMIn (むーみん? Development Status: 2 - Pre-Alpha) なる package， 私などにはかなり邪悪なものに思える …… というのも， 例のモデル選択カエル本にべったりな packge で， QAIC だの Akaike weight だのをむやみに計算するしろもので．
• で， いろいろ解読してみたのだけど， まだまだ開発途上のもの， とわかったのでそのあたりメイルで送信．
• car.normal() なやりとりのつづき． GeoBUGS manual をよく読んでみると，

  Since the CAR model defined above is improper (the overall mean of the S_i is not defined), it can only be used as a prior distribution for spatially distributed random effects, and not as a likelihood for data. It is often convenient to assume that suchrandom effects have zero mean. Besag and Kooperberg (1995) show that constraining the random effects to sum to zero and specifying a separate intercept term with a location invariant Uniform(-infty, +infty) prior is equivalent to the unconstrained parameterisation with no separate intercept. WinBUGS 1.4 includes a distribution called dflat() which corresponds to an improper (flat) prior on the whole real line. This prior must always be used for the intercept term in a model including CAR random effects.

  WinBUGS-1.4 には dflat() なる improper な「平らな」事前分布があるから， car.normal() を random effects にするような一般化線形なモデルの 「切片」 にはそれを使え， といったもの．
• え， とおどろいで， car.normal() の例題 をやりなおしてみたけど， 事後分布はとくに変わらなかった． ひとあんしん．
• 冷静になって考えてみれば， dflat() なる事前分布 (一様分布みたいなもの) は (例題の解説文などで使っていた) dnorm(0.0, 1.0E-4) とかと本質的に異なるわけではないな． つまり， improper な事前分布ってのは 「ハシからハシまで積分しても 1 にならない，尤度が評価されない」 ぐらい意味で， なにか抜本的に異なるわけではない． dflat() の実装は， おそらくシステムの許容する最大・最小の実数を範囲とする一様乱数， 対数尤度はつねに 0， といったところではないかな?
• 本日の A 棟耐震工事疎開準備のかたづけ …… やらなかった． なんだかばたばたした夕方だったので．
• 1900 研究室発． 1915 帰宅． ぐったり． 晩飯の準備． 晩飯．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0840): 米麦 0.6 合． ネギ・卵の炒めもの．
  • 昼 (1330): 研究室． 米麦 0.8 合． 煮干． リンゴ．
  • 晩 (2140): 米麦 0.8 合． 麻婆豆腐． ダイコン・ニンジン・タコの煮物． ヨーグルト．