ぎょーむ日誌 2007-12-01

2007 年 12 月 01 日 (土)

• 0920 起床． 朝飯． コーヒー． 洗濯． 怠業． 1140 自宅発． 曇． 1200 研究室着．
• さて， 昨晩かえりがけに「投げて」おいた MCMC 計算だが …… 2000 MCMC step (トリプルチェイン) に 5.7 時間費やしていたか． 収束はいうまでもなく悪い．
• log(体積) → 定数 + log(長さ) + 2 log(太さ)， なる分割は意味不明ぎみだけど， つじつまはあっている …… などと検討しつつ昼飯．
• とりあえず (現時点では諸悪の根源と断定せざるをえない) 軸の太さ・長さの推定をヤメて， つまり問題をすごく簡単化して挙動を調べる試験運転 …… 初期値のあたえかたとかに改善の余地アリとわかった．
• 次に無情報 (に近い) 事前分布で軸の太さ・長さの事後分布も いっしょに推定させてみると …… 2000 MCMC step 計算時間 750 秒ほどで終了． 樹木個体数 500． ここまでは， まあ， こんなモンか， と．
• 次に軸の太さ・長さに関して軸重量の制約をうけるようにしてみると (つまり材密度の推定なるアヤしげな計算が介在してくる) …… 同じ条件で計算時間はイッキに 1800 秒弱まで増大． これは階層ベイズモデルを表現する DAG (directed acyclic graph) が複雑化というか deterministic node の経路が長くなったため， プラスその他の理由で事後分布にかかる制約が増えたために Gibbs sampling に時間がかかるようになった， といったところなんだろう．
• ということで， DAG 上で deterministic node の経路をぶった切ることで MCMC 計算は高速化されるであろう， と期待される …… そこで「軸の体積」なるアヤしげな潜在変数を log.v[i] ~ dnorm(mean.log.v[i], tau[1]) などと stochastic node 化してやると (そのココロは「軸の体積の個体差」といったところで) …… ほーら， 計算時間は 860 秒ちょいまで短縮された．
• だいたい挙動がわかってきたので 6000 MCMC step ほど計算させてみることに． 待ち時間は 40-45 分とゆーところか? メイル雑用とか．
• 2540 秒計算終了． うーむ， いくつかのパラメーターが収束してない．
• 葉の食害まわりが悪いのかしらん? …… とそのあたりをベータ分布化してみるけど， そうなると食害 rate 小さなところで slice sampler がコケる， と． 呪われぎみな dunif(...) にもどす． 状態変数 damge[i] に初期値を与えてないことに気づいたのでてきとーに与えてみる． 2000 MCMC step 試験運転からやりなおし．
• ほへー， qmail (+ その他 DJB software あれこれ) が public domain software に …… じつは私は qmail server 管理者だったりするわけだが (研究室の mailing list なんかは qmail + ezmlm で動いている)． どうなることやら．
• damage[i] の初期値を与えても改善は微々たるものだな …… ということで， 「同化部・非同化部間の重量分割の個体差」 をいれてみる． BUGS code はまた汚くなったわけだが ……

  mean.logit.pf[i] <- (
  	bb[3] + bs[Sp[i], 3]
  	+ (bb[4] + bs[Sp[i], 4]) * log.height[i]
  	+ (bb[5] + bs[Sp[i], 5]) * (log.w[i] - MeanLogW)
  )
  logit.pf[i] ~ dnorm(mean.logit.pf[i], tau[7])
  logit(pf[i]) <- logit.pf[i]
  log.wf[i] <- log.w[i] + log(pf[i]) + log(1 - damage[i])
  log.ws[i] <- log.w[i] + log(1 - pf[i])
  

  しかし， まあこれでまた一個 deterministic node の経路を切断できたわけで． MCMC 小人さんたちも働きやすくなったことだろう．
• MCMC 小人さんたちから「やりやすくなった!」 との報告． いままで 860 秒を費やしてた 2000 MCMC step 計算が 430 秒に．
• このあたり， 「えーい， どうせこの観測データでは個体差と測定誤差の区別はつかんからなぁ ……」 とこの両者をごちゃまぜにしてたモデリングに， 「この部分は個体差がきいてることにしよう」 と変更したおかげなのである． つまり適切なる random effects 的なパラメーターの導入は WinBUGS の Gibbs sampling を (速度面だけでなく) 改善する， と．
• とうぜんながら 不適切なる random effects 的なパラメーターの導入は， たとえば「同じことを二つのパラメーターであらわす」 状態になったりして， 収束はめちゃくちゃになる， と．
• このあたりのワザの術理は 9 月すなわち あの「悪夢の ECOAS 下うけ計算」の渦中で (シグルイふうに言えば) つかんだものなんだよね …… まあ， どんな練習も無駄にはならない， と申しますか (すみません， リアル 7 巻の戸川スリーポイントシュートに感銘うけたので)．
• ついでにもういっちょ， 非同化部の中の軸 + 葉柄の重量分割も個体差的な random effects をいれてやろう …… 計算は 430 秒から 360 秒に． そしてトリプルチェインの収束は格段に改善された! よしよし， これなら 4000 ぐらいに計算 step 数を延伸してやれば十分だろう． 12 分ぐらいで終了するはず．
• 予想どーり， 718 秒で終わった． 今日はこんなところか? ついでに log(体積) → 定数 + log(長さ) + 2 log(太さ)， 分割も「個体差」を許容してみることに．
• 計算結果はすぐに得られるだろうけど， 作図とかにはまりこんでしまいそうなので撤退． 1835 研究室発． 1855 帰宅． 晩飯の準備． 晩飯．
• 帰りがけにやらせておいた計算， さらに時間が短縮されて 635 秒 (4000 MCMC step) で終了， 収束もさらに改善された． いやー， 2000 MCMC step に 5.7 時間かかったときはどうしようかと思ったよ …… 本日は計算速度を 65 倍に高速化できたわけだ． MCMC 計算って面白くてオソろしいものですねえ ……
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (1000): 米麦 0.7 合． ダイコン・ニンジンの味噌汁．
  • 昼 (1320): 研究室お茶部屋． 納豆巻． リンゴ．
  • 晩 (2130): 米麦 0.8 合． アサリのすまし汁． ナス・ピーマン・タマネギ・エリンギ・鶏肉の味噌炒め．