ぎょーむ日誌 2009-05-26

2009 年 05 月 26 日 (火)

• 0600 起床． 朝飯． コーヒー． 0725 自宅発． 晴． 0740 研究室着．
• セミナーの予習．
• 1030 より 研究室セミナー， 本日は川久保さんで火山噴出物・種子散布・窒素固定植物・地温の 4 論文の紹介． このうち 3 つは Journal of Vegetation Science なのだが …… 空間構造・時系列・多種系というもっとも複雑なる生態学データを とっておきながら， 可能なかぎりヘンてこにデータをあつかう論文を掲載するぢゃーなる …… とゆー気がしてきた． 火山とか極端環境でうまい研究をするのは， 案外むずかしいのではないかな．
• 昼飯．
• 一次元空間モデリング， 試行錯誤のつづき． 脱・Gaussian random field (GRF) のため有限状態の Markovian Random Field (MRF) を検討してみる (注: 時系列データ解析方面のモデルと混同しそうでややこしいのだが …… MRF の定義は「各地点の状態数が有限」ではなく， ある地点に影響を与える「近傍集団」が有限の範囲におさまる， ということ． したがって WinBUGS の car.normal() な GRF もまた MRF の一種である)．
• 出発点としては， 各地点において {0, 1} の二値をとる MRF を考える． これを観測データにあわせるためには， 「観測をまちがえる確率」 みたいなものを導入しなければならない …… そしてこれは観測データから推定することは不可能， とわかってきたので P.err = 0.001 などと決めてやらねばならない ……
• で， とうぜんながらというか， 観測データと対応させようとすると， 各地点の値はほとんど動かなくなり (遷移は上述の P.err にほぼ支配されている)， あんまし確率場 (random field) ってかんぢではなくなる ……
• そこで， べるぬーい (dbern(p[i])) からばいのみある (dbbin(p[i], N.size)) に拡張してみる ……
• こういう新しい試みはなかなかうまくいかない …… とりあえず， 途中で一時間ほど藤沼君と現状について検討してみる． 他人と検討してみると， 少しアタマが整理される．
• また別のデータ解析こんさる 40 分間ほど．
• 来週の土日の北大祭雑用， 静かなる混乱というべき状況にあるようだ． 担当者と事務方面 (全体ならびに各専攻) がばらばらに動いている． 全体の事務から必要なら出勤とどけをだせ， というメイルがきたんだけど， とくに仕事は assign されていないので， そういう紙きれは出さないことにする． まあ， どうせ土日は大学にいるだろうけど …… 土曜日は勉強会もあるし．
• 二項分布な MRF， 試行錯誤しながらその性質を理解しようとしている …… だんだんわかってきた．
• 1845 研究室発北大構内走． てけてけ走っているとアタマの中が整理されてくる …… 1930 研究室もどる．
• ということで， これでうまくいくだろう， と思いきってパラメーター数を減らしてみたら， うまくいった! 翌日以降， またいろいろと変更するかもしれないので， とりあえずここに記録しとこう (BUGS code， 事後分布表)．

• さて， この程度のことが GRF ではうまくできないのか， と再確認してみたくなった． まず Distance 依存性がない最単純モデルによる推定 (BUGS code)． この程度の計算は， さすがに収束する．

• ちょっと横道にそれるが， 私の場合， こういう単純なモデルによって GRF のような (じつはけっこう難しい) モノの性質の理解が深まる． たとえば上の binomial なやつとこの GRF なモデルの相違は， intcpt なる「切片」パラメーターの有無にある． GRF なほうでは intcpt ナシではうまく収束せず， binomial MRF なやつだとこういう「切片」パラメーターは ただのジャマもの (ということに北大構内走の最中に気づいた)． intrinsic Gaussian CAR normal なモデルは計算はカンタンなんだけど (条件つき局所確率は単純な正規分布にしたがうから)， その「実態」の理解は私などにはなかなか難しい (確率場全体を簡単な多変量確率分布で定義するのが難しそう) といったシロモノなのだが …… もしこの car.normal() で生成される GRF をひとつの多変量正規分布 (?) で書けるとすると， その「各地点の平均はことごとくゼロ」 と設定されているはずだ， とわかった (ただし分散共分散行列のほうはナゾで， 定義できないのではないかと思う …… まあこのあたり variogram と対応可能な GRF あたりから勉強しなおすべきかも)．
• で， 上の単純な GRF モデルで気になるのは， 気になるのは「zero inflated な地帯」ともうしますか， ひたすら -0-0-0-0-0-……と続いてる場所において (ひたすら -1-1-1-1-1-…… でも同様だろうが)， 状態変数が -10 から 0 ぐらいの範囲をぶんぶん 振動しまくっていること． これはまったく妥当な事後分布なのだが ……
• しかし， このあたりが原因 (?) で， モデルをちょっと複雑にしたときに収束しなくなるのかも， という気もする． 今度は Distance 依存性も推定させようとしたもの (BUGS code， 事後分布表)． 上の binomial なやつではまがりなりにも距離依存性を 推定できていたわけだが， こちらの GRF なやつだとかくのごとくコケてしまいます， と．

• 上のように計算しながら， 院生のデータ解析コンサル． R による作図と glm() によるやや複雑なモデリング …… 「あろめとりっく」なハナシなので， 説明変数を log.x としながら， family = gaussian(link = "log") とするわざ (今年 3 月の盛岡での GLM 自由集会 の久保投影資料を参照)．
• もちろん， こんなごまかしわざ (盛岡での粕谷さんのポスター発表で この「x と y に『誤差』がある場合の回帰の問題点」 が指摘されてましたね) ではなく階層ベイズモデル化してケリをつけたいところだが …… まあ， 院生たちはまだゐんばぐすは修行中， ということで． しかし， 推定結果と観測データを比較してみたら， やっぱり「個体差」がすごかった．
• 2145 研究室発． 2200 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • 北大構内走 45 分間． ラジオ体操． ストレッチ．
• [今日の食卓]
  • 朝 (0700): 米麦 0.6 合． キャベツ・ネギのサラダ．
  • 昼 (1230): 研究室． 食パン． チーズ． キャベツ・ネギのサラダ．
  • 晩 (2230): 米麦 0.8 合． キャベツ・ニンジン・エノキダケ・鶏肉の煮物． サバの電磁波酒蒸し． ヨーグルト．