KuboLog 2007-09-18

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

2007 年 09 月 18 日 (火)

0910 起床 …… うーむ， ECOAS 下うけ労働による生活周期撹乱だ． 0925 自宅発．晴． 0945 研究室着．朝飯．コーヒー．

昨晩，帰りがけに命じておいた投げやり長時間計算， 30000 MCMC step に 3200 秒弱を費やして終了していた．やはり，というか，ほとんどのパラメーターはキレイに収束しているんだけど， deviance がいまだに変なオチかたをしていてキモチ悪い． tau.log.potential と deviance の比較から，おそらくゐんばぐすの Gibbs sampler ちゃんが random effects まわりの探索的なぢりぢり改良に熱中しているのだろう，と憶測される．いやはや，ってかんぢ．

いまさらながら zero-inflated 型のベイズモデルに変更してみるかな …… (いる・いない) + (バイオマス) の混合分布なモデルってことで．おそらくバイオマスに関する標高依存性は多くの species でスっ飛んでしまうだろうけど．実際のところ，観測データはかくのごとく zero-inflated なんで．

Zero-inflated 型へのモデルの改良， Trap 窓がでまくりで苦闘する．本日の教訓: 「いる・いない」問題だからといって，何でもかンでも logistic なモデルにすればいいってもんぢゃない．

つまり，潜在的バイオマス = 0 なら存在する確率はゼロにきまっているんだから，たとえば p = 1 - exp(-係数 * バイオマス) にする，と．これはぢつは先日の ECOAS プランクトンデータでも使ったモデリングなんだけど，自分のやったことを忘れていましたよ．

で，この「係数」とやらが常に非負でなければならぬ，と．そうか，こういう制約をあたえるなら logistic でもうまくいったかもしれんな．しかし， 1 - exp(...) のほうがパラメーターが少なくっていいや．

モデルを複雑にしてしまった罰として，バカみたいに計算時間が長くなった． 2000 MCMC step に 560 秒．たぶん， p = 1 - exp(...) などと確率を定義しちゃったから，ゐんばぐす内部で Gibbs sampling ができなくなって Metropolis 法にきりかわったんだろう．

計算まち時間， Nature ながめる．手もとにある最新号 (vol. 449 num. 7159) は集団・個体 level の植物生態学ものが三つのっている．どれも現象論と観念論の中間的なものだけど， scale-free network だとか scaling 則だとか持ちだして， Nature 誌面の occupancy に成功している．

(表紙にもなっている) その半乾燥地帯の植物のパッチサイズ分布 (例によってベキ分布) 生成する集団モデル論文のほうに対応する Solé による News & Views では昔なつかしのカタストロフ図なんて出てるし．やっぱ，皆さん，無限大集団の性質を知りたいのかしらん? というか，無限大集団でないと集団の性質を調べるのがめんどうってコトなんだろね．

ECOAS 八甲田湿原群バイオマスデータの統計モデリング，「観測誤差」に該当する部分の事前分布にやや強引な制約を与えることで収束よくなった …… 6000 MCMC step (およそ 2400 秒) の計算結果をみると，安全のため 8000 ぐらいに延伸して計算やりなおしたほうがよさそうだ．

計算まち時間にこの zero-inflated なモデルに生じているヘンな非対象性について検討しつづけている． Zero-inflated な統計モデルは「観測されたバイオマス (あるいは個体数でもなんでも) がゼロだからといってホントにその場にいないと断定していいの? サンプリング面積とか努力が足りないだけなんでは?」といった状況に対処するものだ．ここでちょっと気になる非対称性があり，

ある場所で観測されなかった生物のあつかい: これは「いるかいないかわからない」ということで zero-inflated なモデルでうまくあつかえそう
ある場所で観測されてしまった生物のあつかい: これはもう「いる」ということが確定してしまっていて，統計モデルの中の「いるかもしれないし，いないかもしれない」という不確定性の部分が影響してくるのは，じゃま

…… と整理してみると，やはりもうちょいモデルの修正が必要だな．えーい．「いるんだから『いる』んです」というふうに観測データ使わねばならぬ気色わるさよ．

ということで，ゐんばぐすのけったいな関数 step() を使って，統計モデルをこう定義する．

for (i in 1:N.quad) {
	for (j in 1:N.sp) {
		MassA[i, j] ~ dnorm(massA[i, j], tau[1])
		massA[i, j] <- occ[i, j] * bm[i, j]
		occ[i, j] <- step(Occurrence[i, j] + occurrence[i, j] - 1)
		# occurrence
		Occurrence[i, j] ~ dbern(prob[i, j]) # observation
		occurrence[i, j] ~ dbern(prob[i, j]) # hidden variable
		prob[i, j] <- 1 - exp(-exp(betaB[5]) * bm[i, j])
		# biomass after competition
		bm[i, j] <- alpha[i, j] * alpha[i, j] / sum.alpha[i]
	}
}
(以下略)

step(x) は x ≥ 0 (ゼロ以上であることに注意!) のときに 1 で，そうでないときはゼロを返す関数である．

ということで，計算を途中で止めて最初からやりなおし．

8000 MCMC step (6000 step burn-in) の計算結果 (事後分布表， BUGS code)． deviance の収束はよい．潜在変数 occurrence[,] を導入したご利益か? beta* については …… まあ，こんなものかな．

隠れ変数 occurrence[,] はどうでもよい …… と言いたいところだが，あとで事後分布を調べるときに「occurrence[,] == 1 という条件のもとでの事後分布は?」といった計算をしなければならないかもしれないので，とりあえず出力させているだけである．

うーむ， ECOAS 的な …… と書いていて ``ECOAS'' が何の意味なのかわからなくなった．驚くべきことに Google 検索しても自分の HDD 検索しても，まったく意味がわからん．環境省の地球環境研究総合推進費なる予算の平成 17 年度新規採択研究課題であることにはまちがいないんだが …… ECOAS っていまや obsolete になったお役所用語か何かだったのか?

まあ，いいや …… ともかく何か ECOAS 的な (「環境省のヒモつき研究にありがちな」ぐらいの意) 推定計算みたいなことをやるときには，各コドラートの random effects の事後分布もやっぱり出力しといたほうが便利だね …… ということでまた計算やりなおし． 55 分ほど要するはず．

さて，今月はじめごろから着手した ECOAS データ (高山湖沼のプランクトン群集 + 八甲田湿原群の植物群集) の統計モデリング，とりあえずひとくぎりつきつつあるように思う …… 「とりあえずひとくぎり」とは「簡単にできそうなコトだけ (ECOAS 的な最低限の必要条件にもおそらく留意しつつ) やったけど，その先にある面倒には手をだしてない，やったらタイヘンなことになるから」状況とでもいいますか．

まあ，タイヘンといっても WinBUGS ではもう計算できないような統計モデリングになるってことぐらいで …… プランクトンに関してはもう少し WinBUGS で押せるかも．いずれにせよ，当方のすけぢゅーりんぐ破綻を回避するべく，これ以上に発展したモデリングはやらない，と．

…… あ，いかん，ベイズ統計学本などながめて，より面倒な問題を解決する方法など検討していた．撤退てったーい．

ECOAS PukiWiki 上にはごく簡単な計算報告はいろいろアップロードしているが (本日のぶんはまだだけど) …… この wiki 上の解説を詳しくするべきか，それともまた別に報告書を書くべきか．

そういや，パラメーターなんかの事後分布は推定したけど，これを使った予測計算なんてのはぜんぜんやってないな．というか事後分布図なんかも作ってないし．

などとうだうだやってるうちに計算終了 (事後分布表， BUGS code は先ほどと同じ)．なぜか beta* の収束が先ほどより良い．まあ，実質的には何もかわっていないんだろうけど．

ちゃんとした報告を書くのはしんどいので，とりあえず PukiWiki 上に図とか説明文とか並べてみることにするか …… 近ごろはゐんばぐす専用機と化しつつある A801 院生部屋の Dell 機からは荷物はまとめて仮想的にというか ssh 的に撤収．

やっぱり図は描いてみるもんだなぁ …… このモデルでは種間競争がキビしすぎるみたい，とわかった．モデリング & 計算やりなおしになりそうだな．

[キビしすぎる競争??]

この観測データ中，もっとも現存量の多い植物のひとつヌマガヤ (spc23) の例．横軸は「単独種で実現したバイオマス」でタテ軸は「現場で観測されたバイオマス」．たとえば他種がいなければヌマガヤだけで 1000 g 超になるはずだったプロットで， (他種との競争の結果) 現実には 382 g でオワった (ちなみにこのコドラートの合計バイオマスは 694 g)．いくら何でもこのコドラートの最優占種がここまで競争疲れでヤセてしまうものだろうか? 二番手三番手は spc36 (217 g) と spc01 (68 g)．

上の例，数値的にはあってる，といいますか，たしかにこうなるように (現象論的な) 競争モデルを定義した．

改善策あれこれアタマに去来するのだが，どれも計算時間がさらに長くなりそうなモノばかりで …… 何もせぬまま撤退． 1910 研究室発． 1930 帰宅．晩飯の準備．晩飯．

[今日の運動]

今日もうんどう休養日

[今日の食卓]

朝 (1000): 研究室．セイコーマートのラスク．リンゴ．
昼 (1300): 研究室．セイコーマートのラスク．ヨーグルト．リンゴ．
晩 (2030): 米麦 0.7 合．ニラいり麻婆豆腐．ダイコンの中華ふうサラダ．

ぎょーむ日誌 2007-09-18

2007 年 09 月 18 日 (火)