ぎょーむ日誌 2007-04-04

2007 年 04 月 04 日 (水)

• 0810 起床． ねむい …… コーヒー． 朝飯． 1040 自宅発． 曇． 1055 研究室着．
• 後藤さんトドマツ問題のつづき． 昨日までの元データ図表化でだいたいの傾向がわかったので， 今日はそいつを階層ベイズモデル化， BUGS 言語でそれを表現してみる …… あらら， けっこう簡単な問題なのかも? しかし眠くてへろへろぎみ．
• で， 例によって R2WinBUGS で計算させてみようとしたんだけど …… えーい， Vine Linux 3.2 → 4.1 upgrade にともなって Wine の再コンパイルが必要になった． これは 30 分以上かかる．
• …… で， そのあたりなんとかなって， トドマツ母子距離データの 階層ベイズモデルの MCMC 計算 なんかもわりと簡単に収束したんだが …… なに? この平均種子散布距離 200 m とかいうのは?
• ぢつはこういう結果が得られるかも， という予想があった． 昨日までの元データ図表化によって， わかっていたんだけど，
  • 母子間距離は母樹によってさまざま
  • 近距離にたくさん子供がいる場合もある
  • 子供は少ないんだけど， それがけっこう遠距離にいたりする
  といったパターンがある．
• ところが単純な統計モデルではこのあたりの区別がつかないので， あたかも一個体が近距離・遠距離散布の両方をやっているかのように みなされてしまう …… もちろんそれが正しい場合もある． しかしながら， よく考えずにわりと複雑な (というか，ときとしてかなり奇怪になりがちな) 関数形をもつ「散布カーネル」工夫の曲芸に走りがちなんだよね， このへんのヒトたちは ……
• 遠近両用カーネルモデルの場合， 「子だくさん」母樹こそが近傍にも遠方にも自分の子供を配置しているはず， なんだけどこの富良野トドマツデータでは必ずしもそうはなっていない …… ということで， 私はあえて単純な散布カーネルを仮定し， これが個体によって異なる (この「個体差」の正体はいろいろと議論可能) というモデルを作ってみたつもりなんだが ……
• しかしながら， 推定された散布カーネルは「べたっ」とつぶれたもので， しかも「個体差」はない …… くそう， 計算をさぼっていやがるのか? いやいや当方のモデリングに何か問題あるにちがいない．
• で， 沈思黙考して到達した結論は， 私の凡ミスというべきしくじりの結果として， こうなってしまったのだとわかった． ある母樹からみて「子供がいる倒木」だけに着目してたんだけど， そうではなくて「(たとえば 200 m ぐらい遠方にある) 子供がいない倒木」 もふくめてパラメーター推定しないといけない， ということだ．
• えーと， だとすると 全母樹 (子なしも含む) × 全倒木 (トドマツ稚樹なしも含む) × 全 (稚樹の) 年齢 のデータ行列を作らんといかんわけ? こりゃー， 100 万要素数ぐらいの莫大なモノになってしまうな．
• とりあえずデータ変換 …… えーい， 距離行列の計算がおかしいと思ったら …… 元データゑくせる上の入力まちがいぢゃん． たった一文字のちがいでここまでぼろぼろになるとは． はい， ゑくせる上でのデータ操作はやめましょう．
• ためしにその巨大データをゐんばぐすに食わせて計算させてみたんだが …… やれやれ， 完全に沈黙というか， すごくすごく絶望的に時間がかかる．
• ここから「計算可能」な世界にむかっての凄惨な撤退戦闘が始まる． まずは「子なし母樹」は除くことにしよう …… 「子なし」といってもたまたま plot 内に「母子鑑定ずみ」 子供がみつからなかった， というだけのことなんだけど． すると 母樹 × 全倒木 (トドマツ稚樹なしも含む) の組合せは 8848 まで減少する． しかもほとんどは「倒木上に子供なし」状態だ (いっぽうで， ある倒木上には 33 個体の子供が乗っていたりする)．

  > summary(as.factor(n.obs))
     0    1    2    3    4    5    6    7   11   14   15   20   24   25   27   33 
  8678  125   14   12    6    1    2    2    1    1    1    1    1    1    1    1 
  

• これでも実生・稚樹の齢分布までチェックするモデルは計算が遅すぎて 使いものにならん …… で， 今のモデルの仮定について再検討してみると， 各母樹の種子生産量には年変動はないといった仮定をおいているので， これをうまく悪用できそうだと気づいた． つまり (ここでは幾何分布で表現されていた) 齢分布を考慮しなくても， 1 パラメーターだけで「倒木のよしあし」が表現できる， ということだ．
• むろん「年変動がない」というのがヘンなのは当方ですら承知しているんだけど， この観測データはそういう年変動が見えるような情報はほとんど含まれていない． なぜなら
  • ほとんどの母樹は子供数がすくない
  • 子だくさん母樹は少数いるけれど， 子供の各コホートは個体数が少ない
  といった理由で「年変動」が見えないためだ．
• ということで， そのあたりも「くりこみ」したモデルに変更して …… ぎりぎりの後退に後退を重ねた末に敵さんを 当方の得意ワザでもって始末できるところまで引きこんで撃退， というところか …… 2045 よーやくにして， マトモな計算結果が得られた． というかこれってなかなかスゴイよ!

• 計算 step 数が 500 なんでまだ収束よろしくないんだけど (翌朝もっと長 step 計算やって上のとおりでだいたいまちがいない， とわかった) …… log.sd なる散布カーネルの幅， 「平均的」には exp(3.4) = 30 m ぐらいで， これが母樹によって (log.sdMt) そうとうに違っていそう， という推定結果だ． つまり複雑怪奇な散布カーネルなんぞいらないんじゃないの， というもの．
• まあ， 遠近両用 (ここでは種子散布の) 散布カーネル不要， というのはちょっと極端に思うかもしれないけれど …… じつはここでは各母樹の散布カーネルパラメーターの事後分布の 「もと」になる事前分布 (これも階層ベイズモデルによって観測データから特定されたもの) が重要なんだよね． 散布カーネル「母樹間共通部分」 をこの事前分布で積分するとおもしろい「散布カーネル」 が得られるだろうし ……
• つまり， ひとつの母樹の「散布モード」が遠近の両モードがあって， 年によってそれが変わる， みたいな解釈も可能かもよ， と．
• 計算うまくいってよしよし & 計算ばてでへとへと状態で 2110 研究室発． 疲れたので外で晩飯． 2210 帰宅．
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0930): 米麦 0.5 合． ネギ・ブナシメジ・豆腐の味噌汁． 一昨日ののこりもの， つまり …… ネギ・タマネギ・ショウガ・ニンニク・ブナシメジ・ 豆腐のオイスターソース炒め． コマツナごまあえ．
  • 昼 (1330): 研究室お茶部屋． 米麦 0.5 合． ワカメスープ．
  • 晩 (2130): 地環研最近傍のおいしい和定食屋 ごはんやデンスケ の「今晩のおすすめ」， 750 円．