ぎょーむ日誌 2006-01-09

2006 年 01 月 09 日 (月)

• 0900 起床． 朝飯． コーヒー． 1025 自宅発． 晴． 1040 研究室着．

[雪の札幌]

現在こんなかんぢでつもっております． 北大地環研前．

• 屋久島計算， 50 MCMC step とばしでもダメか …… 逆分散 hyperparameter scale = 10， シュート伸長の休眠・二度伸びは昨日のものを使った (こっちは 10 とばしで十分なのに)．

葉面積あたり窒素量モデル

葉重量あたり窒素量モデル

• それじゃ， 今度は 100 MCMC step とばしで． 計算終了予定時刻は本日の真夜中． 生残確率一億分の一問題も少し手をくわえる．
• なんかシュート伸長の休眠・二度伸びという「隠れ状態」たちが このへんの計算を遅く難しくしている …… ような気がする． ややこしいことに， 葉っぱ生き死にプロセスもこの確率に関与してくるからだ． 休眠回数を増やしたければ， 死亡確率は低くなければならない …… うーむ， この尤度評価方法がまちがいなんだろうか． いや， 「観測されなかった休眠」 の存在を推定するには， こうするのが正しいはず．
• このべいぢあんな状況でモデル選択みたいなことが， やはり必要となるのか? ちょっとそれはあとまわしにしたいんだけどな ……
• 締め切り もせまってきたので生態学会新潟大会要旨の作文． 30 分ほどかかる． 小林さんと加藤さんにこれでよろしいでしょうかと送信．
• 新潟大会といえば …… 「個体差データ解析」自由集会． えーと， 今回は自由集会の時間が 2 時間． 私が 50 分間 (45 分ぐらいにせんといかんな) ほど入門的なことをハナして， あとは他の 4 人の皆さんがそれぞれ 15 分間ぐらいかな．
• そして集会後の「うちあげ」とかもまた考えねばなりませんな． たぶん大阪のときよりは選択の余地が少ないだろうけど．
• 航空券はまだ取れず． 私はおそらく 3/24 (金) 午前出発の 3/30 (木) 夕方にむこう発， というかんぢかな． 3/29-3/30 はくだんの「統計合宿 in 新潟」なんで．
• と作業してるうちにも屋久島常緑樹葉っぱ生き死に推定 MCMC 計算 (100 MCMC step とばし) 進行中． 面白いことに， このデータの葉重量窒素量が葉寿命に与える影響は 「長命化」と推定されることが多いんだけど， じつは「短命化」の方向にも「尤度の高原」があって， そちらに拘束されることもときどきある． まあ， 数万 MCMC step もしないうちにそこから「落っこちる」 (「高原が崩壊する」と考えるべきなのか?) 程度の 「狭さ」みたいなんだけど． これに対して， 「長命化」高原のほうは数百万ぐらいは落ち着いていられるよーなので， こちらのほうが尤 (もっと) もらしいのである．
• いづれにせよ， 窒素量なる測定値はかなりノイズ的で， あまり信用できない測定量ということなのかもしれない (ということで「モデル選択が必要かなあ」とボヤいてた次第)． もうひとつの可能性としては， 樹種によって窒素量が多いと長命になったり短命になったりする， という樹種差が連続的というより むしろ長命化・短命化の二グループに分断されているかも， というもの． これは標本数が多くないこの観測データから明示してみせるのは難しいだろう．
• しかしこの対数尤度が数百とかそれ以上も一気に悪化する 「ジャンプ」 は心臓に悪い …… と言ってるうちにもまた「ジャンプ」して 「窒素たくさんあれば短命化の高原」が崩壊した (あるいは高原から滑落した)． メトロポリス法って， たしかに原理的にはこういう「一気に劣悪化」もアリなんだろうけど ……
• 窒素はたしかにノイズ的だ． 定数項パラメーターとすごく相関がありそう． 負の相関 …… つまり両者を足すとだいたい一定になる， というような． そこまで極端ではないのかな?
• 昼飯．
• おっと， 屋久島計算が暴走ぎみだな． しまった， 生残確率百兆分の一問題には対処してたけど， その逆の死亡確率百兆分の一問題への対処は不充分だった． SWIG から呼び出す尤度計算の関数， かなり作り変える． つまり天文学的な確率がでたら 対数尤度 をすごく悪くしてやる， という方式で． これでどんなむちゃくちゃな確率をつきつけられても 問題なく MCMC 計算できるだろう． また計算やりなおし． まあ， 明朝までには終わるでしょう．
• 上のように改良したら， 値の変動がえらくぎょーぎよくなってしまたような気がする． しかもよく変動してくれてるよーな． 単なる気のせいかもしれん． そしていつまた， わけわからん状態に遷移するかもしれん …… 油断禁物．
• と書いたら (書くからいかんのか?)， またジャンプした． しかし対数尤度の悪化は 400 ぐらい． とうぜん， もとの方向に回復していくわけだが …… これには 10000 MCMC step ぐらいを費してるよね．
• 「じゃんぷ」考察つづく． 対数尤度でいっきに -1000 悪くなるとか， これはありえないことだ． で， 「みかけ上」そうなっているのは， ここでいう対数尤度とは事前分布のそれもふくんだ量になっているから． これに対して (先週じたばたと考えていたように) この MCMC 計算に関するメトロポリス法では事前分布は除外して考える …… ということで， 確率論的モデルのデータへのあてはまりが劇的に悪くなっているわけではない， はず．
• たとえば R で min(log(runif(1000000, 0, 1))) とすると -14 ぐらいの値になる (極値分布ですね)． この程度のはずなんだが ……
• メトロポリスのところで一様乱数の対数をとって， 対数尤度差 (事前分布ぬき) と比較してるのがまずいのかしらん? とりあえずここは修正して再計算 …… あまり変わらんような気がする．
• あるいは事前分布と提案関数が同じ (そしてこれは独立サンプラーになっている)， という MCMC 計算がとてつもなく効率悪いとか? …… うーむ， Dell 機では現状のまま計算させといて， 手もとの計算機でちょっと調べてみるか．
• 独立サンプラーやめるんなら対称サンプラーだよな …… この問題に関する効率よさそうな Gibbs sampler はいまだに思いつかんので． ありきたりながら， がうしあんカーネルにするか． とびハバはてきとーに …… というか， この私の自作 Perl 階層ベイズ推定プログラムにおいて， 切り替えるなら何もかも対称サンプラーにせんと面倒になる， とわかった． HyperpriorGamma に生成させる分散逆数は …… かけることの exp(平均ゼロの正規乱数) とでもするか． これでも対称サンプラーだ．
• そして対称サンプラーへの変更にともなって， Metropolis-Hastings 法の評価方法も変える， と． 独立サンプラーでは事前分布は除外していたが (提案関数が事前分布と同じであったため)， 対称サンプラーでは事前分布を勘案しなければならぬ．
• 夕方， 石橋君のミズキモデリングこんさる． そこにトレイドオフらしきものがあったか …… ということで， ミズキ主幹成長すなわち高さ成長の概要はわかってきたよーな気分．
• 屋久島計算， 対称サンプラー改造後の試験運転． またしても小バグ撃滅に苦闘する．
• だんだん動作がそれっぽくなってきた． なんかまた「流れ」がよくなったのか? ふーむ， ここ三日間ほどの
  • 正しくめとろぽりすする
  • 可変平均な事前分布の導入
  • 独立サンプラーやめて対称サンプラーにする
  といった改造によって， ますます 「さらさら MCMC 計算」 になりつつあるような …… さーらさらゆくよ♪
• 試験運転の計算待ち時間にお茶部屋雑談あれこれ． 大学院生たちは常に私を驚かせる …… しかし世の中のヒトたちは， まさか理系大学院生のあいだでもかくも 「理科ばなれ」 (少なくとも高校・大学で教えられる科学ぎらい) が蔓延してるとは思わないだろうな …… 理科教育って破綻してんのかしらん?
• 屋久島計算， 問題なく稼働してるよーだ …… 20 MCMC step とばしの計算を命じて 2110 研究室発． 買いもの． 2140 帰宅． 晩飯．
• 2330 対称サンプラーを使った計算終了した． 20 MCMC step (休眠・二度伸びは 10 MCMC step) とばしとは思えんほど「すっきり」した結果だ． ここ数日の苦闘は何だったのだ ……

シュート伸長の休眠・二度伸び

葉面積あたり窒素量モデル

葉重量あたり窒素量モデル

• 検算も兼ねて， 50 MCMC step とばし計算を命じておく． どうなることやら．
• [今日の運動]
  • またうんどう休養日
• [今日の食卓]
  • 朝 (0920): 米麦 0.6 合． ダイコン・ニンジン・ジャガイモ・ナメコ・油揚・サケの味噌汁．
  • 昼 (1330): 研究室お茶部屋． 米麦 0.5 合． ダイコン・ニンジン・ジャガイモ・ナメコ・油揚・サケの味噌汁．
  • 晩 (2240): 米麦 0.7 合． タマネギ・ニンジン・セロリ・ブナシメジ・タチのシチュー．