ぎょーむ日誌 2004-05-04

2004 年 05 月 04 日 (火)

• 1020 起床． いかん． 生活周期がまたずれてきている． 朝飯． コーヒー． 今日も雨か．
• 苫小牧直径成長モデリングの推定計算， R だけではどうしようもなく計算が遅いから， Rから他言語利用 による拡張関数かきの続き． とりあえず， 尤度関数評価のたびに全樹木のデータをやりとりするのは ばかげているので， あらかじめ C の関数に樹木データを「記憶」 してもらうことにする．
• これはどうすればよいか? …… おそらく static を使えばよろしかろう， ということでテストコードを書いてみる． 何の問題もなく動いた． ということで， 混合な確率論的モデルの尤度計算関数もこの計算方式に切り替えてみる．
• 計算速度は pure R に比べると， 10 倍ぐらいの向上かな． まだ遅いんだけど …… というのも， 尤度関数の評価以外はいぜんとして R でやってるからかも．
• 昨日は最大化された対数尤度がみょーな値になる， というバグが修繕できなかった． しかし本日みなおしてみると， C の「配列の添字」問題という基本的なミスだったことに気づく． いやはや． で， そこを修整してみると pure R で計算させていたのと， だいたい同じ答えを返してくるようになった． が， まだ完ペキではない． ちょっとだけバグが残存してるかんぢだ． 何か値のうけわたしに失敗してるよーな ……
• 1400 よーやくバグとれた． R コードの中で発生している計算ミスによって， C で書かれた関数の中でヘンな場所への数値参照が発生していた． つまり例によって配列問題． さすがは元祖呪われ言語 …… ああ， しかしながら， ともあれこれでよーやっとまともな計算結果が， それなりの速度で得られるようにはなった．
• これ以上の「R → C 置き換え」 は苦労がおおいわりには報われない気がする． たとえば optim() つまりここでは Nelder-Mead 法による パラメーター探索だけど， C 版 (Numerical Recipe 掲載のものとか) あるいは C++ 版 (久保自作) に取り換えたところで使いづらくなるわりには たいして速くならないだろう． ということで， 速度性能の追求はこのあたりで妥協せんといかんのでは， という気がする． そもそもモデルだってまだまだ変更されうるわけで．
• ThinkPad で計算試運転させといて， 1430 北大構内走発． いつのまにか晴． ここしばらく走ってないんで， ゆっくりとしたペイスを維持する． 1525 帰宅． 体重 73.6kg．
• 計算のほうは …… おお! 一時間に 180 セット以上も計算がススんでいるか． いやー， すばらしい． ThinkPad X31 (PentiumM 1.6GHz) でこうなら研究室の Dell 機 (Pentium4 2.4GHZ) ならかなりの範囲を探索できそうだな．
• おそい昼飯． 推定計算の経過をながめつつ． まだ狭い範囲しか調べてなんだけど， いったん止める． その時点での最良モデル．

  > best.model
  $aic
     1275.83
  
  $par
    scale.log      const    dbh.log dbh.log.sq    ps.prev        ct6 
     -1.01933   -9.40252    5.11672   -0.74374   -0.20741    0.29254 
  
  temp.selector: best=18.5/width=8.0
  rain: factor=0.20/time=0.30
  ppfd: factor=0.01
  temp: factor=0.05
  

  …… まだまだ改善の余地がいろいろあるような．
• それはともかく， とりあえず 1700 自宅発． 晴． 1710 研究室着． Dell 機で動かしてた計算はまだ 200 ちょっとしかいってないわけで …… まあ， 一日ついやして「C による拡張わざ」 を修得する価値はあった， と言いますか．
• 推定計算プログラムの細かいところを直す． 計算する場所を Dell 機に移す．
• で， ぢりぢり計算させてるわけだが …… もし今つかってる「前年の光合成関数」 が正しいとすると (この仮定がアヤしいんだけど)
  • もし気温依存性の幅が狭いとすると 影響はプラス
  • もし気温依存性の幅が広いとすると 影響はマイナス
  なる奇妙な性質がありそう， ということがわかってきた． いやー， 後者が「正しい」とすると， 繁殖とかの影響があるにちがいない． どうなんだろ．
• とはいえ， 「気温・降水量・明るさを見て ある年の生産速度を憶測する」 というモデルはこれから大変更を要求されるよーな気がする …… まあ， R における速度問題について一定の解決がついてる パラメーター推定・モデル選択システムがもはやできてしまったので， モデルの改善とかコワくはないわけだが．
• 昨日ご教示いただいた計算問題について調べてみる． 丹後俊郎 「統計モデル入門」 …… えー， この本は「入門」などとゆー題がついてるけど， ふつーの統計ユーザーが見てもほとんど理解不能な内容である． 私もまだとても全部はわからん． ときどき見直すと以前は意味不明だったところが 「こういうことだったのか」 と理解できるようになってる， といったじつに不思議な著作であったりする．
• この助言は フロリダ大学 の奥山さんからいただきました． いつもありがたいことです．
• で， その p.211-212 から示されてる Poisson * Gamma たたみこみ計算と申しますか， この本のながれでは emprical Bayes 推定用のモデルの構築， というような …… それはともかく Poisson 分布の混合で負の二項分布が導出されてる． つまり「こんなの積分できん」と騒いでたわけだが， Poisson パラメーターが積のカタチになっていてもやはり 負の二項分布になりますよ， というやや脱力な結着のよーで． 式変形の肝要な部分がトバされてるのは， いかにもこの著者らしい． 「こんなのあたりまえだろ」 というような． たぶん Gamma 関数もしくは Gamma 分布まわりの積分の例題にあたれば解明できそうだ．
• ついでにぱらぱらとその前後をながめてみる …… ふーむ， Gamma 分布は Poisson 分布の共役事前分布として良い， というのも事後分布も Gamma 分布になるから， とか書いてあるぞ． ここには書いてないけどネット上で調べてみると 自然共役事前分布 (natural conjugate prior distribution) とかいうそうで． で， 「ランダム変量の条件つき期待値」 などと私が呼んでる推定値は， この事後分布と関連づけて説明されてる． なるほど． たぶん今後はこういう Bayes な説明のやりかたが普及していくんでしょうなあ．
• 気象値パラメーターまわりを少しなおしてから Dell 機に再計算を命じておく． たぶん数時間ぐらいで計算終了するだろう． 2130 研究室発． 2140 帰宅． 晩飯．
• うげ． 肩のストレッチングしてたら， イヤな音をたてて肩が亜脱臼してしまった． ひさしぶりだ． なんというか， あたかも知恵の輪がはずれるように 腕の骨がずれてしまうんだよねえ． 「いてて」 と言いつつはめなおして冷蔵庫の製氷室 (当家の冷蔵庫には冷凍庫がない!) にあった冷却材で肩を冷やす．
• [今日の運動]
  • 北大構内走 1430-1520． ストレッチング．
  • 腹筋運動 30 × 3 回． 腕立ふせ 3 × 3 回．
• [今日の食卓]
  • 朝 (1000): 米麦 0.8 合． チンゲンサイ・タマネギ・エノキダケ・ホタテの炒めもの． ブナシメジ・豆腐・煮干の味噌汁．
  • 昼 (1600): 米麦 0.8 合． チンゲンサイ・タマネギ・エノキダケ・ホタテの炒めもの． ブナシメジ・豆腐・煮干の味噌汁．
  • 晩 (): 米麦 0.7 合． チンゲンサイ・タマネギ・エノキダケ・ホタテの炒めもの． ブナシメジ・豆腐・煮干の味噌汁．