ぎょーむ日誌 2006-01-04

2006 年 01 月 04 日 (水)

• 0910 起床． いかん， あいかわらず時差が． 朝飯． コーヒー．
• A 棟 8F でやらせている計算をみると ヘンなパラメーター値まわりでとらっぷされている (例によって例のごとく面積あたり窒素量のほう)． これは初期値で「すごく」むちゃくちゃな値 (例: 一年間に葉っぱが全部死ぬ確率は 99.99%) が設定されてしまうと， その状態から MCMC 計算では脱出できなくなることがあるため． いやはや．
• とりあえず， 計算を最初から全部やりなおすよう命じる …… とうぜんやりなおし計算も 「それっぽいところ」から逸脱した初期状態からスタートするわけだが， この程度のズレなら計算をススめていくうちに「補正」できそうだ．
• こういうのを MCMC 計算に自力復旧させようとするなら …… 拡張アンサンブル法 (たとえば parallel tempering) かしらん． このあたりについては 伊庭さん の「ベイズ統計と統計物理」 などみてください．
• そして自分の屋久島計算コード (in Perl + SWIG) をよく見てみると， これにパラレルテンパリング法を適用すんのはとてもラク， とわかった． 初期化のところで $hyperspecies を一個ではなく二個準備してやって， 片方では (メトロポリス用の) 「温度」を上げてやって， あとは「てきとー」に入れ換えていくだけ， だよね．
• しかしこっちの方向にはまりこんだりすると， またどんなワナがまちかまえているかもわからんので， まあ慎重に． 今月中旬は授業の準備とかで時間をとられるので， そのあいまにでも試運転させとけばいいよな．
• てなことをうだうだと考えていて， 1015 自宅発． 晴． 1030 研究室着．
• 昨年 8 月に Ecological Modelling にアカマツ論文投稿したんだけど， いまだに返信ないので年末に問い合わせてみた． そしたらこれまた返事がおそくて， しかも 「reviewer 二人とも忙しすぎてまだ査読できてない」 だそーで． Ecological Modelling， 掲載された論文とかみてもすごく時間かかってるんだよね． たぶん， 地球環境なんちゃらモデリングのたぐいの投稿とか増えて， 対処しきれていないんではないかな．
• 甲山さんがふらふらとやってきたので， 年末年始の修論進行状況 (屋久島 & ミズキ) を説明 …… なンか私が修論生みたいなんですけど．
• そしてお茶部屋に逃げるとまた亀山さんとツガザクラ雑談モードに． これは面白くて， ついつい時間を使ってしまう．
• 本日から開店している北大生協で昼飯を調達してきて昼飯． 院生たち徐々に研究室にあらわれる． 修論生から 「久保さん， ヒマそうですね， ヒマですよね」 とさぐりをいれられて， ふるえあがる．
• 昨日からやらせている屋久島計算 (この計算の発注者の牛原さんはまだ帰省中である)， やっぱり値のうつりかわりがぢりぢりしてる， というような． やや遅い． 特定の値にはまりこむとかなり動きにくくなる．
• そこで (他の修論生から新たなる下うけ仕事を押しつけられる前に) パラレルテンパリングな MCMC 計算のドライヴァーを書いてしまう． これは上述のように簡単で， 15 分ぐらいでできてしまった．
• と思って， 「ベイズ統計と統計物理」 をみなおしたら， 系列を交換する Metropolis-Hastings 法のところがヘンだった． そうか「温度」が違う chain を比較するときは， 温度で規格化みたいなことやらんといかんわけで …… と修正．
• で， 試運転してみたんだけど， 「温度」の設定がむずかしいなぁ． 温度二倍とかすると， ぜんぜんあてはまりが悪くて高温・定温の系列が交換してくれない． しかし温度を一倍に近づけるとパラメーターの値が高温系列でも あまり変わらなくて何のためのパラレルテンパリングかわからない …… しかも温度比を 1 に近づけてもなかなか交換が発生しない． 高温・低音系列二本だけではだめなんだろうか．
• などとうだうだ勉強してるうちに …… ぢつは Metropolis-Hastings 法を 根本的にまちがっていたかも疑惑 なるものが浮かびあがってきて， かなりアセる． 数時間いろいろ調べての結論だけ先に述べておくと， (ほとんど奇蹟のような偶然のおかげで) 「ほぼ正しい」 けどびみょーにまちがっている， というものだった …… なんだか収束が遅かったりするのもこのへんに関係あるのかもしれない．
• むろん計算はまた全部やりなおす必要ある． シュート伸長の休眠・二度伸びもふくめて．
• ところでアセってるあいだに， (びみょーにまちがってると判明した) いつもの計算結果もでてきたので， まあ図示しておこう． 葉面積あたり窒素量モデル / 葉重量あたり窒素量モデルの順． 階層プロットは 501:1500 のサンプルを使用．

• 上の結果， いまいちヘンだと思ってたんだけど， 原因はめとろぽりすまわりだけではなかった …… 真夜中すぎになってその原因わかった (ああ，修論おいこまれ生活)． 線形予測子の構成がヘンだった． 窒素と加齢の積になってて， ですね ……
• さて， Metropolis-Hastings 法まわりの問題だが …… いろいろ調べてみたんだけど， ここでは 丹後本 の説明を借用してみる．
• まず私のまちがい M-H 法から始める． 実際に屋久島計算で使っているのは多数のパラメーターを ひとつづつ変化させていく single-component M-H 法である． ここでは説明のため一パラメーター α だけからなる問題におきかえる． 尤度を パラメーター α の事前分布を としたときに， と変化する確率を
  
  としていた． これはまちがいである．
• 正しくは α のもとで α' がサンプルされる条件つき確率 なんかで下のように「規格化」 (というか条件つけというか) してやらねばならない．
  
• しかしながら， このときに α' が といった独立サンプラー (independence sampler) で生成され， さらにその関数型が事前分布と同じ， つまり であるなら， と変化する確率は
  
  と簡単化される． つまり， 私のもとの「まちがい式」では分母分子の 事前分布 などがよけーなのである． 二重に を勘定してることになるのかな．
• ところで， この屋久島モデリングの葉っぱ生き死に速度の推定の部分では， もうひとつ考慮しなければならぬ部分がある． それはシュート伸長の休眠・二度伸びにともなう 芽鱗痕での「みかけの葉齢」補正効果の部分だ．
• ぢつはこの部分もたまたま (とゆーか本能的にとゆーか) 問題の少ない方法を選んでいた． 各芽鱗痕での「みかけの葉齢」補正量はある MCMC step において ランダムに ±1 増減する「提案関数」になっている． これは対称サンプラーあるいは Metroplis sampler (random Metroplis sampler) と呼ばれてるらしい． 私は M-H 法つかえば何でも Metroplis sampler と呼ぶのかと思っていた．
• ともあれ補正量 1 → 2 と変化するのも 2 → 1 と変化するのも同じ確率 (0.5) なので， 尤度比だけ見てればよい， ということだ …… つまり従来どおりでよい． だいたいの部分においては．
• ただしこれが通用しない部分が一ヶ所だけある． それは各芽鱗痕の「(根元方向の) みかけの葉齢補正値」 として取りうる値は

  -1, 0, 1, 2, 3, ...

  となっているためである． 二度伸びをしたときには -1 補正が必要なんだけど， これより小さくなることはない． 現状では
  • 0 → -1 と変化する確率は 0.5
  • -1 → 0 と変化する確率は 1
  となっていて対称ではない．
• まあ， これを一番簡単に修正するには …… うーむ， 「-1 → 0 と変化する確率」を 0.5 にするかな． そして残り 0.5 の確率で -1 にとどまる， と． そしてこのときは事前分布をいれて Metroplis-Hastings してやらねばならない．
• 計算プログラムをいろいろ修正して， 計算全部やりなおしを命じる． まずはシュート伸長の休眠・二度伸びモデルのほうから．
• ついでながら …… 昨年 8 月ごろ にやってた， 母子里林冠推定の MCMC 計算 (これはベイズ推定ではない， つもり)， このときもあまりよく考えずに計算してたんだけど， これまたいかなる偶然なのか問題なさそう． というのも， ランダムに箱を選んでランダムに葉っぱを増減させる， とゆー「提案関数」が対称サンプラーのカタチになっているからだ．
• A 棟 8F メイルサーヴァーの qmail が壊れてしまったので再構築． いやはや． 原因はおそらく先月の OS 更新 (Vine Linux 3.2 に) したこと． /var/log/maillog なんかにも

  Jan  4 21:06:38 hosho qmail: 1136376398.511651 delivery 148: success:
  Incorrectly_built_binary_which_accesses_errno_or_h_errno_directly.
  _Needs_to_be_fixed./Incorrectly_built_binary_which_accesses_errno_or_h_errno_directly.
  _Needs_to_be_fixed./Incorrectly_built_binary_which_accesses_errno_or_h_errno_directly.
  _Needs_to_be_fixed./did_0+0+2/
  

  といった呪詛でうめつくされているんで．
• こちら や こちら を参考にしつつ， qmail を 5 年ぶりぐらいに再構築． とどこおっていたメイル配信が再開された．
• しかし上記エラーはとまらん． 調べてみると， てきとーぺえじ にたいへん役立つメモあり． error.h の extern int errno; を #include <errno.h> と変更してまた再構築 (といっても 10 秒ぐらいで終わるんだが)． どうやら， これでひとまず問題解決したもよう．
• 計算やりなおし， ぢりぢりと進捗中． 予想どーり， 結果の定性的な部分はこの修正であまり影響うけてなさそう …… というか， このシュート伸長データの場合， そもそもそれほど複雑な結果がでるはずもないのだが．
• あれこれチェックしてから撤退． 2230 研究室発． 2245 帰宅． ばてばて． 晩飯の準備．
• 100 分ほど費してシュート伸長の休眠・二度伸び確率の推定計算終了． A 棟 8F から結果もってくる． 今回は 50 MCMC step とばし 1500 sample． 見た目はほとんど変わっていない． めでたしめでたし．

• 晩飯． すでに真夜中すぎ． 時差が．
• [今日の運動]
  • 今日もまたうんどう休養日か …… 昼間にどうこうする， といった方向で考えるか?
• [今日の食卓]
  • 朝 (0920): 米麦 0.5 合． カレー．
  • 昼 (1340): 研究室お茶部屋． 北大生協の巻きずし．
  • 晩 (2410): マカロニ． ハクサイ・ニンジン・ブナシメジ・豆腐のシチュー．