ぎょーむ日誌 2005-12-31

2005 年 12 月 31 日 (土)

• 0800 起床． 朝飯． コーヒー． 洗濯． 怠業 ……
• …… ぎみなんだけど， よゆーがないのでネットから北大にもぐりこんだり， 昨日の計算結果とかを調べ直したり．
• 通販で発注してたパンツ (ズボン) 三本とどく． 通販といってもゆにくろですが． いやはや， 数すくない手もちがことごとく穴あき状態なんで． スソとかに穴あいてるのは機能上はさほど不便はないんだけど， ヒトがみると 「何か特別な主義主張にもとづいて穴あきのままなのか?」 と誤解されがちなんで (あるいは大学院生たちから 「ちゃんとしてくださいよっ」 と譴責されるんで)， 今回の発注におよんだ次第．
• 素材は綿でも悪くはないけど， ひとつぐらいはウールにしとけばよかった …… あ， いま思いだしたんだけど， ウールパンツというとついつい 「野外でも使える良いものを」 とゆー厳しい規準で選ぼうとしてしまい， ゆにくろにはそれがなかったんだよね．
• 三本の中の一番の安物であるポリエステルものをはいて 1250 自宅発． 晴れてるけど雪ふり． 化繊でも意外と感触は悪くないかな． ちょい冷えるような気がするけど． 1305 研究室着．
• さて， 屋久島常緑樹の葉っぱ生き死に問題， おそらく最後に残された難問である (最後であってほしい， 最後ぢゃなきゃイヤだ) シュート伸長の休眠 & 二度伸びの効果の具体的なモデリングについて朝から検討している．
• なんでこんなものなんぞを考慮するのかといえば， ぢつは葉齢なる観測値は「芽鱗痕」を利用して憶測されるからだ． ある枝において芽鱗痕が毎年毎年ひとつづつ増えていくなら， 何も問題なく葉齢は特定できる． しかしながら， 年によって芽鱗痕がまったく増えなかったり， ふたつも増えたりすると「みかけの葉齢」と「ホントの葉齢」 は食いちがってくるのである．
• これまでの常緑樹の葉齢研究では， こういう「みかけ」と「ホント」のずれみたいな 都合の悪いことは「なかった」ことにされている (もちろん根性だしてマークつけ & 5 年以上もの追跡やるとゆー戦略もあるわけだが) …… そのくせ (別の文脈では) シュート伸長の休眠があーだこーだと言いたがるのは不思議なコトなんだけど．
• 伸長の休眠と二度伸びに関する hyperspecies - species - tree の階層ある階層ベイズモデルの推定はすでに終了している． おそらく一番首尾一貫してんのは， こういう休眠推定と葉死亡確率推定を同時にやることだろう． しかしながら， これはたいへんだと思うんで， 今回は見送り．
• で， 休眠まわりと葉死亡確率を分離する， という方針でいくなら …… うーん， 樹種ごとの事後分布の平均値を使う， という方向でいくのが無難そう． 個体差はどこにいったんだというハナシもあるが， これはおそらく 「休眠とかの個体差」と「葉っぱ死亡の個体差」 をこの観測データからは識別不可能だろうから， 葉っぱ死亡確率の個体差でもって代用すればよい， と．
• さて， これをどう MCMC 計算に組み込めばよいか …… と考えてみたんだけど， いやはや またしても Gibbs sampler が導出できないので， とりあえず Metropolis-Hastings 法つまりメトロポリスサンプラーで．
• いま， 「ある芽のある一年のふるまい」として (「確率」は事後分布の平均値である) ……
  1. 二度伸びする (根元方向のみかけの葉齢に与える影響 d = -1) : 確率 p
  2. ふつーに成長 (d = 0) : 確率 1 - p - q
  3. 休眠 (d = 1) : 確率 q
     
  がわかっていたとしよう． 休眠・二度伸び現象の発生は毎年独立であるとする．
• 今度は逆に 「ある芽鱗痕を選んだすると，それはどういう状態にあるのか」 を考えてみよう …… こういう確率分布で与えられる:
  1. 二度伸びでできた (d = -1) : 確率 p
     (二度伸びの場合は d = -1 しかありえない)
  2. ふつーの成長でできた (d = 0) : 確率 1 - p - q
  3. k 回休眠している (d = k > 0) : 確率 q^k * (1 - q)
     (k = 1 ..∞の和をとると q になる)
  これらを休眠・二度伸び現象に関する尤度計算につかえばよい …… これは事前分布 (prior distribution) ではない (この部分に関しては事前分布なし?)． たぶん．
• さて， よくわからないところが各芽鱗痕の状態遷移だ． 愚直にやるなら， 各 step において 全芽鱗痕ひとつひとつについて (d をランダムに ±1 増減させて) メトロポリスサンプリングすることになる …… これは計算時間から考えて不可能．
• 次におもいつくのが， いくつかの芽鱗痕をもっている「枝ごとに」めとろぽりすする， という策なんだが …… このときに， 新しく「提案」するパターンを生成するときに 全芽鱗痕で d を ±1 変動させたらほとんど 「却下」 されるような気がする (やってみないとわからないけど)．
• そこでこれを改善する妥協策として， ある枝を選び (枝といっても長さはさまざまなんだけど)， その中のどれかひとつの芽鱗痕での d を ±1 変動させたものを 「提案」してメトロポロス殿に審査していただく， という計算方法がありえるんではないかな， と考えている．
• これのまずそうなところは枝の長短 (正確には枝がもつ芽鱗痕数の大小) になンか影響されそうなあたり， かなあ …… まあ， こんなんで許されるんではなかろーか． ひとまず昼飯．
• お， Dell 機にやらせてた計算が終了している． 昨日の SWIG による高速化で 計算時間は 10 時間 40 分から 7 時間 30 分に「短縮」されてるわけだが …… あまり短縮感ないなぁ．
• 図の保管場所としてここが便利なんで， またしても陳列しておく (501-1500 番目のサンプルのみを使った事後分布図) …… 毎日毎日同じよーな図ばかりを， と思われるかもしれないが， これは昨日とまったく同じ条件 (ただし今回は SWIG 使ってるけど) なのに， 「だいたい同じだけど，びみょーに違う」 というところが気になっている． 観測データの交絡が原因なのか， MCMC サンプリングのやりかたに問題があるのか， それとも尤度まわりの定式化と数値計算にばぐが残存しているのか．

• 上の図で， Hyperspecies 階層の Narea なんてのは， おそらく MCMC 計算がうまくいってなくて (0.4 あたりに「ひっかかっている」ように見え， これは明らかに hyperspecies のパラメーター値として適切ではない)， その原因は …… 観測データとモデルの両方? つまり「よくわからん」ということです．
• よくわからんまま， Hyperspecies の層におけるパラメーターの遷移 (Narea を使ったモデル) の図を置いとく．

• 上の図をみると， step 200 あたりですごい「大ジャンプ」が生じてるよな …… これって MCMC 計算もうちょい工夫しろってことかしらん． パラレルテンパリング法 (拡張アンサンブル法)， とか． 計算がますます遅くなるな．
• パラレルテンパリングはしんどいので， simulated annealing 法もどき， とか? 最初は「温度」を高めにしておいて， 徐々に下げていく， とか． いや， これではあまり改善されんような気がするな．
• 「ごくまれに」 M-H ルールを無視させてみる? アヤしげではあるけど ……
• ついでに， 全体の fixed effects の事前分布の分散逆数の遷移 (Narea を使ったモデルの 501-1500 番目のサンプルのみ) の図も． ふーむ ……

• よくわからんけど， 事前分布の分散に関しては， hyperspecies → species → tree と「もぐる」につれて「分散の分散」がどんどん小さくなる． Narea もそうなんだけど， tree (個体) の層において他の fixed effects の事前分布より 「分散の分散」が大きい． これは何を意味するのだろう? Hyperspecies の層でのあてはまりが悪い， ためか?
• どうも挙動がよくわからんし， 休眠まわりのプログラム書きをやってるあいだは計算がススまんので …… ヒマそうにしてる A801 室の Dell 機に 「確率 0.01 で M-H 法の遷移ルールを無視する MCMC 計算」 なるものを試運転させてみる …… さーて， どうなることやら．
• で， 放置しとくつもりがついつい結果を見てしまう． ふらふらして面白い …… ってゆーか， 「確率 0.01」はほぼ「無秩序」なる状態におちいってしまうよーで (1 MCMC step の中で 200 回ちかくもめとろぽりすなんぞするからだ) …… だったら「確率 0.001」にしてみるか． お， 今度はちょっと見た目にはふつーの計算と変わらんな． ということで， このまま放置．
• しかしやはりついつい見てしまう． まだ 50 sample ぐらいしか取れてないけど (つまり定常分布にもなっていない)， なーんか当方に都合の良さそうな値になってるよな (明るさも窒素も効果は平均的にはゼロ)．
• しばらく見つづけてしまったんだけど …… 確率 0.001 でもかなり無法感が高い． それじゃ確率 1.0e-4 で …… いやこれは小さすぎるかなというこで確率 5.0e-4 ならどーだっ (きりがない)． これで RNA virus の突然変異率ぐらいになったか?
• プログラム改造のほうはススまん …… すでにアタマがばてきったか? ひとまず休憩． 外は大雪． 本日は院生たちも早めに撤退 (というかすでのたいはんが帰省)．
• 確率 5000 分の 1 で M-H 法を無視する実験的な MCMC 計算 …… 確率 5000 分の 1 でもかなり「ぼろぼろ」になるとわかった． つまり実用にはならない． しかしながら， いろいろなことがわかるので， とてもおもしろい． たとえば， この問題の尤度の「地形」が「ひどく複雑」である， といったあたりとか． 窒素なんてどうでもよさそー， というのはだいたいにおいて正しいような気がする．
• さて， ここに休眠とか入れたらどうなるか， だ． ようやくコードとむきあう． 時刻はすでに 1830．
• 上の「各枝の中のひとつの芽鱗痕をひとつランダムに選んでメトロポリス」 という方式でいくか． ついでに， (今までいいかげんにやってた) 「欠測末端」 も MCMC で選ばせる， と． 年内にどうにかならんかな． 残り 5 時間とすこし．
• とりあえず SWIG 整数配列の copy constructor みたいなものを作る． 単なる copy メソッドか．
• 上で検討していたような休眠・二度伸び・欠測の MCMC 部分 (ほとんど SWIG で計算) は簡単に片づいた．
• しかし意外な面倒に遭遇し， かなり最悪な状況に． (きちんと対応づけする必要のある) 休眠データシートと葉齢分布データシートの 学名とその略号が異なるのである …… と発見してしまった． 学名が異なるってのは参照もと文献を変えたからそうなったんだろうが …… 嗚呼， ひどい．
• たとえば， マテバシイが前者では Lied (Lithocarpus edulis) となっているのに， 後者では Pe (Pasania edulis) である． 1/3 ぐらいの樹種で変更されてるような気がする …… こんなのどうしろっていうの? とっくに帰省した牛原さんが仕掛けていったワナにちがいない． うう．
• ネットで和名を検索して複数の学名の対応をとっていく． 最終的に 22 樹種中 5 種が属名から変えられていて， その影響なんかでアルファベット順が変わったのは 9 樹種だとわかった． これが修正した対応表だ． 脱力．

        旧      新    (和名)
   1    Caja    Caj   tubaki
   2    Casa    Cc    sazannka
   3    Casi    Cs    sudazii
   4    Clja    Clj   sakaki
   5    Euja    Ej    hisakaki
   6    Ilan    Ia    sikimi
   7    Ilro    Ir    kuroganemochi
   8    Lied    Pe    matebashii
   9    Maja    Pj    hosobatabu
  10    Math    Pt    tabunoki
  11    Myse    Ms    taiminntachibana
  12    Neaci   Na    inugashi
  13    Neacu   La    baribarinoki
  14    Pona    Pn    nagi
  15    Qust    Qs    urazirogashi
  16    Rhta    Rt    sakuratutuzi
  17    Sybu    Sb    adeku
  18    Sygl    Sg    mimizubai
  19    Sylu    Sl    kuroki
  20    Sypr    Sp    kurobai
  21    Syta    St    onikuroki
  22    Viod    Vo    sanngozyu
  

• 時刻はすでに 2030． またひとやすみ． お， 外の雪はやんだ． 「確率 5000 分の 1 で逸脱」 実験のほうはあいかわらず「ひっそりと狂ってる」かんぢで興味ぶかい． 窒素と明るさの交絡問題， ありそうだな． 休眠とかいれても解決しない， ってことかしらん．
• プログラム改造， すすめてみる．
• 休眠・二度伸びのパラメーターに関しては …… 「和の事後分布の平均」 と 「事後分布の平均の和」 にはほとんど差がない．
• ざっとバグとり終了． ひとまずできた． 時刻は 2300． 100 MCMC step に要する時間， 35.3 秒から 39.0 秒に増大してしまった．
• 尤度が nan とか -Inf になったりしてすごくおかしい． SWIG 経由で呼び出してる C 関数内， 休眠まわりにバグが残存していた．
• ぽいんと: 枝の一番先端の芽鱗痕 (とゆーか翌年以降に伸長する茅) において「休眠」という状態はありえるけれど， 「二度伸びで形成された」 はありえない …… というところが正しく実装されてなかった， と．
• このあたり修正してみたら， ひとまずヘンな尤度は出なくなったので， (またヘンなものを見る前に) さっさと撤収． 2350 研究室発． 年内に大学脱出． 2405 帰宅． そして年明けに帰宅． ばてた． 晩飯．
• [今日の運動]
  • すごく運動不足な年末になってしまった ……
• [今日の食卓]
  • 朝 (0850): 米麦 0.6 合． ハクサイ・ネギ・エリンギ・卵の炒めもの．
  • 昼 (1430): 研究室お茶部屋． 米麦 0.6 合． ハクサイ・ネギ・エリンギ・卵の炒めもの．
  • 晩 (2420): 米麦 0.7 合． レトルトパウチドカレー． ハクサイ・ネギ・エリンギ・卵の炒めもの．