ぎょーむ日誌 2007-09-(11-20)
2007 年 09 月 11 日 (火)
-
0855 起床．
だめだ．
ぷらんくとんばて．
0910 自宅発．
曇．
0925 研究室着．
朝飯．
コーヒー．
-
なぜか PukiWiki を更新する作業，
とか．
-
ECOAS プランクトンデータ解析したうけ，
動物プランクトン作図のつづき．
なーンかヘンだな
……
と思ったら
……
札幌市郊外の空沼 (標高 1460 m)，
動物プランクトンの現存量ゼロ!
なのか．
植物プランクトンはいるのにね．
-
昨日，
species ごとの plot というの作ってみたんだけど，
よくわからんかった．
そこで，
動物プランクトンのバイオマス組成が環境によってどう変わるか，
の R 作図．
右の列は不名誉なる割算値．
しかし，
群集データって何とも難解だねえ．
-
やっぱりよくわからんので，
環境経度にそって「だいたい等サイズ」
になるよう自動グループわけしてみる．
この作図プログラミングは
ちょっとばかり
めんどうだった．
> x <- tapply(b$biomass, b$spc, sum)
> x
spc01 spc02 spc03 spc04 spc05 spc06 spc07 spc08
1.7361e+03 4.6890e-01 2.9974e+02 1.6278e+01 6.6488e+01 6.6609e-03 5.9887e+02 4.2688e-02
spc09 spc10 spc11 spc12 spc13 spc14 spc15 spc16
2.5548e+01 4.0857e+00 6.4783e+00 8.2282e-03 1.9133e+02 1.3868e+01 1.2451e+01 1.9800e+01
spc17 spc18 spc19 spc20 spc21 spc22 spc23 spc24
2.0878e+02 2.4613e-02 4.8735e-02 3.0160e-03 1.9432e+03 6.4767e+00 8.2480e+01 5.8819e-01
spc25 spc26 spc27 spc28 spc29 spc30 spc31 spc32
3.6433e-01 5.0648e+00 2.0127e+01 2.1812e+02 2.1050e+02 6.7776e+01 2.3986e+01 8.1380e-02
spc33 spc34 spc35 spc36 spc37 spc38 spc39 spc40
6.1722e-01 2.7188e+00 6.0130e+02 5.6498e-01 1.7581e+00 2.0585e-01 6.8390e+00 1.2842e+01
spc41 spc42 spc43 spc44 spc45 spc46 spc47 spc48
1.5507e-04 1.2979e-02 4.8535e-01 2.8068e-04 2.7740e-04 8.7592e-02 1.0852e-03 8.0166e-02
spc49 spc50 spc51 spc52 spc53 spc54 spc55 spc56
1.2452e+00 1.9715e-02 8.5896e-03 1.2170e-03 1.6018e-01 1.5544e+00 2.6706e-01 1.1038e-01
spc57 spc58 spc59 spc60 spc61 spc62 spc63 spc64
1.4924e-04 1.4111e-01 1.3283e-04 1.2423e-04 4.7581e-03 2.9081e-05 2.0023e-03 9.6510e-03
> x[x > 10]
spc01 spc03 spc04 spc05 spc07 spc09 spc13 spc14 spc15 spc16
1736.133 299.740 16.278 66.488 598.868 25.548 191.333 13.868 12.451 19.800
spc17 spc21 spc23 spc27 spc28 spc29 spc30 spc31 spc35 spc40
208.779 1943.213 82.480 20.127 218.120 210.497 67.776 23.986 601.300 12.842
> unique(b[b$spc %in% names(x[x > 10]), "spc.name"])
[1] D.dentifera D.galeata D.cristata Ceriodaphnia Bosmina
[6] Chydorus Diaphanosoma Sida Leptodora Polyphemus
[11] Holopedium Acanthodiaptomus Eodiaptomus Megacyclops Thermocyclops
[16] Cyclops Diacyclops Mesocyclops Nauprius Asplanchna
ふーむ，
20 種になったか
……
などと検討してると，
屋久島めんどう作図問題が
……
library(grid)
の復習をして，
なんとか作図プログラムを変更する．
とりあえず 20 種の動物プランクトンを対象とする階層ベイズモデル作り．
さて，
ここで昨日解明した WinBUGS の
dmnorm()
と
diwsh()
を駆使して coding していくわけだが
……
はてさて，
この 20x20 の分散共分散行列の Gibbs sampling
ってどれぐらい時間かかるものなのかしらん?
……
わずか 300 MCMC step 計算させるのに 880 秒もかかった!
うーむ，
これが，
これがゐんばぐすの限界なのか?
だとしたら，
あとさき考えずに
「とにかく種数を増やせ!」
なる自己破滅的なスローガンをかかげて
何かに取り憑かれたかのようにひたすら何百種も
サンプリングしたがるあのアタマのおかしい群集生態学者たちのデータの解析，
こういう道具だてでは無理ってことなのかしらん?
よくわからんけど，
R2WinBUGS
戦闘ドクトリンの限界をさぐるべく，
16000 MCMC step の計算を命じる．
およそ 13 時間後に終了する予定．
ちょー長時間計算命じたものの，
あまりうまくいくとは期待してない．
種間のほとんどの組みあわせにおいて相関ゼロだろうし．
いやはや，
何たる推定計算だろうねえ
……
1930 研究室発．
帰宅途中に
「ゐんばぐすがノロい理由，
ひょっとしたら推定計算には必要のない 20x20 逆行列計算
で苦闘してるのかしらん?」
と考えついて思わず研究室にもどりかけたけど，
「まー，いっか」
と放置することに．
1945 帰宅．
体重 68.6 kg．
うんどう．
晩飯．
[今日の運動]
[今日の食卓]
- 朝 (0930):
研究室．
豆パン．
- 昼 (1425):
研究室．
豆パン．
リンゴ．
早くも当地では 1 個 100 円リンゴがでまわっているので，
一日 2-3 個たべてます．
- 晩 (2200):
米麦 0.8 合．
キムチ・マイタケ・卵の炒めもの．
コマツナのゴマあえ．
キュウリ．
2007 年 09 月 12 日 (水)
-
0850 起床．
うーむ，
あいかわらず，
やられているな
……
朝飯．
コーヒー．
1020 自宅発．
曇．
1035 研究室着．
-
13 時間かかるとおもってた ECOAS 主要プランクトン
20 種データの階層ベイズモデルの計算，
実際のところは 12 時間で終了してた．
収束はぜんぜんよくない
……
やはり 20x20 の分散共分散行列の推定はこのゐんばぐすでは無理なのか，
と考えたんだけど，
じつは当方も一ヶ所モデリングをしくじってることに気づいた．
これって二重に random effects が入ってる部分がある!
-
……
と思ったけど，
モデルをよくよく整理してみると，
log(密度) = (池の効果) + (種の効果) + (池・種の効果)
とモデリングしてよい，
とわかった．
あ，
でも他の環境要因な線形予測子とは区別して書かんといかんわ．
-
ゐんばぐすの意味不明えらー対策．
初期設定ファイル inits*.txt の読みこみで
expected c operator
うんぬんとでてきて計算続行してたら，
これは初期化の部分で一次元の vector をラベルつきで渡している，
ということ．
この場合は初期値 vector に重畳にも as.vector(...)
つけると解決する．
-
昨晩帰宅時に気にしてた「よけーな逆行列計算」
あってもなくても計算時間にはほとんど影響しない
(わずか 300 MCMC step に 880 秒ぐらい)，
とわかった．
20x20 の逆行列計算ぐらいは数百分の 1 秒ぐらいでできる，
ということなのか．
-
ちなみに分散共分散行列はヤメて，
動物プランクトン種間に相関がなく random effects がある，
というモデルに変更すると 300 MCMC step の計算はわずか
24 秒
で終了しちゃうんだよね．
いやはや．
-
ハラがたつので，
この無相関モデルで少し計算してみるか．
-
12000 MCMC step の計算やらせても
640 秒
で終了する．
それはいいんだけど，
収束がむちゃくちゃだ．
うーむ
……
-
それから夜まで階層ベイズモデリングの試行錯誤
ひたすらひたすら続く．
モデルはぢりぢりと改良されてて収束も良くなってはいるんだけど，
数千 MCMC step ぐらいではなかなかキレイにそろってくれない
……
-
ああ，
なんか，
私にとっては
こういう「水中」「組成」モノはけっこう鬼門だなぁ．
推定がしんどい理由はおそらく，
サンプルごとの random effects が大きい，
にもかかわらず (ホントの) 測定時誤差の統計モデリング
をどうすべきなのかよくわからない，
これらの区別が難しくて MCMC 計算がなかなか定常状態に到達しない，
といったあたりか?
-
そして東北地方からはいろいろとメイルが
……
-
へろへろと撤退．
2000 研究室発．
ばてたので外で晩飯．
2130 帰宅．
-
R
こんさる．
apply()
系関数と文字列操作は R 中級者にむけての昇段試験だねえ．
> # たとえばこういう data.frame があったとする
> d
x y
1 1.6 1.1
2 1.1 0.6
3 0.4 -0.7
4 -0.3 0.8
5 0.2 -0.1
> # apply() と sprintf() を組みあわせて新しい label 列が作れる
> apply(d, 1, function(r) sprintf("xy=%g:%g", r["x"], r["y"]))
[1] "xy=1.6:1.1" "xy=1.1:0.6" "xy=0.4:-0.7" "xy=-0.3:0.8" "xy=0.2:-0.1"
> # あるいは paste() を使ってもよい
> d$label <- apply(d, 1, function(r) paste("xy", r["x"], r["y"]))
> d
x y label
1 1.6 1.1 xy 1.6 1.1
2 1.1 0.6 xy 1.1 0.6
3 0.4 -0.7 xy 0.4 -0.7
4 -0.3 0.8 xy -0.3 0.8
5 0.2 -0.1 xy 0.2 -0.1
だめだ，
夜は仕事がススまん
……
[今日の運動]
[今日の食卓]
- 朝 (0940):
米麦 0.6 合．
ワカメ・豆腐の味噌汁．
ブナシメジ・卵の炒めもの．
- 昼 (1530):
研究室．
米麦 0.6 合．
クノールのカップカレースープ．
- 晩 (2030):
グロビュールのスープカレー．
ここの自家製ソーセイジ
(豚肉あるいはエビ)
がたいへんおいしかった
……
2007 年 09 月 13 日 (木)
多様性生物学基礎論（生物多様性論I）
場所：地球環境科学研究院：C−１０４室（去年と同じです）
時間: ５限目（16：30〜18：00）
１ 陸上植物（高橋英樹：博物館） ２ 回（10月1，3日）
２ 土壌微生物との相互関係（春木雅寛：地球環境）１ 回（10月10日）
３ 藻類（小亀一弘：理学） １ 回（10月15日）
４ 陸上動物
哺乳類（鈴木仁：地球環境） １ 回（10月17日）
鳥類（綿貫豊：水産） １ 回（10月22日）
魚類（前川光司：北方圏フィールドセンター）１ 回（10月24日）
昆虫（大原昌宏：博物館） ２ 回（10月29，31日）
５ 陸水生物（岩熊敏夫：地球環境） １ 回（11月5日）
６ 海洋生物（馬渡駿介：理学） ２ 回（11月7，12日）
７ 遺伝的多様性（木村正人：地球環境） ２ 回（11月14，19日）
生態学基礎論 （生物多様性論II）
１ 多種共存機構
植物（高田壮則：地球環境） ２回（11月26，28日）
動物（野田隆史：地球環境） ２回（12月3，4日）
２ 進化と種分化
植物・藻類（堀口健雄：理学） １回 （12月10日）
昆虫１（片倉晴雄：理学） １回 （12月12日）
昆虫２（秋元信一：農学） １回 （12月17日）
哺乳類（増田隆一：先端研） １回 （12月19日）
３ 進化と発生
脊椎動物（栃内新：理学） １回 （1月7日）
昆虫（三浦徹：地球環境） １回 （1月9日）
４ 生物の絶滅（大原雅：地球環境） １回 （1月16日）
５ 生物多様性解析法：統計モデリングの基礎（久保拓弥：地球環境）
２回 （1月21，23日）
八甲田湿原メイル，
そういえば ECOAS PukiWiki 上にも整理してたなぁ，
と読みかえしてみる．
「共存」
(coexistence)
ってよく考えずに使われてる生態学用語の一例だよねえ．
端的にいって，
共存なるものはひどく観念論的な何かで，
現実ときちんと対応づけることができない．
観測データをいくら解析しても「共存」とやらは見えてこない．
では共存はどこにあるのか?
というとそれはすごく抽象化された
数理モデルの中にだけある概念
なんだよね．
無限の個体数をもつ集団
(無限集団)
の数理モデルの世界なら共存は定義可能．
にもかかわらず，
そこらへんに
A と B という生物がいっしょにいるのを見かけたら，
「A と B は共存してる」
とか安易に言っちゃうヒトもけっこういるわけで．
ともあれ，
アタマを少し層別刈り取り植物群落データにもどして，
メイルかきとか．
夕方までメイルやりとりあれこれが続く．
日没時刻ぐらいから ECOAS PukiWiki 上にうだうだと
動物プランクトン「組成」モデルの報告かき．
ひととーり書けたのでアップロードしてメイル連絡して撤退．
1925 研究室発．
1940 帰宅．
体重 69.0 kg．
ぷらんくとん太り?
うんどう．
晩飯．
[今日の運動]
[今日の食卓]
- 朝 (0900):
米麦 0.6 合．
ブナシメジ・ワカメ・豆腐の味噌汁．
- 昼 (1320):
研究室お茶部屋．
米麦 0.6 合．
レトルトパウチドの麻婆豆腐．
- 晩 (2220):
米麦 0.7 合．
ブナシメジ・ワカメ・豆腐の味噌汁．
サンマのショウガ煮．
2007 年 09 月 14 日 (金)
-
0710 起床．
ねむい．
朝飯．
コーヒー．
0930 自宅発．
晴．
0945 研究室着．
-
神山さんから八甲田湿原群の植物群集データおくっていただく．
すでに修論の段階で R をつかったデータ解析しておられるので，
データ読みこみは簡単にできた．
いつでもこの段階で苦闘するので，
すごく助かる．
-
さっさとこのデータの作図作業をススめたいところなんだが
……
わけわからん再査読依頼がまた来ている．
これもまた，
こんなヘンな論文なんか通さなくてもいいのに，
というものだ．
ああ，
どうしようかねえ．
-
イヤな案件の蓄積を圧縮するべく査読ぎょーむにとりくむ．
-
昼飯をはさんで査読報告かき．
うう．
1453 査読報告提出終了．
いやはや．
-
A802 院生部屋の廊下がわのドアのしまらなくなった，
とのことで調べてみたんだが
……
どうも
door closer
(wikipedia,
google image)
が不調のよーで．
アームのねじをゆるめる方向に回したら，
とりあえずしまるようになった．
-
イヤなことセン滅ドクトリンを徹底するべく
(いつのまにか過激な方向に!?)，
もいっこのイヤな再査読にもとりくんでみるか．
-
と思ったら数セミ記事 (階層ベイズモデル解説 6 ペイジ)
の校正がきたのでそちらから着手．
-
それがひと区切りついたと思ったら，
また別のメイルやりとり．
-
お，
先日の
apply() + はりつけわざ
に関して
RjpWiki
では
apply()
なしの paste()
だけわざが紹介されてるな
……
じつは私はこういう「異なる長さの vector を引数にする関数よびだし」
をあまり理解してないというか，
理解はしてると思うけどなかなか使えないんだよねえ．
まあ，
R の場合，
こういうふうにお手軽に実験できるんで，
不審に思ったらすぐに挙動を調べればいいだけなんだけど．
-
再査読などイヤなこと抹消する元気もなくなってきたので，
R 上で八甲田湿原群落データをながめてみたり．
ふーむ
……
じつは統計モデリングも何も，
対数の世界でみたらだいたいもう直線のようにみえるんだけどねえ
……
-
だとすると，
この植物群集データのモデリングって，
高山湖沼プランクトンデータと同じく単なる群集バイオマスと
その組成の統計モデルの構築になってしまうのかしらん?
だとしたらハナシは比較的簡単なんだけど，
どうなんだろう?
1935 研究室発．
1950 帰宅．
晩飯の準備．
晩飯．
-
電気代 1728 円支払う．
いつもより 300-400 円ぐらい高いような．
新型洗濯機のせいかしらん?
静粛かつ powerful なやつなんだけど，
いささか電気を食うのかもね．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0830):
米麦 0.7 合．
ブナシメジ・ワカメ・豆腐の味噌汁．
サンマのショウガ煮．
- 昼 (1310):
研究室お茶部屋．
米麦 0.5 合．
フリーズドライ野菜スープ．
リンゴ．
- 晩 (2040):
米麦 0.8 合．
キャベツ・ニンジン・ニラ・ショウガ・ニンニク・
干しエビ・豚肉の炒めもの．
2007 年 09 月 15 日 (土)
-
0900 起床．
うう．
朝飯．
コーヒー．
洗濯．
怠業．
-
いつのまにか昼すぎだったので昼飯．
-
ECOAS 八甲田湿原群の植物群集データにとりくむ
……
が送っていただいたデータを R でいろいろとひねくりまわしていると，
何だかヘンなところに気づいた．
質問メイルおくって，
解析は中断．
-
とはいえ何もしないのもナンなので，
data.frame を再整備する R プログラムを書き
(問題のある種名は一時的にてきとーに改竄)，
library(lattice)
の
dotplot()
でてきとーなる図を作ってみる．
うーむ
……
「標高差」，
ねえ?
-
雨ふってるので研究室にいく気になれず．
買いもの．
晩飯の準備．
晩飯．
-
夜はひたすら八甲田湿原群の植物群落モデリングについて検討．
ECOAS 的なオトシどころとしては，
植物であれプランクトンであれ，
気温変化があった場合に群集の組成や
系のバイオマスの増減がどう応答するか，
に回答できればよいわけで
……
と考えると，
この湿原植物データの統計モデリングにおいて，
「獲得光量 (← モデルによる瞬間値の予測)) /
その植物種のバイオマス」
なる割算値指標はとくに役にたたないわけで
……
-
競争の結果として実現する種ごとのバイオマスは
Gibbs 分布にしたがう，
というベイズモデル straightforward の世界，
だよな．
競争力なるモノというか隠れ変数 (hidden status variable) が，
気温とかそういったあれこれに依存していて，
個体差とか調査地ごとの random effects とかにふりまわされている，
と考えればよい．
-
てなふうに考えがマトまってきた，
と思ったらとっくの昔に真夜中すぎてましたよ．
ECOAS は生活周期をめちゃくちゃにしてくれる．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0940):
イソップベイカリーのクロワッサン．
キャベツ・ニンジン・ニラ・ショウガ・ニンニク・
干しエビ・豚肉の炒めもの．
キャベツ・ニンジン・ネギ・パセリ・ハムのサラダ．
- 昼 (1330):
スパゲッティー．
レトルトパウチドのカルボナーラソース．
イソップベイカリーのバゲット．
キャベツ・ニンジン・ネギ・パセリ・ハムのサラダ．
- 晩 (1840):
米麦 0.8 合．
ネギ・ブナシメジの味噌汁．
サンマのショウガ煮．
コマツナのゴマあえ．
キャベツ・ニンジン・ネギ・パセリ・ハムのサラダ．
2007 年 09 月 16 日 (日)
-
0740 起床．
朝飯．
コーヒー．
デイパックの背負いバンドの修理．
0935 自宅発．
曇．
雨はやんでる．
0950 研究室着．
-
ECOAS 八甲田湿原群の植物データ解析のつづき．
昨日，
なぜかコドラート内で重複しているナゾの植物種コード
mikaya
でナヤんでいたんだが，
神山さんに教えていただいたところ，
これは、非常に見分けるのが困難な2種（ミヤマイヌノハナヒゲとミカヅキグサ）あるためです。解析に用いるデータとしては適切でなかったかもしれません。この2種が一緒に出現するコドラートでは、穂がついているものは2種の見分けがつきますが、付いていないものは分からないので、「ミヤマイヌノハナヒゲ」「ミカヅキグサ」「ミヤマイヌノハナヒゲかミカヅキグサ」の3タイプに分けてサンプル処理をしました。しかし、非常に似ているため、これらをまとめて「mikaya」としました。
とのことなので，
私の独断で乱暴にもまとめてしまうことに．
-
そしてデータ全体をながめるいろいろな作図作業がつづく．
図たちは ECOAS PukiWiki にアップロード．
-
この多種データからはそれほどコマかいことが
言えない
だろうな，
と見当がついてきたので，
ごく大雑把な
(パラメーター数が少ない)
階層ベイズモデリングから始めるとするか．
-
コーディングはすぐにできるんだけど，
ナゾのゐんばぐすエラーに苦闘する．
まあ，
こういう場合は，
ほぼまちがいなく，
当方の統計モデリングがかなりオカしいわけだが．
-
昼飯たべつつ調節つづく．
ECOAS プランクトンデータの場合と同じく，
random effects がかなり大きいのに
(ホントの)
測定誤差と区別がつかん!
というあたり．
-
どうやってもうまくいかんので，
とりあえずバイオマスの測定誤差の分散逆数を 10
などと固定してみる．
-
試行錯誤計算ススめつつ，
この統計モデリングのメイル書き．
毎度のことながら，
こういう説明メイル書きは時間がかかる．
-
ほとんどのパラメーターはだいたい数千 MCMC step
で収束してるように見えるんだけど，
deviance の減少に時間かかる．
やっぱり，
モデルがどこかまずいんだろうなぁ．
-
とりあえず 22000 MCMC step の計算結果
(事後分布表,
BUGS code)．
計算時間は，
えーとたしか 20 分ぐらい．
うーむ，
一番下の deviance がびみょー
……
-
なーんか理由があるはず，
なんだけどよくわからない．
-
わからぬまま撤退．
2000 研究室発．
2020 帰宅．
うんどう．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0830):
米麦 0.6 合．
キャベツ・ニンジン・ネギ・パセリ・ハムのサラダ．
- 昼 (1320):
研究室．
イソップベイカリーのバゲット．
ヨーグルト．
リンゴ．
- 晩 (2200):
米麦 0.6 合．
ネギ・ブナシメジ・卵炒飯．
コマツナのゴマあえ．
豆腐．
2007 年 09 月 17 日 (月)
-
0740 起床．
朝飯．
コーヒー．
1000 自宅発
……
しかし廊下の鏡をみるとヒゲがむさくるしかったんで，
いったん戻ってひげそり．
ここ数日，
カミソリまけしてたので放置してたんだけど，
あまりよくなってないような．
1010 自宅発．
曇．
1025 研究室着．
-
ECOAS 八甲田湿原植物データ解析のつづき．
さっさと「脱」割算値計算やれ，
という指示がきたので，
バイオマス分布問題はあとまわしにして，
「脱」割算にとりくむ．
-
これは基本的に簡単すぎる問題なんだが
……
とうっかりしたまま BUGS coding などやってると，
ついついまちがえてしまった．
ここまでの ECOAS データ解析 (動物プランクトンも含む)
は
biomass[調査値, species]
などなどと指定するデータ構造だったんで，
今回も惰性でそうしちゃったんだけど，
今回はサンプルごとの通し番号でやらねばならぬのだった
(「欠測」 species 処理のため)．
-
で，
割算値 4 個ばかりを「脱」割算な階層ベイズモデルで
自動的にどんどん計算させるしくみを作る．
さいわいなことにひとつの計算は 4000 MCMC step
(10 秒)
程度で終了するので 1 分もかからずに全 MCMC 計算終了，
と．
-
ひとやすみして昼飯
……
の前にお茶部屋かたづけをせねばならぬようだ．
昼飯．
-
計算結果のアップロード，
説明かき，
メイルかき
……
で一区切りついたのが 1500 すぎ．
-
なんだかよくわからない別件のデータ解析こんさるメイルがきてたので，
返事かいてみる．
-
ひと休み．
ふう．
A 棟 8F，
無人状態で平和だ．
-
ふーむ，
R
の package に
library(bootstepAIC)
なるものがあるな
(CRAN,
help)．
AIC のばらつきを評価するもの．
-
また八甲田湿原バイオマスの統計モデリングに復帰．
種間競争の関して，
もう一個パラメーターが必要，
という気がしてきた．
現象論的な種間競争モデルなんだけど
……
-
試行錯誤のあいまに断続的な計算まち時間が生じてしまう．
Tilman とかのやってた実験って，
種数が多いほど生産力だったか現存量 (収量)
が高くなるといった結果だったっけ?
うーむ，
そういう「すきまをうめる」ようなハナシって，
どうモデリングしたらよいのかわからんなぁ．
-
……
はっと気づくと，
手がとまっていた．
どうもアタマがばててきたよーなので，
砂糖いり紅茶を飲んで休憩．
すでに夜．
-
いろいろ試してみたんだけど
……
どうも非負のベキ数なんかを Gibbs sampling すんのは難しいね．
この案，
しばらく廃棄．
-
なげやりな長時間計算を命じておいて，
2005 研究室発．
雨．
2020 帰宅．
体重 68.4 kg．
-
リンゴなど食って休憩してれば，
うんどうする気になるかしらん，
と思ったんだけど
……
どうもアタマがばててるのか，
そういう気になれづ．
あきらめて晩飯の準備．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0900):
米麦 0.6 合．
ワカメ・煮干の味噌汁．
キムチ．
- 昼 (1310):
研究室お茶部屋．
米麦 0.6 合．
ワカメスープ．
リンゴ．
- 晩 (2130):
米麦 0.7 合．
ネギ・ニンニク・卵炒飯．
ワカメ・煮干の味噌汁．
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 月 19 日 (水)
-
0730 起床．
朝飯．
コーヒー．
湿原植物の競争モデルについての再検討，
いろいろと．
いくつか簡単な定式化の方法は思いつくんだけど
……
0930 自宅発．
曇．
0945 研究室着．
-
さっそく ECOAS 八甲田湿原群のバイオマス分布モデルの改善．
たとえば，
「自分以上の potential をもつ species にだけ悪影響をうける」
という単純化した競争モデルをゐんばぐすの
step()
で単純に実装してみると
……
やはり，
というか計算時間が
(冗談ではなく)
むちゃくちゃ
に長くなってしまいましたよ．
200 MCMC step に 1942 秒!
……
といった次第で，
三次元配列などというバカげた実装はたちまち却下される．
-
BUGS 言語の限界か
……
これはプログラミング言語のようなみかけをしているけれど，
実際のところは
「loop を展開して数式を並べる」
だけの機能しかないからなぁ．
-
また「potential のベキ」の方向で探索してみるか
……
みこみウスいけど．
-
2000 MCMC step 1130 秒でオワった．
収束はまだまだだけど，
もうちょい，
このあたり調べてみるか，
という気分にはなった．
potentialpower
で power > 1 を強制してるんだけど
(すると優占種が有利になる)，
けっこう差がつく方向に流れはじめている．
-
また計算まち時間．
今回は 50 分ちょいか?
-
もしこれでうまくいかなかったら
……
そうだなぁ，
また別の策としては，
ちょっとバカげたかんぢではあるけれど，
個体群ではなく群落レヴェルで「収量一定の法則」
みたいなのが (局所的に!) 成立していて，
群落バイオマスはその種組成には「あまり」依存しない，
というアイデアはどうかしらん?
そうすると，
先日のプランクトンモデルと同じく，
その限られたバイオマスを植物たちがあたかもぶんどりあいをしている，
といった描像で
……
まあ，
これならわれらがマイティーゐんばぐすでも問題なく計算できるよね．
-
ゐんばぐすは便利だが，
こいつの性能限界に生態学データ解析の
統計モデリングの発想が制約されるのは何ともはや，
な状況だ
……
ということで，
MCMC 本などながめてみる．
ふーむ，
slice sampler か
……
アルゴリズムわかったけど，
これって既知の確率分布でないときは計算たいへんそう．
-
HBC
(Hierarchical Bayes Compiler)
って Haskell + C で書かれているんだ．
それはいいんだけど，
これって複雑な階層ベイズモデルを定義できるのかしらん．
演算子
``
<-''
がないような気がする．
-
3300 秒ほどで上記の 6000 MCMCM step 計算終了したんだが
……
うーむ，
昨日よりはマシになったとはいえ，
種間競争は依然としてキビしい．
[まだキビしすぎる競争??]
またヌマガヤ (spc23) の例．
たとえば他種がいなければヌマガヤだけで
800 g 超になるはずだったプロットで，
(他種との競争の結果)
現実には 382 g でオワった (白丸)．
ちなみにこのコドラートの合計バイオマスは 694 g (黒丸)．
-
で，
A 棟 8F の湿原植物専門家たる小山さん・江川さん
(西村さんは不在だった)
におそるおそるお尋ねしてみたところ
……
久保質問が不明瞭で何とも回答に困るんだけど，
いくらなんでもヌマガヤだけになったからといって
現存量 600 → 800 g はヘンでしょう
……
とのことで，
がっくりしつつすごすごと撤退してきた．
-
ということで，
モデルをバイオマス「ぶんどりあい」型に変更して計算やりなおし．
その間に昼飯．
-
うう，
挙動もよくみないでいきなり長時間計算してもうまくいかないわけで
……
いろいろと統計モデルの改善のための試行錯誤つづく．
-
どうも
latent な
「いる・いない」変数のあつかいがマズくて，
MCMC 計算を乱しているらしい，
とわかってきた．
そのあたり修正してまた 2000 MCMC step の試運転．
収束のみこみが見えてきた!
計算時間も 620 → 500 秒と短縮されたし．
6000 MCMC step 計算を命じる．
25 分ぐらいかな．
-
次の生態学会大会
(ESJ55)
の案内くる．
3/14 (金) - 3/17 (月)，
最後の日が懇親会，
と今までとはちょっと異なる．
しめきりあれこれ:
-
企画集会，自由集会申し込み 10/26 (金) 17 時
-
一般講演申し込み 11/16 (金) 17 時
-
講演要旨登録 2008 年 1/8 (火) 17 時
なんだか気が重い．
ポスター発表は「アリの噛みつき確率」
かなあ
……
自由集会
はどうしよう．
またベイズものになるんだろうけど．
うーん，
やっぱりベイズの基本的なところをやったほうがいいのかしらん?
-
今回は粕谷さんが大会の実行委員長なんで，
すごくおいそがしそうな気がする．
-
6000 MCMC step 計算，
だいたいいいかんじなんだけど，
パラメーターのうちいくつかは 5000 あたりから定常になってる
……
ということで 8000 MCMC step
(burn-in 6000 MCMC step)
に延伸してやりなおし．
二重にコドラート random effects が入ってるので
(コドラートのバイオマス合計と各 species のぶんどり能力)，
苦闘してるのかしらん?
-
1800 秒弱で終了．
収束は問題ナシ
……
しかし統計モデルがダメだ．
標高によって変化するはずの競争力に関する種差がまったくでない．
そのあたりのめんどうをぜんぶ random effects に押しつけてるな
……
-
よくよく BUGS code を見直してみると，
これは二重 random effects うんぬん以前の問題として，
私の定式化がマズかった，
と気づかされた．
例によって「いる・いない」問題で，
やはりここは観察された
Occurrence[i, j]
と latent な occurrence[i, j] の確率を区別して
dbern(...)
せんといかんわけね
……
じゃないと，
species ごとの
「もうこの高度では生きていけません」
がうまく表現できないわけだ．
植物がどうやって標高を感知してるのかはまったくのナゾ
ではあるけれど．
-
そのあたり修正して，
また再計算を命じる．
-
ダメだ．
種差はでてるけど，
でかたがおかしい．
低地にしかいないやつの標高依存性がプラスだもんなぁ．
しかし random effects ということにして無理矢理補正してるよ．
-
このあと 23 時ごろまで試行錯誤の泥沼にはまることになる．
-
最初に，
「そもそも『いる・いない』のモデリングなんて，
もう不要じゃなかろうか?」
と考えてこれを削除してみた．
これは正解だったよーで，
パラメーターと deviance はすっきりすばやく収束するようになった．
-
それはいいんだけど，
列強な species
(各コドラートの上位 3-4 種)
たちが chain ごとに「談合」のごとき状態に入ってしまって
……
つまり，
おたがい相談してるかのように
「なるほどあんたが競争力 xxx という値をとるなら私は yyy
ぐらいの競争力があればいいんでしょう」
てなふうに xxx だの yyy だの値
(これはおもに random effects で決まっている)
が chain ごとにばらばらに「固定」されてしまって，
混交しなくなってしまった．
-
このあたり事前分布の変更で改善できないかしらんと工夫してみたり，
いろいろやってみたんだが
……
やればやるほど状況は悪化するばかり．
-
しかも夜になってアタマがへたばってくると，
凡ミスも増えて時間の浪費速度がすさまじく加速される．
-
とりあえず，
状況を完全に把握できるところまで撤退して本日は終了．
2300 研究室発．
2317 帰宅．
体重 68.8 kg．
ECOAS うんどう不足ふとり．
つくづく健康によろしくない．
晩飯．
-
[今日の運動]
-
うんどう休養日
……
そろそろ生活をたてなおさなくては
-
[今日の食卓]
- 朝 (0830):
米麦 0.6 合．
ニラ卵炒飯．
ダイコンの中華ふうサラダ．
- 昼 (1320):
研究室お茶部屋．
米麦 0.6 合．
インスタント味噌汁．
リンゴ．
- 晩 (2340):
明太子スパゲッティー．
キムチ & 豆腐．
2007 年 09 月 20 日 (木)
-
0750 起床．
朝飯．
コーヒー．
0920 自宅発．
曇．
0935 研究室着．
-
また ECOAS 八甲田湿原植物データ解析のつづき．
さすがに昨日あれだけいろいろ調べまわると，
何をどうやればいいのか，
少しは見当がつくようになってきた．
とはいえ，
やはり試行錯誤をつみかさねていくことでしか，
この統計モデルは改善されないわけだが．
-
計算まち時間に ECOAS 動物プランクトンに関する解析の結果をみなおす．
「脱」 ECOAS のためには，
-
得られた事後分布の表示が改善できないか検討する
-
この推定結果をもちいて，
何か ECOAS 的な予測というか憶測計算を示してみる
といったことが必要なのだろう．
-
八甲田湿原植物の試行錯誤．
競争力に上限・下限をもつような関数型にするとずいぶんとマシにはなるけど，
(昨日かいたような)
「談合で実現する種間競争力」
になってしまうところがあるねえ．
ここでは種間競争が完全に相対評価になってしまってるからなぁ
……
-
とはいえ，
絶対評価のがちんこ勝負にすると
「種間競争キビしすぎ」
な状況になる，
と (種多様性が高いほど現存量がどんどん減る
……
実際にそうなってる群落があっても不思議ではないけど)．
とすると，
バイオマス獲得の限界を種間競争だけではなく種内競争にも原因が，
と何らかの形で表現できればいいのか?
……
これは意味不明ぎみだな．
-
奇策は回避して，
相対的な種間競争に power を導入してみることに
……
これはよけーに事態を悪化させるとわかった．
-
いろいろ試行錯誤しつつ昼飯．
まずいなあ，
計算試行錯誤に時間をとられすぎてるよ
……
-
群落バイオマス分割，
ディリクレ分布
p[i, 1:N.sp] ~ ddirch(alpha[i,])
つかえばいいぢゃん，
と今さらながら気づいたんだが
……
unable to choose updater for node p[1,1]
なぜだー
……
dbern(...)
なんかとの組み合わせで無ければ動いてくれない，
ってコトかしらん?
共役事前分布でなきゃイヤイヤってのも横着なハナシだなぁ．
-
だったら自分で
ぢりくれ分布
つくればいいぢゃん，
と気づいたので，
alpha[i, j] ~ dgamma(mu.alpha[i, j], 1)
あれこれなどなどとしてみたんだが
cannot bracket slice for node alpha[12,11]
ゐんばぐす内の slice sampler がコケてるみたいだねえ
……
-
alpha[,]
に初期値 (1 とか) を与えたら計算するようにはなった
……
しかし 4000 MCMC step に 1165 秒も費やして
(slice sampler とかが苦闘してたんだろう)，
結果はやはり同じ，
と．
パラメーターとかの収束はいいのにねえ．
うーむ，
とはいえ，
やはり，
統計モデルがかなり良くない，
ってことかしらん?
-
で，
またまた今さらながら，
このコドラート内バイオマス「分割」モデルのよけーな計算に気づいた
……
「分割」だからって割算とかする必要ないぢゃん!
で，
憎むべき割算を取りさると計算は数倍速くなったんだが
(割算のせいで sampler によけーな心労・負担をかけていた，
ということだ)
……
うん，
列強談合体質はまだ完全には取り去れないね．
ちょー長時間計算すればいつかは一ヶ所にまとまるのかしらん?
-
どーしよーかなーとお茶部屋に逃げると，
そこにも湿原植物問題が．
といっても八甲田ではなくサロベツの湿原なんだけど．
江川さんから今シーズンの調査とその後について教えていただく
……
あいかわらずサロベツ的な状況
(A 棟 7-8F まわりでは「過酷・熾烈な重労働」を意味する)
のようで．
しかし昨年の春に播種した 3 種のうちヨシがめでたく全滅して，
残りはヌマガヤ・ミカヅキグサだけになったみたいだから，
解析したうけ業者としてはラクになったな．
-
1700 から研究院アワーに参加する．
気候力学分野の PD 奥西さんによる
「ニューラルネットワークを利用したマイワシの産卵回遊の再現」，
人工ニューラルネット (ANN) の脳をもつ
マイワシが北太平洋から土佐湾に帰ってくるように
人為選択 (genetic algorithm による)
かけてみました，
というもの．
わりとあらっぽい ANN でも行動を自動制御できたりするもんだな
……
-
1840 研究室発．
1855 帰宅．
晩飯の準備．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0820):
ヨーグルト．
リンゴ．
- 昼 (1330):
研究室お茶部屋．
食パン．
春雨スープ．
リンゴ．
- 晩 (2030):
米麦 0.7 合．
ジャガイモ・タマネギ・ニンジン・トマト・牛スネ肉のシチュー．
キャベツ・ニンジン・ネギのサラダ．
