KuboLog 2008-06-25

苦情・お叱りは，たいへんお手数かけて恐縮ですが，久保 (kubo@ees.hokudai.ac.jp) までお知らせください．

2008 年 06 月 25 日 (水)

0730 起床．朝飯．コーヒー．帯広開花モデリングについて検討してみる． 0940 自宅発．晴． 0955 研究室着．

あらら …… ThinkPad X31 の LCD が割れてしまった ……

[割れた液晶パネル]

今朝，デイパックにいれるときにちょっと落してしまったんだけど，まさか LCD が割れていたとは．しょうがないので自力で修繕すべくネット上で中古品をさがす．ありがたいことに，いまの時代なら送料こみ 5000 円ぐらいで X31 用中古品 (LCD パネルだけ) がすぐにみつかる．

LCD 中古品を発注．まあ研究室では ssh 経由で使ってる (たとえば，こうやってぎょーむ日誌かいてる) ので LCD 関係ないし，自宅でもこの程度の画面がみえてればメイルよんだりするぐらいはできるだろう．

気をとりなおして R-2.7.1 をインストール．これは昨晩かえりがけに Dell 機の中で RPM package 化させたやつだ．

帯広樹木開花データちぇっくのつづき．サイズ (DBH) の変化が異常 (と言いますか，こういうのはデータ記録・入力まちがい) の個体を自動的に除去するようにしてみる …… たいへん乱暴なやりかただけど，これをいちいち確認しながら手作業するのはたいへんなので． 13085 サンプル中 147 サンプルを除去．つまり全体の中の 1% のかわりものだ．

ねんのため書いとくと，このように，つまり「ダメそうなデータをてきとーに除去する」という方法はこんにちの統計学では推奨されてない．ではどうするかといえば，「人間が記録まちがいとかする」事象 (これは単純な正規分布なんかでは表現できない) まで統計モデリングに含める，というのがひとつの理想だと考えられてる．とはいえ，なかなかこう理想どーりにやるのはいろいろたいへんなので，とういわけで ……

データの性質をみるには統計モデリングしてみるほかあるまい …… ということで，「とりかかりの統計モデリング」として個体差とか場所差はいれてないけどいろいろな樹種差は考慮している開花の階層ベイズモデルを構築してみる．

だいたいできたので，計算させつつ昼飯．

計算時間 30 分ぐらいかかった．データ数が多いせいか， random effects が入ってないせいか …… ともあれ，収束は 2000 MCMC step ていどでだいたい OK? 結果はなかなかきょうみぶかい．

よくわからんけど個体差とかいれてみるか． (欠側部分のない) 全樹木個体数は

	> length(levels(factor(d$id)))
	[1] 2213

で， 2001-2004 年のあいだに何回観測されたか，は

	> table(table(d$id))

	  1   2   3   4
	779 793 514 127

つまり 1434 個体が開花状況を 2 回以上観測されてるわけで， random effects 的な「個体差」がそれなりに重要になってくるはず．

ついでながら，観測した調査地の数は ……

	> length(levels(factor(d$site)))
	[1] 20

うひゃー

「年差」はわけわからなくなるから，とりあえずやめとこうかな? random effects くみこみ版を作って計算やりなおし．これまた時間がかかりそうだ．

しかしながら 850 秒ほどで終了．計算時間が短くなったのは random effects を巧妙にいれたため，ということで．しかしパラメーター数が増えたせいか，収束がちょっと悪くなった．もっと長く sampling せんといかんな．

うーむ 3000 でもよくないようで …… いろいろ試行錯誤してみる．

計算まち時間はメイルをかいたり， R 自由集会記事を修正したり．

8000 MCMC step で 3100 秒ぐらい．うーむ，たぶん，まだモデルがよくないんだろなぁ．なかなかおもしろい結果ではあるけれど． BUGS code と事後分布表．

この統計モデリングは手ぬきなところがいろいろあるんだけど …… たとえば時系列データなんだけど， latent variable を使ったきちんとした時系列モデリングになってない，といったあたりだろうか．しかしこれをまじめにやろうとすると，開花だけでなく樹木のサイズや「明るさ」なんかも latent variable にする必要があるんで，なんというかたいへんなんだよねえ．データ整理整頓だの入力まちがいさいがしだのも今まで以上に苦闘することになりそうだし．

しかし時系列データモデリング的なところはぜんぜんうまくモデリングできてるような気がしない．昨日つくった図から，おそらく何種かは「去年開花したから今年は繁殖しない」 (しかも個体間で同期) というのがいたりするように思うんだが …… うーむ，やはり fixed effects と random effects がうまくあつかえてないような． (後記: 翌朝になってこの BUGS code で記述されてるようなモデルではダメで，「個体差」を状態変数としたモデリングが必要，とよーやく理解できた)

紺野さんからメイルいただき，開花ではなく結実のモデリングするなら雌雄異株な樹種の場合は注意すべきでしょう，とのこと．結実まではとてもテがまわらない予感．そして個体の雌雄データのあつかいは難しい．雌雄同株のやつもいるわけだし ……

とりあえず，また帯広にデータ整理と計算報告のメイルをかいて，本日のデータ解析したうけぎょーむは終了． 1845 研究室発． 1900 帰宅． 1930 自宅発北大構内走． 2010 帰宅．晩飯の準備．晩飯．

LCD 壊れた ThinkPad X31 もちかえってみたものの，メイルよみぐらいにしか使えん ……

[今日の運動]

北大構内走 40 分間．ストレッチ．

[今日の食卓]

朝 (0800): 米麦 0.6 合．ニラ卵炒飯．
昼 (1320): 研究室お茶部屋．米麦 0.8 合．レタス・キャベツ・ニンジン・キュウリ・ミズナ・ネギのサラダ．コマツナ．
晩 (2100): 米麦 0.7 合．ネギ・ホッケの電磁波酒蒸し．チンゲンサイ・ブナシメジの炒めもの．

ぎょーむ日誌 2008-06-25

2008 年 06 月 25 日 (水)