ぎょーむ日誌 2007-09-08

2007 年 09 月 08 日 (土)

• 0800 起床． 朝飯． コーヒー． 洗濯． そうじ． 怠業． 1120 自宅発． 台風 9 号通過したあとの晴． 今日は気温たかくなるもよう． 1145 研究室着．
• そういえば昨日， 荒木飛呂彦先生のイラストが Cell の表紙に! というけっこうショッキングなニュースが駆けめぐりましたなぁ …… 牧先生とか発奮しておられるにちがいない． 私もせいぜい良い論文を書いて， 吉田戦車画伯に表紙を描いてもらおう．
• また ECOAS プランクトンデータの階層ベイズモデル化にとりくむ． 一晩ねてアタマがすっきりしたおかげなのか， BUGS コーディングの改善点がわかったような気がした …… で， 試してみるとうまくいった． 最終的な結果を得るためにはやはり 10000 MCMC step 超の計算が必要だろうけど， 混交は格段によくなっている．
• 改善点 (らしきもの) の要点はこうだ． 統計モデリングをするときに fixed effects と random effects のパラメーターをわけて考えるクセがついていると， こういうコーディングになりがちなんだけど

  	# coding 1
  	for (i in 1:N.x) {
  		...
  		x[i] <- beta0 + beta1 * z[i] + re[i]
  		re[i] ~ dnorm(0, tau.re)
  		...
  	}
  

  もしデータ構造が下のような書きかたを許すのであれば (つまりひとつの対象 i から 1 回だけデータを得た場合だ)，

  	# coding 2
  	for (i in 1:N.x) {
  		...
  		x[i] ~ dnorm(mu.x[i], tau.re)
  
  		mu.x[i] <- beta0 + beta1 * z[i]
  		...
  	}
  

  このように書いたほうが混交が改善されるように思う． 少なくとも， 計算速度はかなり改善される． ただしモデルとデータ構造によっては こういうカタチに書きなおすことができないので， 注意．
• …… しかし Trap 窓は忘れたころにやってくる …… と思ったらこれは単純なコトで， 負の値の log を計算させてた．
• 計算させつつ昼飯．
• 階層ベイズモデル， ぢりぢりと改善されつつある． パラメーター数が多いので定常状態に到達するまで かなり時間かかるのはしょうがない， のかしらん? …… ぢつは， ゐんばぐすを用いた計算でこの時間を短縮するために有効なワザは 「説明変数の 中央化 (説明変数の平均がゼロに近いこと)」 なのである． うん， なんだかいつも理不尽に思えてしまうんだけど， そうなっているのだ． そして統計モデル中で生成される説明変数 (隠れ変数) をうまくあつかわねばならぬ， というあたり， よくよく注意する必要あるわけで．
• しかしこのあたり， 小細工にハシりすぎると計算がどんどんおかしくなるわけで．
• 21000 MCMC step で 390 秒ぐらい， 混交はもうちょい …… 魚いる・いないモデルから「水温」なる説明変数をぬくか．
• …… といったコトどもが問題なのではなく， ここによーやくにして 「モデルが複雑になると MCMC 計算がふらふらしなくなるのか?」 がわかってきたよーな気がする． うーむ …… おそらく， なのだけど， 説明変数は可能なかぎり確率変数化せよ， といった原則があるように思えてきた．
• WinBUGS の中では変数は確率変数かそうでないか， のふたとおりある． 推定すべきパラメーター数が少ないときは， とくに何も考えずにモデリングしても， Gibbs sampler はしゃかしゃか走りまわってくれる．
• ところが， パラメーター数が増えて， 変数のあちこちに (交互作用ではなく) 相互作用があるような場合には …… これらの変数は (上述のように random effects と関連づける， など何かうまい口実をもうけて) 確率変数化する必要がある (BUGS 言語で言えば ~ で生成される変数)．
• というのも， 決定論的な変数 (<- で生成される変数) が説明変数としていくつも含まれているようなモデルは「カタい」と申しますか， あちこちにシバりがあって Gibbs sampler が自由には動きまわれぬ， と見えるんだよね．
• ともかく， 今回のモデリングに関してはこういうふうに変更してみると， 「どうしてココは動かんのだろねぇ」 というところがクビキをとかれたようにしゃかしゃかしはじめた． いやー， 今まで気づかんかった ……
• そして魚いる・いないモデリングを logistic な定式化から， ポアソン分布にもとづくもの (コーディングでは dcat(p.fish[i,]) をもちいる) に変更してみる． いる・いない問題だからといっても何でもかんでも logistic にすればいいってモンじゃないしね． 得られた結果もわかりやすくなってるように思う．
• いろいろ試行錯誤して， まっとうそうな計算結果得られた． ssh 経由で使ってる A801 の Dell desktop 機で 160 秒．

Inference for Bugs model at "/home/kubo/ecoasZP/winbugs/model.bug.txt", fit using winbugs,
 3 chains, each with 16000 iterations (first 10000 discarded), n.thin = 30
 n.sims = 600 iterations saved
                mean     sd    2.5%     25%     50%     75%   97.5% Rhat n.eff
beta.f[1]      -0.71   0.45   -1.53   -1.04   -0.73   -0.41    0.18 1.00   600
beta.f[2]       0.89   0.25    0.40    0.73    0.89    1.04    1.41 1.02   130
beta.f[3]      -0.16   0.15   -0.46   -0.24   -0.16   -0.06    0.11 1.01   450
beta.f[4]       1.34   0.13    1.10    1.26    1.34    1.43    1.57 1.01   600
beta.z[1]      -0.12   0.00   -0.12   -0.12   -0.12   -0.12   -0.12 1.00     1
beta.z[2]       2.43   0.00    2.43    2.43    2.43    2.43    2.43 1.07    64
beta.z[3]       0.50   0.36   -0.25    0.28    0.51    0.73    1.15 1.03   130
beta.z[4]       0.20   0.12   -0.02    0.12    0.21    0.28    0.43 1.03    94
beta.c[1]      -0.10   0.21   -0.59   -0.21   -0.08    0.04    0.21 1.02   170
beta.c[2]       0.64   0.54   -0.24    0.24    0.57    0.93    1.81 1.01   280
beta.c[3]       0.25   0.07    0.14    0.20    0.25    0.30    0.41 1.01   390
beta.c[4]       0.00   0.01   -0.02    0.00    0.00    0.01    0.02 1.04    66
beta.c[5]       2.72   0.72    1.33    2.29    2.75    3.14    4.18 1.01   190
tau.err[1]     30.19  33.70    1.13    6.32   17.88   41.49  127.30 1.00   600
tau.err[2]    941.65 829.95   93.97  341.90  683.70 1308.25 3069.92 1.00   600
tau.err[3]     54.72 175.42    0.05    0.24    2.75   19.75  665.88 1.02   110
tau.err[4]     54.34 223.75    0.00    0.00    0.11    8.58  537.04 1.21    14
tau.err[5]    150.44 334.49    3.28   14.22   41.55  129.35  910.75 1.00   600
tau.delta[1]    0.14   0.03    0.08    0.11    0.13    0.16    0.20 1.00   600
tau.delta[2]    0.33   0.15    0.17    0.24    0.29    0.37    0.64 1.02   420
tau.delta[3]    5.59  24.88    0.61    0.94    1.49    3.36   25.40 1.03   120
tau.delta[4]  108.36 258.86    0.51    1.19    5.65   75.45  907.25 1.03   110
deviance     -250.71 140.08 -495.42 -350.30 -263.80 -159.12    9.43 1.01   320

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

pD = 282.3 and DIC = 31.6 (using the rule, pD = Dbar-Dhat)
DIC is an estimate of expected predictive error (lower deviance is better).

• へろへろと撤退． 1850 研究室発． はてさて， 上のモデルは因果関係で 魚-動物プランクトン-植物プランクトンの 観測値間の関係を説明しようとしたわけだけど (BUGS code) …… いろいろと他にも改善すべきところはあるんだけど， じつはこのあたりこそが (例の呪われ Wishart 分布から生成されてしまった) 分散共分散行列に駆動される多変量正規分布な事前分布 (その平均値は種間相互作用ナシ条件化のもの) を使って生成するべきなのかも， という気がしてきた． つまりこのデータからは因果関係は推定不可能であり (原理的に)， 相関関係だけが推定可能なのかも， というハナシで …… 北 12 生協で買いもの． 1910 帰宅． 晩飯の準備． 晩飯．
• [今日の運動]
  • またうんどう休養日になってしまった
• [今日の食卓]
  • 朝 (0900): クルミパン． オレンジ． ナシ．
  • 昼 (1310): 研究室． 北 12 生協の納豆巻． ミディトマト．
  • 晩 (2030): 米麦 0.8 合． ネギ・豆腐の味噌汁． 煮魚． コマツナのゴマあえ． キャベツ・ニンジン・キュウリ・ネギのサラダ．