ぎょーむ日誌 2005-11-29

2005 年 11 月 29 日 (火)

• 0905 起床． 生活周期， さっさともどさなくては． 朝飯． コーヒー． 1010 自宅発． あられまじりの雨． 1025 研究室着．
• 毎週火曜日 1030 から研究室セミナーなんだけど， 予想どーりというか， 今回担当の某隠匿道師匠が (たぶん 「このあまりにも貴重なる観測データを 無知矇昧な連中の目から護りとおすには， もはやわが身ごと消え去るよりほかにすべなし」 といった理由で) 逐電してしまったので， 本日はナシ (ということは昨日わかってたのだが)． そういや， いかなる偶然か， 本日は雨の中から雷鳴が聞こえるな．
• 屋久島シュートモデリング， いちばん簡単なところからとりくむ． すなわち， シュート伸長の休眠現象の推定． さて， こいつをべいぢあんふうに定式化しようとすると …… 最初の一歩でさっそく足どめをくらう． 「樹種差」と「個体差」， この両者をどうやってうまく推定してみせるか， という問題だ．
• 問題じたいは簡単で， 無情報事前分布を用いてそぼくに推定させると なにかしら面倒なことになるかも， という気分になってきた． 次の状態は果たして区別できるのだろうか， ということ．
  • 樹種差はあるけど個体差はない
  • 樹種差はそこそこ個体差もそこそこ
  • 樹種差はないけど個体差はある
  さーて ……
• よくわからぬままふりーづぎみになっていると， 実吉さんからケムリ発芽実験種子データのさしかえ問題， というデータ解析こんさる． うーむ， これってどうだったっけ， と二人でふりーづしつつ作業をススめる．
• やったことはミヤマハンノキの死亡種子をとりのぞく， というもの． 呪われゑくせるファイルから CSV ファイル出力で データ救出して vimdiff してみると …… たしかにミヤマハンノキだけが違う， と確認． 新旧データを merge する Perl スクリプトを書き， あとは R で全処理 (ほぼ) 全自動化されてるんで， あっというまに解析・作図やりなおし終了．
• 昼飯． 女性大学院生密度高い． そして来週月曜日の某氏の予備審うちあげを口実に 「おいしいものを食べましょう」 欲求が表面化し嵐のようにお茶部屋でうずまいている． こそこそと逃げる．
• また屋久島シュート事前分布問題． よくわからない問題は， いつものごとく基本にたちかえって考えてみよう． たとえば， 樹種差・個体差の無情報事前分布をどちらも 独立同分布 (IID) の正規分布だとおいてみよう． そして (個体の状態) = (樹種差) + (個体差) だとすると， {樹種差, 個体の状態} は (定義から) 相関 sqrt(0.5) = 0.71 ぐらいの二次元 (二変量) 正規分布になる (ただし「個体の状態」のほうは variance 二倍)． なるほど．
• このことから， 樹種差と個体差をそれぞれ IID の正規乱数として生成させないですませるやりかたがある， とわかる． つまり (相関を自由に指定できる) 二次元正規分布を無情報事前分布としてみる， としてみたほうがすっきりしてるのかもしれない． この (まさに「個体差」の重要性を定める) 相関は超事前分布 (hyperprior) で指定してもよいし， 「無情報」などと称しつつもてきとうに都合の良い値 (> 0.71?) とかを与えてしまう， という方策もありえる．

[相関 0.71 の二次元正規分布の例]

横軸は標準正規乱数 (x とする)． 縦軸は x にもひとつ標準正規乱数くわえて 0.71 倍したもの． これは R の library(MASS) の mvrnorm() 使った場合の mvrnorm(n, c(0, 0), matrix(c(1, 0.71 * 1.41, 0.71 * 1.41, 2) ,2 ,2)) と同じ．

• まてよ …… ひとつの樹種はいくつかの個体から構成されるわけだよな． とすると， 個体数 n の場合は (n + 1) 次元正規分布となるわけか． まあ， 乱数生成に関しては二次元正規分布が n 個あるのとほぼ同じだが．
• そして， 無情報事前分布のごとく 「分散でかけりゃそれでいいさ」 といった投げやりとしか言いようのない状況では …… ハナシは結局のところ最初の IID 正規乱数の和， にたちかえってくることになる． なぜかというと， 「無情報」のもとでは樹種差・個体差どちらの寄与がでかいのか， それがわからないからだ． ここでの無情報相関は 0.71 なのである …… どんな事後分布がでるかは観測データまかせ， それがイヤなら相関を事前に指定せよ， ということだな．
• ということで， この場合の無情報なる {樹種差, 個体の状態} の事前分布とは， 正の相関がかなり高いように思える二次元正規分布なんで (つまりそういうふうに制約されてしまっている) …… まあいきなりわやくちゃな推定にはならいかも， と気分的におちつく (実際の事後分布はどうなるのかわかったもんではないが)．
• R こんさる． 小山さんと R で「サロベツ最近傍やちぼうず地図」を作る． 地図つくりはいつも楽しい．
• 屋久島シュート問題， 夜中までかかってシュート伸長休眠まわりの MCMC 計算による 推定部分の中核部を書いてしまう． 面倒な部分は， やはり現象の階層的な部分 (個体 in 樹種 in 樹木全体) であった． Metropolis-Hastings 法でパラメーターなんかの値を更新するときに， 尤度の評価が必要になるんだけど， このときに「自分より下」だけで尤度計算するようなしくみが必要となる．
• 今回も Perl で書いてる． R にはこみいったプログラミングが難しいという側面が， というか計算コストが重要になってくる場合の書きかた， というのを私がまだ理解してないんで．
• 重要な部分は書けた． しかし面倒が全部かたづいたわけではない． たとえば入出力まわりはまだ書いてないんだけど， これは退屈だけど面倒な部分だ． しかし本日は撤退． 2120 研究室発． まだ小雨ぱらつく． 2150 帰宅． 晩飯．
• 生活周期ずれぎみなんで， なんとか復旧したいところなんだけど …… 真夜中すぎに寝ようとして， 雑用放置を思い出してしまった． 12/9-12/11 の東京出張の宿泊・航空券の予約．
• で， ばててるのになぜかなかなか眠れん．
• [今日の運動]
  • 腹筋運動 30 × 3 回． 腕立ふせ 5 × 3 回．
• [今日の食卓]
  • 朝 (0920): 米麦 0.7 合． キャベツ・タマネギ・ブナシメジ・豆腐・エビすりみのシチュー．
  • 昼 (1350): 研究室お茶部屋． 米麦 0.5 合． キャベツ・タマネギ・ブナシメジ・豆腐・エビすりみのシチュー．
  • 晩 (2230): 米麦 0.7 合． ハクサイ・ネギ・ワカサギの味噌汁．