ぎょーむ日誌 2006-10-(21-31)
2006 年 10 月 21 日 (土)
-
0715 起床．
朝飯．
コーヒー．
怠業．
洗濯．
メイルかき．
-
1030 自宅発．
晴．
1045 研究室着．
-
またメイルかき
……
ひたすらデータ解析こんさる，
というか私自身の統計学勉強というか．
-
また
昨年と同じく
12 月 9 日 (土) に統数研でアヤしげな雑談やることになりそう．
今回は
「ふつーのユーザーが R
まわりで MCMC 計算やろうとすると，
どういった事態におちいるのか」
といった内容になりそう．
-
とあるメイルでいただいた冗談:
日本で一番 R を使ってるのは久保
これは正しくは
日本で一番 R を使った下うけぎょーむやってる一人は久保
といったところだろう．
-
で，
その
R
を使ったサロベツデータ整理のつづき．
昨日の作図関数を汎用化して個体数の変化も図示できるようにする．
やっぱりカウントデータはいいねえ．
これに対して「カタチの変化」は悪夢だ．
-
1300 ごろ，
北大で開催されてる
日本科学哲学会第 39 回大会
で講演される
三中さん
がこちらに立ち寄ってくださった．
40 分ほど楽しく雑談．
当方が危惧してたとある企画は順調に suspend 中とのことで，
ひと安心．
-
来週からまた調査にでる江川さんと，
サロベツ統計モデリング相談．
なンか他にもいろいろデータあるわけね
……
そして階層的な構造 + 季節変化の個体数変化の推定か．
図上で数の変化をおってみると，
fixed effects は密度依存性 + サイズ依存性，
といったところ
(季節変化効果はサイズ変化と区別がつかないだろうからとりあえず無視)．
random effects はいつものごとく複数サイト内の plot 内の
……
うーむ，
「個体差」も必要とされてるのか?
-
いずれにせよ
Bayesian
だ．
ここ地環研 A 棟 7-8F では師弟ともども
「解析のこととか考えないで，
とりあえずデータとりましょう」
というハナシが多いんだけど
(で，
私自身も事前に「こういうふうにデータとるべきだ」
とか言えるほど現場も対象生物も研究目的も理解できてないし)
……
-
ともあれ，
まあ，
こういうふうにとられたデータは Gibbs sampler
で MCMC 計算シャカシャカやって，
パターンを説明できる統計モデルのとりうる「範囲」
みたいなのをゆるゆると見定めてやるほうがよろしいのでせう．
これに対して，
「ホントによい推定値を探しにゆきます」
最尤推定ってのは
ぎりぎりの曲芸のような数値計算ワザ駆使したりすることもあるわけで，
これはもっと random effects 少ない観測データでやるべきだよね
……
と言いつつ今までは無理矢理に最尤推定してたわけだが．
-
1850 研究室発．
1900 帰宅．
体重 72.6kg．
晩飯．
-
[今日の運動]
-
まあ，
明日はまた山行になりそう
(札幌岳)，
ということで
……
-
[今日の食卓]
- 朝 (0800):
食パン．
- 昼 (1230):
研究室お茶部屋．
食パン．
- 晩 (1930):
スパゲッティー．
豆腐トマトソース．
2006 年 10 月 22 日 (日)
-
0540 起床．
さて，
天気わるくなりそうだけど，
ホントに
札幌岳に行くのか
……?
0555 自宅発．
晴．
0610 研究室着．
朝飯．
ウスくいれたインスタントコーヒー．
-
0700 札幌岳登山隊出発．
甲山さん・川合さん・私．
また甲山さんのクルマで．
札幌市南区の温泉街・定山渓の横から豊平峡のほうに入る．
0815 札幌岳の登山口着．
[08:24 登山口出発]
標高約 450m．
曇．
札幌市外縁のヤマ，
ということで，
すでに朝から車が多数とまっている
(写真にはうつってないが)．
しかしいったん山に入るとなかなかヒトにあわない
……
本州・九州のヤマと比べて，
という意味だが．
[08:49 倒木 ver.2004]
最初の一時間ほどは冷水 (ひやみず) 沢ぞいを進む．
2004 年台風によるもの，
と思われる倒木があちこち．
植林されたトドマツとか．
まぬけなことに，
この小屋をみてようやくかつ突然に思い出したんだけど
……
私は一度札幌岳に来たことがある!
ときに
2002 年 11 月 3 日，
腰までうまる雪に行くてをはばまれ山頂目前にしつつも吹雪の中の劇的なる撤退
……
いやはやー，
ホントに忘れてました．
[10:04 急登の尾根，雪]
本日は北海道上空に強い寒気が入りこみ，
冷水小屋あたりから雪が降りだした．
そして小屋うらから
(2002 年次にはラッセルで苦闘した)
急登が 20 分ほどつづく．
甲山さんいわく，
「阿寒づかれ」．
[10:09 トラヴァース & 斜上]
急登がおわるとしばらく等高線ぞいにトラヴァース．
ここらへんはよく覚えてるなあ．
そして雪はいよいよ激しく．
次に，
急登尾根と主稜線のあいだのゆるい斜面を頂上めざすことになる．
[10:25 山頂直下ののぼり]
ゆるい斜面上の登山道は溝というか沢のようになってしまっている．
そしてここも倒木多数．
前回はこのあたりで撤退したはずだ．
[10:47 札幌岳山頂]
そして今回は雪のふるなか標高 1293m の山頂に到達することができた．
しかし，
二週間もすればこのあたりは腰までの積雪になっているんだろうな．
かなり寒いなかで，
昼飯．
ガスがかかっており展望わるい．
[11:47 下山中]
1105 寒い山頂発．
高度をさげるにつれ気温高くなることが実感できる．
降雪も少しよわまってきた．
[12:37 冷水沢渡渉]
前回
の記録にもあるように，
ここは何度も冷水沢の右岸左岸をいったりきたりするわけで．
[12:56 下山]
登山口駐車場のむこうにある標高 715m ぐらいの名もなき山こそが
「いかにも紅葉の山」ってかんぢでした
……
とはいえ，
この写真はいつものごとくヘボいな．
札幌岳は冬山になりつつある．
[13:00 問題の登山靴]
山道具屋では修理してくれそうもなかったので，
自分で強力接着剤つかって修繕した登山靴．
強度的には問題なかったと思う．
しかし接着剤を注入しすぎてシャンクが固くなってしまった．
ホントにびみょーな差なんだけど，
慣れるまで歩きづらかった．
登山口からちょっと下ったところにある温泉．
紅葉シーズンだったので満員だった．
露天風呂で甲山さんと雑談してると，
ついつい「風呂きりあげ時刻」をすぎてしまった
……
という失敗は
(これまたど忘れしてたんだけど)
2004 年 2 月 1 日
とまったく同じだ!
[15:17 温泉本格インド料理]
この温泉なぜか本格的な
北インド料理
がうりものだったりするんだよねえ
……
しかも大盛況なので整理券わたされて一時間ぐらい待った．
また甲山さんに御馳走になりました．
ありがとうございます．
温泉のカレー，
といってもイロものではなく，
ホントにおいしい．
-
1550 ごろ豊平峡温泉発．
札幌市への道は行楽がえりの車で渋滞気味．
1720 北大まで車で送っていただいた．
-
明日はまた新千歳空港まで来客の出迎え．
今回は
John A. Silander
さん．
時刻など確認したり，
また (ハナシのねたがないと困るので)，
Silander べいじあん論文をプリントアウトしてながめてみたり
……
Alan Gelfand
さんと生物分布の統計モデリング，
いろいろやっておられるわけで
(その一部は
Clark & Gelfand 本
にも掲載されている)．
1940 研究室発．
1955 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0625):
食パン．
- 昼 (1050):
札幌岳山頂でローソンにぎり飯 (梅)．
- 晩 (2030):
うどん．
2006 年 10 月 23 日 (月)
-
0710 起床．
0730 自宅発．
晴．
札幌市西側の
手稲山
(標高 1023m)
は昨日の寒気のせいでうっすらと冠雪してるな．
0745 研究室着．
朝飯．
コーヒー．
下は A 棟 8F 非常階段からのながめ．
-
空港まで
John A. Silander
さんを迎えにいくべく
0928 研究室発．
数秒の判断逡巡で
0940 発の快速エアポートをのがす．
発車ベルが鳴った瞬間にプラットホームの表示がきりかわるようで
(駆けこみ乗車防止のため?)，
「これホントに空港いきか?」
と迷ってしまったため．
そして目の前を去りゆく車両の中から手をふる三中さん
(本日つくばにお帰りのようで)
……
-
0955 JR 札幌駅発．
1032 新千歳空港．
下の写真ではよくわからんけど．
ANA の到着出口になぜか人垣
(それも修学旅行途中で空港にいるらしい女子高校生の)
ができていたので「何ごとかしらん?」
といぶかしく思ってたら，
(食べものを前にしたお茶部屋の院生たちのごとく)
興奮ぎみの高校生たちが「シンジョーが，シンジョーが」
と口ばしっておられるので，
日本ハムファイターズ関係か，
と憶測できた．
たしかに Silander さんの東京発便の前は中部国際空港発の便だな
……
[女子高校生の人垣→津波の自律的形成現象]
で，
プロ野球選手たちがこういう一般の到着口から出てくるはずもなく
……
やがていかなる機構なのかわからないんだけど，
「むこうのほうでバスにのってるらしい」という情報が伝達されたらしく，
津波のごとくどどっと移動していった．
彼女たちがどのように情報交換して人垣形成したり津波的移動したのか不明である
(けーたい電話?)．
-
急に閑散とした羽田発 ANA 57 便 (定刻より 10 分早く到着)
到着口で Silander さんと会えた
……
ネット上のフィールド写真みるとかなり長身のヒト
(190cm ぐらい?
あとで聞いたところによるとおじいさんは Sweden 出身とのこと)
だったので，
すぐに判別できた．
1119 新千歳空港発．
Silander さんの近ごろの関心は New England における
外来植物 (とくに日本産の) の侵入拡散過程なんだけど
(で，
生態学の研究には Bayesian 的な現象の理解のしかたが重要，
と言いだした一人で
Clark & Gelfand 本
の contributors の一人にもなっている)，
ハナシをうかがうと何でもやってきたヒトで，
もともとは生態遺伝学 (集団遺伝学) の出身だったり，
あるいは森林動態シミュレイター
SORTIE
の開発なんかも
(これは私も知っていた
……
というか私にとっては Silander さんといえば
「SORTIE なヒト」だったのだが)．
車中ではそのような雑談．
1158 JR 札幌駅着．
-
駅前のアスペンホテルにお連れして荷物を置いてもらう．
それから研究室にいっしょにいったんだけど，
甲山さんはまだ会議からもどっておらず．
Silander さんと地環研前の蕎麦屋で昼飯．
-
研究室にもどってあれこれとどたばた．
-
1430 ごろ夜おそくまでから矢澤さんと Pasoh データ整理・調査準備用
R
プログラム開発．
(清野さんふうに言うなら)
「毎木 1000 本ノック」ならぬ「樹高測定 1000 本ノック」なアレ．
本日のお題は熱帯林 50 ha 調査地の
38 万本毎木調査データが与えられたときに，
-
ある基準でもって
50 樹種合計 1000 本を選抜したものを
1 セットとし，
それを樹種重複なく 8 セット
(実際には 50 樹種 8 セット + 39 樹種はんぱセット)
構成せよ
-
つぎに選ばれた 1 セット 1000 本
(50ha 内のあちこちに点在)
の配置などをよく考慮して，
「調査しやすい」
データシートを自動生成せよ
というもの．
1. はともかく 2. は私の苦手な「ユーザーインターフェイス」にからむもので，
実際かなりたいへんであった．
-
先日，
「R には here document の構文がなくて」
とか文句いってたんだけど，
今回は「Perl 使うな，R だけ使用可」
という発注者からの厳命をうけてしまったので
……
貧すれば鈍する，
とでもいうのか下のような邪悪そうな書きかたをいろいろと
「発明」してしまいましたよ．
cat(sprintf("
\\textbf{\\gtfamily Set%02i(%02i/50): x = %04i-%04i}
\\hspace{10mm} \\textsf{Nov-Dec 2006} \\hspace{130mm} \\textsf{\\thepage/50}
\\hspace*{-20mm}
\\begin{tabular}{rrrrrr|%s}
\\hline
& tree\\# & spc & x & y & DBH & %s \\\\
\\hline",
i, page.number, x0, x0 + dx,
paste(rep("r", length(cols.input)), collapse = "|"),
paste(cols.input, collapse = "&")
))
現在、Ｂ棟２階にて有機溶剤の廃液（ブロロホルム）が漏れたことにより，
刺激臭が発生しています。
揮発による処理のため，館内に臭いが回っていますので，各室，ドアを閉
める，窓を開ける等により対応してください。
ただし，東側の窓を開けて処理していますので，東側の部屋の方は，窓は
開けない方が良いかと思われます（空気より軽いため）。
目の痛み，吐き気等の症状が発生します。
うがい，鼻の洗浄を行うようにして下さい。
軽い症状であれば時間経過により落ち着くようですが，必要に応じて医療
機関の診察を受けるようにしてください。
……
いやはや何がどうなっているのやら．
A 棟 7-8F 内のヒトたちの応答をみると，
この毒ガスにかなり敏感なヒトとそうでないヒトがいるとわかった
……
という実験をやりたかったのか，
コレは?
私はちょっと目がちくちくする程度ですんだんだけど．
そうか，
つまり「カナリア」
として使えるヒトがいる，
ということだ
……
ふーむ．
昼飯．
地衣類繁殖モデリング，
JAGS のデータファイルやコマンドファイルはめんどうなので，
これらを生成する R
プログラムを作る．
model ファイルは自動生成というわけにもいかず，
こういったかんぢになった．
model {
# inverse variances
iv.p ~ dgamma(1.0e-3, 1.0e-3)
iv.prock ~ dgamma(1.0e-3, 1.0e-3);
iv.pface ~ dgamma(1.0e-3, 1.0e-3);
iv.pid ~ dgamma(1.0e-3, 1.0e-3);
# factors
pspcU ~ dnorm(0, iv.p);
pspcL ~ dnorm(0, iv.p);
pareaU ~ dnorm(0, iv.p);
pareaL ~ dnorm(0, iv.p);
pslri1U ~ dnorm(0, iv.p);
pslri1L ~ dnorm(0, iv.p);
pslri2U ~ dnorm(0, iv.p);
pslri2L ~ dnorm(0, iv.p);
for (i in 1:n.rock) {
prock[i] ~ dnorm(0.0, iv.prock);
}
for (i in 1:n.face) {
pface[i] ~ dnorm(0.0, iv.pface);
}
for (i in 1:n.sample) {
pid[i] ~ dnorm(0.0, iv.pid);
apothecia[i] ~ dpois(exp(
prock[rock[i]]
+ pface[face[i]]
+ pspcU + (pspcL - pspcU) * spc[i]
+ pid[i]
+ (pareaU + (pareaL - pareaU) * spc[i]) * log.area[i]
+ (pslri1U + (pslri1L - pslri1U) * spc[i]) * slri1[i]
+ (pslri2U + (pslri2L - pslri2U) * spc[i]) * slri2[i]
));
}
}
計算に時間がかかるので
(といっても数千 MCMC step に 5-10 分ぐらい)，
二ヶ月ぶりに A801 室の Dell 機を起動して
こちらで計算．
上の model 定義はかなりじつはかなりあとになって整ってきたもので，
最初のうちはどうしたもんだか，
と試行錯誤がひたすら続いた．
たとえば「二種類の生物」というのはやっかいな状況である．
三種以上なら問答無用にて階層ベイズモデル化すればよろしい．
しかしながら，
二種の場合はそうすべきなのかびみょーなところだ
……
いろいろためしたあげくに，
上の model file にあるがごとく
Umbilicaria cylindrica
(タカネコゲノリ)
と
Lasallia pensylvanica
(オオイワブスマ)
を「同時かつ別々に」推定していくことに．
はてさて，
どうなることやら．
(後記:
よく見ると上の model
定義のこのあたりは冗長な書きかたになっているので
マネしないよーに)
ひたすら試行錯誤．
どのようなモデルをもちいようとも，
三段階のややこしー階層性をもつ
このパターンを説明するパラメーターの事後分布たちは
数千 MCMC ステップぐらいでは「収束」しない，
とわかってきた．
本日は院生密度が低い．
静かだ．
院生密度が低いと露崎さんの出現頻度が高まるらしい，
という観察が得られた．
ばててきたので計算命じておいて撤退．
1930 研究室発．
JR 札幌駅の本屋によってから帰宅．
1955 帰宅．
[駅で野球観戦するヒトたち]
地元球団だまって応援するたちの異様な迫力がにぢんでいます．
-
そして「栗御飯用のクリ」を研究室におきわすれてきたことに気づいたので，
(速歩) 往復 20 分間を費して取りに行った．
このクリはとある複雑なる事情によって，
昨日「塩づけ」みたいな状態にされてしまい，
さらに「塩ぬき」と称して過去 24 時間水没させられていたもの
……
かかる状況が発生した件に関して私は責任とらされる立場ではないんだけど，
「まあ，
この水びたし塩クリは栗御飯ぐらいにしか使えないだろう」
と当家でひきとることになった次第で
……
お茶部屋ではおりあしく院生の皆さんの晩飯直後で，
私の行動・服装などに関していろいろとキビしく指摘されてしまったので，
すばやく逃げる．
帰宅して晩飯の準備．
晩飯．
-
そして A 棟 8F 闇ネットにもぐりこんで計算状況をしらべてみる
……
-
[今日の運動]
-
うんどうできるかできないか，
これは一日の「どこに」うんどう時間をもってくるか，
にかなり依存してるんだよね
……
また昼ジョギング復活か?
-
[今日の食卓]
- 朝 (0840):
研究室お茶部屋．
食パン．
- 昼 (1320):
研究室お茶部屋．
食パン．
- 晩 (2200):
「写真とれ」
と命じられたので撮っておきました
……
[栗御飯]
塩ぬき塩クリを使っていつもの炊飯用土鍋
(炊飯器が壊れたままなので)
で炊いたもの．
まあ，
このクリはそれなりに良いもので栗御飯もそこそこにおいしかったんだけど，
いかんせんずっと水づけだったんでかなり崩壊しててですね
……
なんというか
「飯にちょっとマッシュポテトを混ぜた状態」
みたいになってしまった．
2006 年 10 月 26 日 (木)
-
0710 起床．
朝飯．
コーヒー．
0840 自宅発．
曇．
0855 研究室着．
-
昨日の JAGS コードはあれこれ冗長なのでまた書きなおす．
やっぱり「二種類の生物」ってのは難しいねえ．
二種といってもサンプル数同じぐらいの場合と，
かなりサンプル数が異なる場合では考えかたを変える必要あるだろうし．
今回は二種が 2:1 ぐらいのサンプル数で，
なかなかナヤましい．
-
地衣類繁殖べいずモデル MCMC 計算のつづき．
このモデリング，
数千ステップでは MCMC 計算が「収束」しない，
ということが確認できた．
burn-in に 10-20K MCMC step が必要で，
sampling も数万 step の規模になる
……
パラメーターによっては
Gibbs sampler が数千 step 周期ぐらいで
「うねる」みたいなんで．
まあ，
モデルに何か問題あるかもしれん，
という可能性も検討してますよ．
-
ふーむ，
(よけーそうなやつも含めた)
全パラメーター使うとやはり「収束」
がよろしくない
……
この程度 (ってどんな程度だ?)
の問題でもパラレルテンパリングとかいった
複数チェイン計算わざが必要とされているのか?
たぶん DIC とかも「収束」してなかったら計算しても意味なさそうだし．
-
まあそのうち，
JAGS
以外の Gibbs sampler，
OpenBUGS だの WinBUGS だのでの挙動も調べてみる必要あるな．
-
で，
説明変数を増やすと下のごとくよたつく，
と．
もっと長い burn-in と sampling range が必要，
ということか．
ついでに記録しとくと，
25K MCMC step でおよそ 25 分の計算だった．
-
ヘンなハナシかもしれんが
……
link 関数 log のもとで，
variance parameter の
超事前分布をホントに「無情報事前分布」にしてしまうと
すごく効率が悪くなるような気もする．
無情報だけど variance に関する超事前分布のばらつきやや小さめ，
といった設定にすると「うねり」が少ないような
……
当方のモデルがヘン，
あるいは JAGS の Gibbs sampling アルゴリズムに
何か問題あるのかもしれんけど．
いやいや，
やはり sampling の長さをもっと長くしないことには
……
-
少なくとも知りたいいくつかのパラメーターの事後分布は，
この超事前分布の「無情報ぐあい」に左右されているところがある
……
左右されるパラメーターのほうが少数派か?
-
ともかく事前に burn-in の長さとかは決めようがない．
出力された結果をみて，
ということになる．
-
会計係より連絡．
昨日のどくガス事件，
の後始末のよーで．
１０月２５日に発生した標記の件について，各研究室内で医療機関を受
診した教職員学生がおりましたら，氏名を会計係までお知らせ下さい。
お忙しいところ，大変恐縮ですが，現時点での判明分だけで結構です
ので，取り急ぎ願います。
また，追加で情報がありましたら都度ご連絡願います。
なお，医療機関を受診した場合には，領収書等を保管するようご指示願
います。
昼飯．
午後もひたすら計算つづく．
burn-in 20K step + sampling 40K step で 40 分間，
か．
パラメーター 2 個増やすと
burn-in 20K step + sampling 40K step で 47 分間．
で，
この長 sampling & パラメーター増計算で，
(正しいかどうかはともかくとして)
「地衣類は暗すぎるところ & 明るすぎるところは好きじゃない」
なる結果が
(個体 / 岩面 / 岩の nest した「個体差」 random effects を考慮してなお)
でてきたのでこのあたりをさらに長 sampling で調べてみることにする．
予想計算時間 80 分間
(後記: 実際には 74 分間だった)．
このように断続する計算まち時間に
OpenBUGS
ペイジを少し整備してみたり
……
ゐんどーづユーザーは
WinBUGS
を使えばよいわけで，
OpenBUGS や JAGS を使うのは
-
ゐんどーづを使わない & 使いたくないヒト
-
WinBUGS を使えない & 使いたくないヒト
に限定される，
というぢつに狭い範囲にしかウケない情報提供になるわけで．
VMware 上の R に
R2WinBUGS
を install.packages()
してみて example(R2WinBUGS)
走らせてさらにその
## Not run:
コード
(help(R2WinBUGS))
から先を動かしてみたら，
じつにじつに簡単に WinBUGS が動作して，
しかも R から自由自在にコントロールできてしまったのでくやしい．
JAGS や OpenBUGS における苦労はいったい何なのか，
と．
Linux 上の VMware (仮想マシン) 上で Windows XP 使うよりはまだしも
敗北感のウスい
Wine
(Linux 上で動作するゐんどーづソフトウェア実行環境;
ゐきぺでぃあ解説)
を導入してみるか
……
で，
まずこの wine のコンパイルがすごくたいへんなモノで，
ThinkPad X31 だと一時間以上を費した
(あとで Dell desktop 機でやってみたけど，
それでも 25 分かかった)．
ちなみに Linux OS 中核部であるカーネルは
5 分もかからずにコンパイルできるんだけど
……
まあ，
wine のコンパイルの様子とか見てると，
graphics まわりに莫大な手間をかけてるみたいだったけどね．
とりあえずの試験運転とて
winecfg
を動かしてみると
(ここで $HOME/.wine/
以下のディレクトリが作られるとわかった)，
こういういかにもゐんどーづ的なイヤらしい GUI の
configurator
が起動した．
fontforge
による無理やりぎみな日本語表示はすごい．
-
さらに試験運転ということで，
(Linux 上で使う予定はないんだけど)
R むけテキストエディター
Tinn-R
をインストールしてみる．
Tinn-R installer をダウンロードしてきて，
wine Tinn-R(ヴァージョン番号).exe
でインストール開始．
その後は
wine ~/.wine/drive_c/Program\ Files/Tinn-R/bin/Tinn-R.exe
などと起動すればよいのか?
下の図は Linux デスクトップ上で，
Tinn-R の installer が動いていたり，
Tinn-R で Linux 上で書いた R コードを開いてみたところ．
-
さてさて，
ゐんばぐすだ．
インストールは簡単でダウンロードしてきたファイルを
wine WinBUGS14.exe
とするだけ．
(後記: 1.4→1.4.1 の patch あてが必要，
翌日のぎょーむ日誌で説明)
-
そして問題はこいつを「R から使う」というところだ．
install.packages("R2WinBUGS")
して，
Rでベイズ統計学
などを参考にしつつ，
ぢたばたしてみる．
-
とりあえず動かせるところまでできた．
library(R2WinBUGS)
して example(bugs)
するとモデル定義ファイルやデータがよみこまれ初期設定がなされる．
そこから先の WinBUGS 呼び出しは
## Not run:
で抑止されている．
ゐんどーづ上なら上記のごとく簡単にその次に書いてある bugs()
できるわけだが，
Linux & Wine ではここで工夫が必要であり，
私のところでは下のようにすると
R から Wine をつうじて WinBUGS に計算を命じることができ，
さらに計算結果を受けとることができた．
WINEPATH <- "/usr/local/bin/winepath" # ここで定義する必要あり
schools.sim <- bugs(
data, inits, parameters, model.file,
n.chains = 3,
n.iter = 5000,
#bugs.directory = "c:/Program Files/WinBUGS14/", # wine 上では不可
bugs.directory = paste( # そこでこのように変更
Sys.getenv("HOME"),
".wine/drive_c/Program\ Files/WinBUGS14/",
sep = "/"
),
working.directory = NULL,
clearWD = TRUE,
program = "WinBugs",
useWINE = TRUE,
newWINE = TRUE,
WINE = "/usr/local/bin/wine", # /usr/local/ にインストールしたので
WINEPATH = WINEPATH
)
で，
うけとった計算結果の作図だが
……
plot(schools.sim)
で自動的に呼び出される bugs オブジェクトの
default plot() (bugs.plot())
は下のごとくいまいち感があるので
……
-
いつものごとき図を描かせるには下のごとくやればよい．
library(coda)
mcmc.ss <- as.mcmc(schools.sim$sims.matrix)
par(mfrow = c(6, 4), mar = c(3, 3, 1, 1), cex = 0.5)
plot(mcmc.ss, auto.layout = FALSE)
-
いやー，
これであっさり WinBUGS が使えるようになってしまった
……
ゐんどーづ使いばかりがそろっている院生たちから
「私の PC でも MCMC 計算つかえるようにしなさいよ」
と命じられたらどうしようかと危惧してたんだよね
……
JAGS
や
OpenBUGS
はゐんどーづ版もあるんだけど，
どちらも R からは直接よびだせないからね
(さらに厳密に言えば，
OpenBUGS を呼びだせないのは Linux 環境でのハナシで，
ゐんどーづでは
library(BRugs)
を使えばよびだせる)．
-
1835 研究室発．
1850 帰宅．
体重 72.6kg．
晩飯の準備．
洗濯．
晩飯．
-
2205 自宅発北大構内走．
寒い．
しかし寒い夜空の下を走っているヒトたちが他にもいた．
2255 帰宅．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0740):
米麦 0.6 合．
栗御飯．
ダイコン・ネギ・煮干の味噌汁．
- 昼 (1340):
研究室お茶部屋．
食パン．
- 晩 (1930):
米麦 0.7 合．
栗御飯．
チンゲンサイ・ネギ・豆腐のシチュー．
2006 年 10 月 27 日 (金)
……
と，
どんどん修正していったあげくにつきつけられた最後通牒とでもいうべきエラー:
educational version cannot do this model
てめー，
このばぐす野郎どもめ
……
と立腹してみせてもしょうがないので，
Registration Form
(登録は無料なんだけど)
をかいて送信．
すぐに educational version 解除 key を暗号化して送ってきた
……
のは良いんだけど，
This key expires on 31st December 2006
とは?
なんでこんなにすぐに無効になるわけ?
もう WinBUGS なヒトたちの思惑は考えないことにして作業をススめる．
上の patch あてと同じ手順で制限解除．
制限解除してもまだ WinBUGS の Gibbs sampler
はあーだこーだと文句ばかりつけてきて計算はじまらん
……
trap window とかが出て，
shape parameter (r) of gamma iv.pface too small -- cannot sample
……
ということで超事前分布 (hyperprior) の
「無情報ぶり」
を調節してみたり
……
といった悪戦苦闘のすえに，
よーやくこの地衣類繁殖モデルに関する試験運転に成功．
モデル定義ファイルはまだ JAGS との互換性を維持している．
すでに正午．
気のせいか「さっさと昼飯くっとけ」
といったようなご指示をうけたような気がしたので，
北大生協で昼飯調達．
その前に 100K MCMC step の試験運転を R2WinBUGS に走らせておく．
お茶部屋でさっさと食ってしまう．
昼飯後もまだ時間的ゆーよを与えられてるようなので，
試験運転的な計算のつづきを．
なんと JAGS だと 74 分かかる 100K MCMC step計算が 42 分で終了．
しかも混交のぐあいも一見したところ良さげに見える
(下の図にはこのあたりは示されていないが)．
収束も速いのかもしれん
(とすれば，
こんなに長い計算は不要ということになる)
……
そして得られた bugs object をそのまま R で
plot()
するとこうなる．
ふーむ
……
-
で，
次に超事前分布をエラーがでる限界ぎりぎりまでひろげてみると
……
(また 40 分間ほど経過)
……
結果がびみょーに変わった．
アヤしいパラメーターはやはりアヤしい，
という結果になったな．
-
収束のぐあいを調べるため，
n.burnin や n.iter
をかなり短く，
そしてチェイン数を 3 にしてまた bugs()．
50 分ぐらいかかるだろう．
-
ぱそー
下うけぎょーむ．
10 月ぶんは本日で終了する予定．
-
野外調査用データシート自動出力の改善
-
野外調査用
5m×5m
局所地図 1000 箇所「地図帖」自動生成
-
50ha 地図
-
R
+
LaTeX 生成するデータシートは以前と同じで列名を変更したり．
-
「地図帖」は
5m×5m
subgrid ごとに 1000 樹木ぶん，
こんなかんぢで 84 ペイジ出力．
これは
R
だけで作れる
(そして
pdf() で出力)．
-
そして 50ha 地図はふつーにかんなかんぢで．
点々は選抜された (第一セットの) 1000 個体．
-
ということで夜までかかって Pasoh 仕事ひとまずは終了．
-
志水さんの大雪山系地衣類ベイズモデルの WinBUGS
Monte Carlo トリプルチェイン計算のほうは
……
ふーむ，
興味ぶかい．
-
2015 研究室発．
2030 帰宅．
体重 72.4kg．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0800):
研究室お茶部屋．
Sunks のサンドイッチ．
- 昼 (1240):
研究室お茶部屋．
北大生協の巻きずし．
389kcal．
- 晩 (2100):
うどん．
チンゲンサイ・ネギ・豆腐のシチュー．
2006 年 10 月 28 日 (土)
-
0740 起床．
朝飯．
コーヒー．
洗濯．
ひたすら怠業．
昼飯．
1330 自宅発．
晴．
1355 研究室着．
-
また大雪山系地衣類繁殖努力のベイズモデル，
R
→
R2WinBUGS
→
Wine
→
WinBUGS
な MCMC 計算のつづき．
WinBUGS を使おうが何をしようが，
この計算には時間がかかることが明らかになりつつある
……
-
R メイリングリストに流れた情報によると，
Thomas Yee さんの R package
VGAM
が
CRAN 上で公開
されたよーで
……
ただし，
CRAN ミラーサイトに行きわたるにはまだ数日を必要とするようだ．
-
地衣類繁殖べいづモデルの計算を試行錯誤しつつ，
R2WinBUGS
めもペイジの整備，
とか．
-
「明るさ」 (0 - 1) みたいなものは log link 関数の線形モデルの中では，
log(明るさ) でいれたほうが解釈しやすいな
……
こうしておけば完全暗黒のもとでは (カウントデータである)
繁殖努力の期待値がゼロになるんで．
-
で，
このあたりをいろいろとひねくっては MCMC 計算
……
という悪しきループにつかまる．
-
まち時間にいろいろと志水さんデータながめてみる．
事後分布にみょーな傾向があったので，
ひとつの岩面に複数種いる事例はどれぐらいあるんだろう，
と調べてみた
……
R
の
tapply()
は便利なわざでこういう集計にむいてる．
これはうっかりまちがいなんだけど，
> tapply(d$species, d$face, function(x) length(x))
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
6 5 5 5 4 7 11 9 8 8 6 7 8 7 7 6 7 5 8 6 4 10 6 9 8 8 11 7 11 6
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
9 11 9 7 11 19 11 12 7 10 6 14 15 15 12 11 10 10 11 12 16 6 11 14 9 19 8 13 13 17
61 62 63
25 13 17
63 面ある岩面でそれぞれ何個体の地衣類を観察したか，
という集計になっている．
各面で何度もサンプリング (4-25 個体) してることがわかる
……
ということで岩面の random effects を考慮しないといけない，
ということで，
すでに「岩 - 岩面 - 個体」と三重に nest した random effect factor
はベイズモデルに組みこまれている．
-
では，
各岩面に「何種の地衣類がいるか?」
を正しく集計させると，
> tapply(d$species, d$face, function(x) length(unique(x)))
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
61 62 63
1 1 1
……
となり，
私にとっては意外な結果だったんだけど，
じつは志水さんの調査では
「ひとつの岩面に複数種がいる」
という場合がひとつもない，
と今さらながらわかった．
みょーな事後分布
(モデルによっては岩面の random effects
と種の効果がみょーに相関してるようにみえる)
はこれで説明されるわけだ．
-
それでは「岩面」 (63 面) ではなく「岩」 (24 個) でみた場合，
「ひとつの岩にいる地衣類の種数は?」
をみると，
> tapply(d$species, d$rock, function(x) length(unique(x)))
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 2 2 1 2 2 1
24 個の岩の中で複数種いる岩は 6 個だけ，
ということか．
ふーむ．
-
ついでに 63 面の岩面が 24 個の岩にどのように配置されてるか，
は
> tapply(d$face, d$rock, function(x) length(unique(sort(x))))
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
2 2 2 2 2 3 2 2 2 2 3 6 4 4 2 2 3 2 3 3 2 2 4 2
こうなる．
-
ベイズモデルの計算はうまくいかぬまま撤退．
2325 研究室発．
2340 帰宅．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0810):
バゲット．
チンゲンサイ・ネギ・豆腐のシチュー．
- 昼 (1230):
うどん．
- 晩 (2145):
研究室お茶部屋．
豚角煮弁当．
500 円．
2006 年 10 月 29 日 (日)
-
0730 起床．
朝飯．
コーヒー．
怠業．
-
昨日から，
というかここしばらくひっかかっているのは，
10/20
にやった「明るさ」なる計算なんだよね
(下の図はそのときに志水さんデータをつかって
R
で作図したもの)．
これは夏至の日の南中時刻のひあたりを相対値化したものである．
[大雪山主稜線上の地衣類]
……
が好む岩面?
志水さんデータを使った作図．
北緯
43°
において日当たりが悪いところはダメそう．
しかしながら，
日当たり良すぎるところもやはり
……
-
まあ，
これで太陽をくるくるまわしたところで，
あまり上の図は変わらないだろう，
とは思うわけだが
……
-
で，
この「明るさ」と呪われ割り算値で表現されてる繁殖努力の関係はこう，
と．
[呪われ割り算値]
ディスプレイ上で作図してすぐに消去する場合は
「うん?
ユルい相関?」
とかいいかげんな感想もつばかりだけど，
こうやってじっくり見ていると，
次第にいろいろ気になってくるわけで
……
-
やっぱりこういうカウントデータを割り算値化するのは
ぢつにバカげた所業ですね，
ということで志水さんの元データを使って図を作りなおす．
やはりこちらのほうがずっと安心感がある．
-
上の図でわかることは，
私がてきとーに三次元計算してみせた
「明るさ」とやらは繁殖努力とあまり関係なく，
やはり個体サイズそのもののほうがよっぽど支配的である，
ということだ．
そしてこの overdispersion は
random effects の大きさを表現している．
-
apothecia 数 200 以下の部分だけを図にすると，
こうなる．
そうすると今度はナリはでかくてもあまり繁殖してないやつがいるわけで
……
これは何か「明るさ」とやらで説明できるのか?
それとも「人間の知ることができる限界」こと
random effects の嵐の中にかき消されてしまうのか?
-
1020 自宅発北大構内走．
使うかどうかはともかく，
走りつつアタマの中で球面座標系を設定して，
そこで太陽を動かしてみる．
「夏至の日の日の出から日没まで何時間か?」
といった問題を考えるときは，
「天動説」的に不動の大地から見あげた
天球上に太陽の軌道を設定してやると考えやすい
……
「地球は自転してるから」
とか考えるとわけわからなくなりがちである．
そしてその世界の「外側」の
「東/西」にたってながめ，
次に「北」からながめれば，
どうやって昼間時間を計算すればよいかわかる，
と．
-
1120 帰宅．
体重 72.2kg．
昼飯．
-
1315 自宅発．
晴．
1330 研究室着．
-
昨晩帰るまえに R2WinBUGS
に命じておいた，
n.burnin = 30000, n.iter = 60000
のトリプル MCMC 計算はとっくにオワってるわけで
……
計算時間はおよそ 5000 秒だったか．
きれいに収束している．
そして fixed な「明るさ」 (今回は logit 変換)
のパラメーターの事後分布に関しては，
まあゼロにひっかかる「妥当な」結果ともうしますか．
-
上の軌道とこの下の図をみくらべると，
意外と夏至の日ってのは北側からも
……
-
本日は (娯楽的な) 太陽軌道の計算だけ．
「明るさ」計算の再検討はまた明日以降，
とゆーことで．
2000 研究室発．
2015 帰宅．
晩飯．
-
輪読会の予習．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0800):
バゲット．
- 昼 (1200):
うどん．
ニラ卵．
- 晩 (2110):
米麦 0.7 合．
チンゲンサイ・ネギ・ショウガ・エノキダケ・豆腐のスープ．
2006 年 10 月 30 日 (月)
-
0740 起床．
朝飯．
コーヒー．
0840 自宅発．
晴．
0855 研究室着．
-
帯広の平林さんが用事のついでにこちらによってくださったので，
R
による作図こんさる．
-
よく見られる R 作図「よくないわざ」に関する注意を一点．
ひとつの
plot()
内に複数の図を描きたいときに
par(new = TRUE)
は使うべきではない方法．
こういうときは，
plot()
したら，
そのあとの線・点の追加には
lines()
だの
points()
だのを使う．
-
1030-1200 研究室の
輪読会，
Silvertown & Charlesworth の ``Plant Population Biology (4th ed.)''．
第一章も Fig. 1.5 (c) とかすごくまちがった図の引用あるわけで
(正しくは Caswell (2000) の行列モデル本 p.61)．
どうなることやら．
-
また昨日の太陽軌道計算プログラムとかひねくる．
-
昼飯．
-
で，
「太陽を動かした場合の岩面の明るさ」
なるアヤしげな数値を R で計算させたわけだが
……
左が以前のもの
(夏至の日の南中時刻だけ)，
右が今回のやつ
(夏至の日の日の出から日の入まで，
0-1 で規格化!!)．
-
なンというか
「岩面なんてどっちむいてたって別にどうでもいいぢゃん
(特に極端な日かげをのぞけば)」
という方面にハナシがすすみつつあるような．
まあ，
それはそれで妥当な気もするが．
とりあえずこの数値でもって R2WinBUGS
な MCMC 計算を命じてみる．
時刻は 1600．
-
ついでに新しく計算しなおした「明るさ」なる数字と，
地衣類の子器数をながめてみる．
まあ，
以前よりはマシな気もするのだが
……
-
で，
繁殖のベイズモデルの MCMC 計算の結果はこうなりました，
と．
以前と変わらずというところか．
-
1840 研究室発．
1905 帰宅．
体重 73.0kg．
晩飯の準備．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0800):
米麦 0.6 合．
チンゲンサイ・ネギ・ショウガ・エノキダケ・豆腐のスープ．
- 昼 (1340):
研究室お茶部屋．
北大生協の巻きずし．
389kcal．
- 晩 (2055):
米麦 0.7 合．
ニラ卵．
ブナシメジ・ワカメ・煮干の味噌汁．
2006 年 10 月 31 日 (火)
-
0750 起床．
朝飯．
コーヒー．
0910 自宅発．
曇．
地衣類繁殖モデルに組みこまれている nest した random effects
について検討しつつ，
そのへんをふらふらと歩く．
0935 研究室着．
-
「個体 in 岩面 in 岩」
という nested random effect 構造の中で，
もしも
『よけーな』
random effect factor があるとすれば，
それは何か?」
……
と考えてみたんだけど，
「これはなくてもいいでしょ」
という部分はなさそう．
-
たとえば，
「個体」の差ということを考えるとすると
……
たとえば一枚の岩面上に複数の個体がいるように見えたとしても，
それってもともと一個体だったんだけど，
いつのまにかちぎれてしまって，
と考えると一個体の複数回観測だよな．
ということで，
試しに「個体」による random effects ぬきで MCMC 計算やってみる．
-
計算まち時間．
とくに意味もなく，
R
の
pairs()．
-
あるいは今日のセミナー予習の論文よみとか．
-
「個体の差」はすごくすごく重要そうじゃない，
と 2500 秒の MCMC 計算によって示された．
これは意外というか
……
つまり一枚の岩面上であっても場所による差だとか，
あるいは地衣類 (あるいはその中のソウ類?)
たちの「性能差 (遺伝子差にもとづく?)」
がかなり重要みたいということなのだろう．
-
1030 より 1.5 時間ほど，
研究室セミナー
で，
塩寺さんの熱帯山地林葉寿命ハナシ．
多種・多変量・季節変化・pseudo replication・形態，
と生態学データ解析をひどく難しくしてしまうあれこれがてんこもり
……
なんだけど，
もし (random effects ぬきの)
比例ハザードモデルが葉寿命なる平均値をうまく推定できている，
とするならば，
-
窒素濃度
(mg g-1)
なる割り算値より
LMA
(g cm-2)
なる割り算値のほうが
現象論的なモデルの説明変数としては適切そう
(このデータセットには)
-
「生活形」内では LMA による従来モデルで
葉寿命は (目的論的に) 説明できていそう
(「コストの高い」葉ほどモトをとるのに時間がかかる)
-
しかしながら「生活形」間で比較すると
来モデルでは説明できない
(いろいろな生活形の中で，
高コストそうな葉っぱをつけてる連中は
ぱかぱかと使い捨ててるように見え，
一方で安そうな葉をもつ連中は
それをけちけちと使いのばしている?)
ぐらいにまとまるのではなかろーか?
3. の「説明できない」はその指摘だけで終わってよいような．
「新理論」を構築できるほどデータもないし
(というか，
これ以上データがあっても難しすぎて解析できそうにないし
……
人類文明の科学の限界?)．
-
地衣類計算，
「個体の差」ほどではないが「岩面の差」も重要みたいだね，
とわかった．
もはや「総あたり」的に，じゃあ「岩の差」はどうなの?
とまた MCMC 計算を命じておいて，
昼飯．
-
母子里林冠とりまとめな作業とか．
-
「岩 random effects」ぬきトリプルちぇいん計算結果でた
(計算時間 3000 秒)．
「個体 in 岩面 in 岩」
という nested random effect 構造の中では，
こいつの重要性がいちばん低そう
……
朝から「重要」だのそうぢゃないだのといった
アヤしげでうさんくさいことばかり述べているが，
これは R から呼び出されている (wine 上の) WinBUGS が吐きだす
DIC
(deviance information criterion)
なる数値をながめてつぶやいているわけで
……
DIC が良いかどうかはともかく，
やはり階層ベイズモデルの世界にもモデル選択規準があったら便利ではあるよな．
単なる平均 deviance でいいぢゃん，
とういことなのかもしれんけど．
-
母子里作文を校閲にだしたり，
Annika さんの email アカウントを申請したり．
-
また計算結果が．
こんどは「明るさ」のいれかたを変えてみた
……
ふーむ，
よけーな変換するなってコトかしらん?
「岩」効果を復活させ，
link 関数は log なので
(ナマの「明るさ」だと意味不明なのでやむをえず)
「せんたりんぐ」した「明るさ」で再計算を命じる
(しかしばかばかしいのでこのあとは log とか logit 変換した値にした)．
ホントにいきあたりばったりだな．
-
昼飯が食パン一枚だと日没のころにはかなり空腹感がせまってくる．
しかしながら
……
この空腹感なるものは
視床下部
の摂食調節機能とかに由来するものとのことで，
じつはいまだにそのあたりのしくみがわかってないらしい．
たとえば，
血糖値は空腹時でも低下するわけではないらしく，
というのも低下したら脳とか動かなくなるわけで．
ともかく，
いますぐ何か食わなければ餓死するわけでもないのに，
「はらへった」
警報を発するのはなんとかならぬものか．
自分の制御プログラムとかも
テキストエディターか何かで簡単に書き換えられればいいんだけど．
-
そして院生の計算機に R インストール作業
(こういう単純作業をやると空腹感は忘れられる)
……
たしかに CRAN の
「どこに最新版 R があるの?」
ってのはわかりにくいねえ．
-
まずは
http://www.r-project.org/
の News で最新の version 番号 (2.4.0 とか)
を確認する
-
CRAN のミラー，
たとえば
(もしあなたがゐんどーづ使いなら)
ftp://ftp.ecc.u-tokyo.ac.jp/CRAN/
の下の
ftp://ftp.ecc.u-tokyo.ac.jp/CRAN/bin/windows/base/
にある最新版 R
(ここでの例だと
R-2.4.0-win32.exe)
をダウンロード，
それを「実行」してインストール
(後記:
いやはや，
当然ながら
ftp:// 経由でつなぐより
http:// でつなぐほうがわかりやすいわけで
……
R のインストール
参照)
-
1920 研究室発．
雨．
1935 帰宅．
体重 72.8kg．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0830):
米麦 0.6 合．
ブナシメジ・ワカメ・煮干の味噌汁．
- 昼 (1340):
研究室お茶部屋．
食パン．
- 晩 (2010):
うどん．
ブナシメジ・ワカメ・煮干の味噌汁．
