ぎょーむ日誌 2005-07-08

2005 年 07 月 08 日 (金)

• 0800 起床． 朝飯． コーヒー． また昨日の Markov Chain Monte Carlo (MCMC) 植生パターン問題にとらっぷされてしまう． 0920 自宅発． 晴． 0935 研究室着．
• この日は結局一日中 「植生分布 MCMC」論文たちというか 空間統計学勉強に没頭してしまった．
• MCMC な問題ってのは， ようするに植物の分布で言えば 「場所 A とその近傍の場所 B での植物 X のいる・いないは独立ではない (相関がある)」 という状況を Gibbs 分布であつかう， ということだ．
  「場所 A にいるんだから B にいても不思議なかろう」
  「いやいや B にいるから A にもいるんだ」
  と何が原因で結果なのか判然としなくなったときであっても， A と B のあいだの「つながりの強さ」とかを なんとかうまく推定しましょう， というもの．
• やはり， というか， じつは植生 MCMC 論文なるものはそれほどあるわけではなさそう． おそらくその最初の論文は Wu & Huffer たちによるもので
  • Wu, H and Huffer, F.W. 1997. Modelling the distribution of plant species using the autologistic regression model. Environmental and Ecological Statistics 4:49-64.
  • Huffer, F.W. and Wu, H. 1998. Markov Chain Monte Carlo for autologistic regression models with application to the distribution of plant species. Biometrics 54:509-524.
  どちらも 1995 年に投稿された論文だ． このあたり Wu の 1994 年の Florida State University での学位論文にもとづいており， たぶん Huffer はその指導教官なのだろう． どちらも統計学部なヒトたちである．
• で， やってることは …… Maximum pseudo likelihood 推定値 (MPLE) の MCMC による補正， つまり
  • Geyer, C.J. and Thompson, E.A. 1992. Constrained Monte Carlo maximum likelihood for dependent data. Journal of the Royal Statistical Society, Series B 54:657-699.
  など Geyer たちの 1990 年代前半の研究によっている． この計算方法は， 私がよくひきあいにだす
  • Hubbell, S., J. Ahumada, R. Condit and R. Foster. 2001. Local neighborhood effects on long-term surivival of individual trees in a neotropical forest. Ecological Research 16:859-875.
  での死亡率推定と同じである (Hubbell たちも Geyer & Thompson (1992) を引用して「こういう補正計算でいいんだ」と正当化している)．
• MCMC による MPLE の補正 (近似最尤推定値を得る)， という計算に関しては Ogata & Tanemura の 1980 年代前半の計算その他とかもあるんだが， なぜかこちらはあまり参照されないよーで (マルコフ確率場ではなく点過程をあつかっていたためか?)． 昨年末 に苦労して解読した (と言いますか， これが Gibbs な統計学勉強の出発点だったわけだが)， 間瀬さんの林冠動態論文では Huang & Ogata (1999) の MPLE 補正した推定量 (approximated maximum likelihood estimator; AMLE) を計算している．
• ともかく， このへん 1970 年代の Besag のわくぐみ， 加えて MPLE 補正計算が 1990 年代なかごろの主軸になっていたよーだ (今年 1 月ごろの マルコフ確率場 (Markov Random Field; MRF) の検討も参照)． つまり Wu & Huffer 論文たちで整理・比較・検討されてるよーに，
  1. coding 法 (Besag)
  2. pseudo likelihood 法 (MPLE 計算; Besag)
  3. MCMC を使った方法
  • MPLE 補正 (Geyer & Thompson)
  となってるわけだ． 詳しくは 間瀬 ・武田 (2001) 「空間データモデリング」 を参照．
• さて， 私の理解ではこの最後の「MCMC を使った方法」 として現在あげられるものとしては
  • MPLE 補正 (Geyer & Thompson, 尾形たち)
  • EM アルゴリズムを使った最尤推定
  • ベイズ推定
  となってるんではなかろーか． 最後のふたつはどちらも「隠れ」変数 (潜在変数) をうまく使った hidden Markov model (HMM) で解決せんといかんのでは， と考えて観測データに対応する尤度方程式をひねくっているわけだ．
• というのも， MPLE 補正は見かけは簡単だけど， 一種の近似計算なんでまあちょっと不安がある， と言いますか． さらに， MPLE 補正で何でもかんでも推定すりゃあいいってもんでもあるまい． この方式は植生分布みたいにとても「ゆっくり」 変化する現象にはむいているだろう． つまり観測データとして一枚の「地図」しか得られないような場合だ．
• …… と書いて「しかし MPLE 補正なんてへっぽこですよ」 と続けようとして， 「地図」が二枚以上ある場合でも， MPLE 補正計算はできるなぁと気づいてしまった． えーい． (後記: 翌日になって冷静に考えてみると， やはり「隠れ」変数つかわん pseudo likelihood 計算だとヘンになるかもしれん)
• ともあれこのあたり， 空間がらみの問題に限定してもなお， 生態学まわりのいろいろな問題に連関するわけで． たとえば， 生態学会大会の 「保全」 session の発表でよくみかける 「分布の推定」とか．
• あとで調べてみたら， そういう研究やはりすでにあった．
  • Argáez, J.A., Christen, J.A., Nakamura, M and Soberón, J. 2005. Prediction of potential areas of species distributions based on presence-only data. Environmental and Ecological Statistics 12:27-44.
  など． 上記のはいる・いないデータから空間相関あり potential 関数を推定， というやりくちで私の考えてる方式に近いよーな気がする． MCMC 法による Bayesian inference．
• そして， MCMC 法の応用は空間統計学に限定されるものではない， と強調しておきたい．
• あらら， もう 1420 すぎか． 昼飯．
• 「助手各位」 なる紙きれがまわってくる． 北大でも「助手」なる職位を 「助教」 (assistant professor)と 「助手」 (research associcate) にわけることになるそーで． 粕谷さん が以前に説明されてた 「助手が助教に」 というやつで， 「すでに「助教」の任期制が」 にあるがごとく， 北大では (その時期は不明ながら) 「新制度移行時」に採用される 助教は最大 5 年 + 再任 2 年という任期つきになるそーで (助手は教育研究支援職員とくっつけて消滅していく)． じゃ， 私みたいなのはどうなるの， ということなんだけど新制度移行時には「助教」 「助手」のどちらかに名前がかわって， その後はどうなるか書いてないわけで．
• なぜか昨日・今日と別の院生たち (赤坂君・森君) に同じ R 質問された． ってことは， 同様の疑問を持ってる「隠れ」(潜在) 院生も多かろう と最尤推定されるので， ここにその回答を書いてみる． たとえば， こういうデータ d

  	> d
  	   spc treatment value
  	1    A         x     3
  	2    B         y     6
  	3    A         z     6
  	4    B         x     2
  	5    A         y     5
  	6    B         z     6
  	7    A         x     3
  	8    B         y    11
  	9    A         z     7
  	10   B         x     2
  	11   A         y     3
  	12   B         z     4
  	13   A         x     3
  	14   B         y     7
  	15   A         z     4
  	16   B         x     7
  	17   A         y     2
  	18   B         z     5
  

  があったとする (R 勉強にはこういう例をさっと作れることがすごく重要で， たとえばこの場合だと d <- data.frame("spc" = rep(c("A", "B"), 9), "treatment" = rep(c("x", "y", "z"), 6), "value" = rpois(18, 5)) などとすればよい)． これの spc かつ treatment ごとに value の平均値を計算したければ

  	> tapply(d$value, d[, c("spc", "treatment")], mean)
  	   treatment
  	spc      x      y      z
  	  A 3.0000 3.3333 5.6667
  	  B 3.6667 8.0000 5.0000
  

  というふうに tapply() すればよい．
• また空間統計学わーるどにはまりこんでしまう．
• てゆーか， 以前に考えてた logistic モデルと Gibbs 分布の対応がわからなくなってしまったんだよね ……
• この問題は帰宅してからの深夜の北大構内走の最中にとけた． ちょっととりっきーかもしれん． というのも， logistic モデルの側から考えるのではなく， 「はじめに Gibbs 分布ありき」 の方向から検討したほうがわかりやすいように思えるんで． マルコフ確率場 (MRF) のある配置 k での各点の値を とする． したがって配置 (configuration) k とは といったパターンのことだ． この配置が出現する確率は Gibbs 分布で与えられ，

ここで は各場所 i ごとの変数で (あ vector だから太字にすべきだったか)， は隣接した場所 j との相互作用の強さをあらわす． 分母はかの悪名たかき 「すべての可能な配置に関する和」 となっている．
• さて， ここで (発想の飛躍というか，とーとつながら) この Gibbs 分布の Metropolis-Hastings sampler を考えてみよう …… いや， ここは Gibbs sampler を考えるべきだな． ある配置 k のもとで場所 i の値が y となる条件つき確率は

と書ける． どうです? だんだん logistic っぽくなってきたでしょ (いうまでもなく exp(0) = 1)． そしてすべての場所に関する上の確率の積は pseudo likelihood 方程式になっていて， そこから MPLE が得られるわけだ． なるほど．
• では， ここでもし場所間での相互作用がなければ? これは とすればよく， そうすると配置 k とは無関係となり，

と logistic model が得られるわけだ． うーむ， なんというつながりかた …… たぶんこのあたりの専門家にとっては 「直観でわかるハナシ」 だから論文とか教科書にはのってないんだよね．
• しかし大学にいるあいだは上の関係が理解できずに撤退． 1850 研究室発． 1900 帰宅． 北大構内走にでるまでうだうだと逡巡． 2240 自宅発北大構内走． 上述の Gibbs 分布 ⇔ logistic 問題はここで解けた． 金曜日の夜は北大構内走者が多いのか? 2335 帰宅． 体重 73.6kg． 晩飯．
• [今日の運動]
  • 北大構内走 55 分間．
• [今日の食卓]
  • 朝 (0900): 米麦 0.7 合． タマネギ・ニラ・ショウガ・ホタテの炒めもの．
  • 昼 (1430): 研究室お茶部屋． 米麦 0.5 合． タマネギ・ニラ・ショウガ・ホタテの炒めもの．
  • 晩 (2355): 米麦 0.7 合． ネギ卵炒飯． 乾燥野菜海藻スープ．