ぎょーむ日誌 2007-09-19

2007 年 09 月 19 日 (水)

• 0730 起床． 朝飯． コーヒー． 湿原植物の競争モデルについての再検討， いろいろと． いくつか簡単な定式化の方法は思いつくんだけど …… 0930 自宅発． 曇． 0945 研究室着．
• さっそく ECOAS 八甲田湿原群のバイオマス分布モデルの改善． たとえば， 「自分以上の potential をもつ species にだけ悪影響をうける」 という単純化した競争モデルをゐんばぐすの step() で単純に実装してみると …… やはり， というか計算時間が (冗談ではなく) むちゃくちゃ に長くなってしまいましたよ． 200 MCMC step に 1942 秒! …… といった次第で， 三次元配列などというバカげた実装はたちまち却下される．
• BUGS 言語の限界か …… これはプログラミング言語のようなみかけをしているけれど， 実際のところは 「loop を展開して数式を並べる」 だけの機能しかないからなぁ．
• また「potential のベキ」の方向で探索してみるか …… みこみウスいけど．
• 2000 MCMC step 1130 秒でオワった． 収束はまだまだだけど， もうちょい， このあたり調べてみるか， という気分にはなった． potential^power で power > 1 を強制してるんだけど (すると優占種が有利になる)， けっこう差がつく方向に流れはじめている．
• また計算まち時間． 今回は 50 分ちょいか?
• もしこれでうまくいかなかったら …… そうだなぁ， また別の策としては， ちょっとバカげたかんぢではあるけれど， 個体群ではなく群落レヴェルで「収量一定の法則」 みたいなのが (局所的に!) 成立していて， 群落バイオマスはその種組成には「あまり」依存しない， というアイデアはどうかしらん? そうすると， 先日のプランクトンモデルと同じく， その限られたバイオマスを植物たちがあたかもぶんどりあいをしている， といった描像で …… まあ， これならわれらがマイティーゐんばぐすでも問題なく計算できるよね．
• ゐんばぐすは便利だが， こいつの性能限界に生態学データ解析の 統計モデリングの発想が制約されるのは何ともはや， な状況だ …… ということで， MCMC 本などながめてみる． ふーむ， slice sampler か …… アルゴリズムわかったけど， これって既知の確率分布でないときは計算たいへんそう．
• HBC (Hierarchical Bayes Compiler) って Haskell + C で書かれているんだ． それはいいんだけど， これって複雑な階層ベイズモデルを定義できるのかしらん． 演算子 ``<-'' がないような気がする．
• 3300 秒ほどで上記の 6000 MCMCM step 計算終了したんだが …… うーむ， 昨日よりはマシになったとはいえ， 種間競争は依然としてキビしい．

[まだキビしすぎる競争??]

またヌマガヤ (spc23) の例． たとえば他種がいなければヌマガヤだけで 800 g 超になるはずだったプロットで， (他種との競争の結果) 現実には 382 g でオワった (白丸)． ちなみにこのコドラートの合計バイオマスは 694 g (黒丸)．

• で， A 棟 8F の湿原植物専門家たる小山さん・江川さん (西村さんは不在だった) におそるおそるお尋ねしてみたところ …… 久保質問が不明瞭で何とも回答に困るんだけど， いくらなんでもヌマガヤだけになったからといって 現存量 600 → 800 g はヘンでしょう …… とのことで， がっくりしつつすごすごと撤退してきた．
• ということで， モデルをバイオマス「ぶんどりあい」型に変更して計算やりなおし． その間に昼飯．
• うう， 挙動もよくみないでいきなり長時間計算してもうまくいかないわけで …… いろいろと統計モデルの改善のための試行錯誤つづく．
• どうも latent な 「いる・いない」変数のあつかいがマズくて， MCMC 計算を乱しているらしい， とわかってきた． そのあたり修正してまた 2000 MCMC step の試運転． 収束のみこみが見えてきた! 計算時間も 620 → 500 秒と短縮されたし． 6000 MCMC step 計算を命じる． 25 分ぐらいかな．
• 次の生態学会大会 (ESJ55) の案内くる． 3/14 (金) - 3/17 (月)， 最後の日が懇親会， と今までとはちょっと異なる． しめきりあれこれ:
  • 企画集会，自由集会申し込み 10/26 (金) 17 時
  • 一般講演申し込み 11/16 (金) 17 時
  • 講演要旨登録 2008 年 1/8 (火) 17 時
  なんだか気が重い． ポスター発表は「アリの噛みつき確率」 かなあ …… 自由集会 はどうしよう． またベイズものになるんだろうけど． うーん， やっぱりベイズの基本的なところをやったほうがいいのかしらん?
• 今回は粕谷さんが大会の実行委員長なんで， すごくおいそがしそうな気がする．
• 6000 MCMC step 計算， だいたいいいかんじなんだけど， パラメーターのうちいくつかは 5000 あたりから定常になってる …… ということで 8000 MCMC step (burn-in 6000 MCMC step) に延伸してやりなおし． 二重にコドラート random effects が入ってるので (コドラートのバイオマス合計と各 species のぶんどり能力)， 苦闘してるのかしらん?
• 1800 秒弱で終了． 収束は問題ナシ …… しかし統計モデルがダメだ． 標高によって変化するはずの競争力に関する種差がまったくでない． そのあたりのめんどうをぜんぶ random effects に押しつけてるな ……
• よくよく BUGS code を見直してみると， これは二重 random effects うんぬん以前の問題として， 私の定式化がマズかった， と気づかされた． 例によって「いる・いない」問題で， やはりここは観察された Occurrence[i, j] と latent な occurrence[i, j] の確率を区別して dbern(...) せんといかんわけね …… じゃないと， species ごとの 「もうこの高度では生きていけません」 がうまく表現できないわけだ． 植物がどうやって標高を感知してるのかはまったくのナゾ ではあるけれど．
• そのあたり修正して， また再計算を命じる．
• ダメだ． 種差はでてるけど， でかたがおかしい． 低地にしかいないやつの標高依存性がプラスだもんなぁ． しかし random effects ということにして無理矢理補正してるよ．
• このあと 23 時ごろまで試行錯誤の泥沼にはまることになる．
• 最初に， 「そもそも『いる・いない』のモデリングなんて， もう不要じゃなかろうか?」 と考えてこれを削除してみた． これは正解だったよーで， パラメーターと deviance はすっきりすばやく収束するようになった．
• それはいいんだけど， 列強な species (各コドラートの上位 3-4 種) たちが chain ごとに「談合」のごとき状態に入ってしまって …… つまり， おたがい相談してるかのように 「なるほどあんたが競争力 xxx という値をとるなら私は yyy ぐらいの競争力があればいいんでしょう」 てなふうに xxx だの yyy だの値 (これはおもに random effects で決まっている) が chain ごとにばらばらに「固定」されてしまって， 混交しなくなってしまった．
• このあたり事前分布の変更で改善できないかしらんと工夫してみたり， いろいろやってみたんだが …… やればやるほど状況は悪化するばかり．
• しかも夜になってアタマがへたばってくると， 凡ミスも増えて時間の浪費速度がすさまじく加速される．
• とりあえず， 状況を完全に把握できるところまで撤退して本日は終了． 2300 研究室発． 2317 帰宅． 体重 68.8 kg． ECOAS うんどう不足ふとり． つくづく健康によろしくない． 晩飯．
• [今日の運動]
  • うんどう休養日 …… そろそろ生活をたてなおさなくては
• [今日の食卓]
  • 朝 (0830): 米麦 0.6 合． ニラ卵炒飯． ダイコンの中華ふうサラダ．
  • 昼 (1320): 研究室お茶部屋． 米麦 0.6 合． インスタント味噌汁． リンゴ．
  • 晩 (2340): 明太子スパゲッティー． キムチ & 豆腐．