ぎょーむ日誌 2006-11-(21-30)
2006 年 11 月 21 日 (火)
-
0820 起床．
うーむ．
朝飯．
コーヒー．
0920 自宅発．
曇．
ちょっと雨．
0935 研究室着．
-
何やらいろいろとメイルかきとか
……
げ，
もう 1100．
A 棟 7-8F だけでなく，
日本のあちこちの皆さんがデータ解析にとりくむ季節になりましたなぁ．
-
昨晩の小山さん解析の WinBUGS ばぐとりに取りくむ．
しゃれにもならん
……
挙動が小山さん PC と私の Linux + wine ではびみょーに異なっており，
Trap 窓とかは出ない．
何なんだろう．
そしてわりとあっさり原因を特定できた．
最強に強まったばかさ加減をこめた必殺の一撃，
というべきか「右カッコがとじてない」という
ちょー脱力系のばぐだった．
コンパイルできないなら「コンパイルできない」と言えよ
……
というか WinBUGS はバカなので無理やりコンパイルしたあげくに
実行して自爆しているのか?
-
そして簡単なモデルなのに，
データが多い
(6000 観測地点，
ただしほとんど場所は植物個体数ゼロ)
という理由で計算に時間かかってるみたいだ．
えーい，
もっとステップ数とか少なくすべきだったか?
-
10 分間まっても終わる気配がないので中止．
n.iter
を 5000 から 500 に変更．
これまた時間かかってる
……
これでも時間かかるか．
ならば 50 なら?
-
私の ThinkPad X31 (PentiumM 1.6GHz)
で 50 MCMC step トリプルチェインで 30 秒弱，
か．
500 でやってみたら 250 秒ぐらい．
モデルが簡単なので，
これぐらいで十分に収束しているけれど．
-
モデルをちょっと改良するためパラメーター数をひとつ増やしたら，
計算時間が 400 秒に増えた．
くそう，
WinBUGS め
……
収束は速いからもっとサンプル数へらしていいかな?
Chain 数 3 だし．
いや，
やっぱりこれより少なくはできんような気もする．
-
まあ，
なんとかなってきたのでちょっと長めの試験運転めいじておいて，
昼飯．
-
昼飯後に小山さんにひきわたし作業
……
またゐんどーづ上で WinBUGS を動かすのに苦闘する．
そして WinBUGS が動作してるあいだは他の作業ができない．
なんともすごい．
-
A801 の共用ゐんどーづ機
(Pentium4 2.8GHz なのでこの研究室内では最も速い計算のひとつ)
で計算環境ととのえる
……
これがまた腐れぎみでいろいろ苦闘する．
なんで私がこんなところでじたばたしてるんだろうね．
-
いろいろ設定したけど R2WinBUGS
から WinBUGS を動かせない
……
と思ったらその原因はひどくゐいんどーづ脱力系のものであった．
R2WinBUGS の
bugs()
はいろいろなファイルへのパスを指定せねばならんのだけど，
その中に Shift-JIS の文字コード
(具体的には「...\デスクトップ\...」)
が混入していると WinBUGS は「ファイルがみつからん」
叫んで仕事を投げだしてしまう，
とわかった．
へっぽこ OS & へっぽこアプリケイションソフトウェア．
ディレクトリ移動で対処．
-
(私は「種ごとの推定」は好きではないんだけど)
とりあえず全種・全 stage の「いる・いない」推定を走らせてみる
……
が二番目の種でコケる．
原因をひとことで言えば，
午前中に私が試行錯誤して作ってたプログラムというか，
その考えかたが抜本的にまちがっていたため．
データがない区画では「ゼロ」データで区画をうめつくすのではなく，
そもそも「そんな区画はない」として扱えばよい．
そうすれば「その区画の『区画差』」の事後分布は
事前分布からの単なるランダムサンプルとなるだけで，
計算全体には何の支障もない．
嗚呼，
われながらへっぽこ統計もでらー，
というかんぢだ．
修正して，
また最初から再実行．
-
とっくに札幌の日没時刻をすぎてる．
いやはや．
-
お，
小山さんところから転送したデータの中に今回つかってる
「やち (ぼうず)
パワー地図」
pdf ファイルがまじっていた．
-
A801 共用 PC 上で全種・全 stage で (とりあえずひととーり)
お試し計算を実行させようとする試行錯誤，
小山さんととりくむ．
一昨日あたりにつくったデータ整理用 R プログラムに
まだ不備が残存していたのを発見・修正したり．
-
「サロベツ調査地の植物データ，いる・いない問題」
の階層ベイズモデル化
(といっても現段階では GLMM でまだ代替可能だが)
WinBUGS 実装は次第に完成しつつある
……
そして次の一手が気になりつつある．
Zero-inflated Poisson (ZIP) モデルであつかうような問題，
BUGS 言語ではどう表現してやればよいのか?
-
そしてたまたまなんだけど，
鳥データ解析の WinBUGS あれこれなどについて
山口さん
からメイルをいただき，
上記 ZIP モデルとの関連が何かありそうだと気づく．
どちらも行き詰まっているところは
-
WinBUGS などで使われている
BUGS 言語はどんな統計モデルでも
自由自在に記述できるわけではない
-
しかしながら，
私の憶測では，
ちょっとややこしげな確率分布をもちいた尤度方程式をもつ
統計モデル (ZIP など) であっても，
何らかの方法
で
「BUGS 言語にあわせた」
等価なモデルとして書き直せるような気がする
-
ただし統計モデルの中には
(尤度方程式は記述できても)
BUGS 言語で記述できないものがあるかもしれない
(あってもぜんぜん不思議ではない)
といったところか．
-
2000 研究室発．
午後はずっと雨ふってたけど，
今はやんでる．
統計モデリングについて考えつつ夜の札幌駅まわりをふらふらと歩く．
ふらふらと
……
なーんだなーんだ簡単なコトじゃないか．
BUGS 言語でもって
隠れマルコフモデル
(もしくは単なる
潜在離散変数モデル)
を実装すればいいだけじゃないか．
それで zero-inflated model も
「個体識別マークつき鳥を見た・見なかった」
時系列観察データもモデル化できる．
ホントに Gibbs sampler が走るかどうかはわからないけど，
定式化は簡単だ．
隠れ変数を悪用して
dbern() や dpois()
が「『実質的に必ず』ゼロをだす」ようしむけてやればいい．
まあ，
明日以降に試してみよう
……
とこれにてめでたく一件落着の結末に到達したので帰還軌道に遷移する．
2050 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0840):
バゲット．
海藻スープ．
- 昼 (1340):
研究室お茶部屋．
食パン．
- 晩 (2220):
米麦 0.7 合．
ダイコン・ネギ・煮干・ワカメの味噌汁．
2006 年 11 月 22 日 (水)
-
0640 起床．
よく眠れんかった．
朝飯．
コーヒー．
0805 自宅発．
雨．
0820 研究室着．
-
農学部は電力まわりでやられてるみたいだけど，
それに関係あるのかないのか，
ここ A 棟 8F でも一回だけ瞬間停電があった．
サーヴァーたちは無停電装置で保護された．
-
データ解析こんさるメイル書き．
-
昨晩おもいついた隠れ変数モデルを BUGS 言語で書いてみて，
ちょっと実行してみる
……
ふーむ，
MCMC 計算の最中に
abundance[i] ~ dpois(0.0)
みたいな状態におちいってもエラーとかは出さないみたいだ．
マトモな計算ができているかはいささか疑問だけど．
-
計算時間は今までの二倍以上かかりそうだな
……
とりあえず 200 MCMC step で 300 秒か．
二倍にはなっていないか．
そして最初の 200 step ぐらいは計算できる，
とわかった．
-
じゃあ，
もっと長くしてみたら破綻したりするかもね，
と 600 MCMC step ぐらいにしてみる．
私の ThinkPad X31 では非力な気がするので
scp
で A801 室の Dell Linux デスクトップ機におくりこみ
ssh で実行を命じる．
-
600 step で 330 秒．
やっぱり速いね．
そして，
これも予想どーりというべきか収束がよろしくない．
というか，
すごく悪い．
興味ぶかい．
-
隠れ変数の悪用方法をちょっと変えてみる．
具体的には
「平均ゼロのポアソン分布」
→
「かぎりなく平均ゼロにちかいポアソン分布」．
いや，
あまり結果は変わらんような気がするな．
よくわからんまままた実行
……
この変更によて計算時間は 370 秒に増大．
そして収束は
……
うーん，
ちょっとマシぐらいのかんぢだろうか?
まだまだぜんぜんよくないけど．
-
こういうのってもっと長く計算すれば解決するものなのかしらん?
ということで
n.iter = 3000
としてみる．
予想計算時間は 1850 秒，
31 分間ぐらいか．
-
またメイルかき．
-
1780 秒．
収束は
……
三重連鎖がかなりヒト山になってきたけど，
まだまだだな．
そして「やちパワー」の影響がかなりあいまいになっている．
ところが「やちパワー」地図でみると，
いま試験運転用にみている植物の当年生実生はかなり
「やち」重視のようなので，
これはやはり私の隠れ変数モデリングに何かへっぽこなところが
まだまだ残存しているのだろう．
基本的なかんちがいとかもあるかも．
-
計算はいったん休止してメイルかきのつづき．
全国的にデータ解析のしーづんなんですかねえ
……
-
さーて，
また統計モデルの修正でも
……
とそちらに没頭してるとまた運動不足になるので，
1210 研究室発，
北大構内走．
晴．
構内は落葉．
と思っていたらいきなり雪がふってきた．
1255 研究室もどる．
昼飯．
-
ZIP を実現する隠れ変数階層ベイズモデルのみなおし
……
うーむ，
一ヶ所 variance がいれちがいになってるところがあった．
しかしこれぐらいで計算がおかしくなるものかな?
とりあえず 300 MCMC step で試運転．
370 秒で終了．
混交は今までよりかなりマシになった．
これでよさそう．
-
さて，
となると「やちパワー」のききかたの問題だけど
……
これって「やちぼうず」直上で最大になるんだけど，
そこには他の植物は入れないんだよね．
やちぼうずを形成しているスゲとかが密集しているから．
そのあたりうまくモデリングせんといかんような気がする．
-
いきなりすごい雪．
窓の外はホワイトアウト．
-
安直ながら
「やちパワー」
の二次関数にしてみる．
2 パラメーター増加．
計算に時間かかってるな
……
520 秒かかってぜんぜんダメ，
とわかった．
考えの足りない安直な複雑化は無意味である，
と．
-
いきづまったので，
小山さんと昨日の計算結果である事後分布の
自動作図ぷろぐらみんぐ．
-
そしてまた階層ベイズ版 ZIP モデル
……
うーむ，
よくわからん．
-
お茶部屋雑談．
本日はめずらしく「食」のハナシではなく，
衣裳に関すること
……
「しゅうかつ」
におけるコートの適切な長さの重要性について．
-
やちパワー階層ベイズ ZIP モデル，
-
Logistic 部: ブロック間
-
Poisson 部: ブロック内小区画間
にばらつき (random effects による)
があると仮定してみるのは?
それほど悪くない仮定だと思うんだが
……
計算時間は 434 秒に増大 (600 MCMC step)．
ぜんぜん収束してない．
-
じゃあ Poisson 部分で random effects がない (!)
としてみるのは?
……
おっとその前にまた定式化の
「ちょっとした」
(などとごまかすにはやや深刻な)
まちがいみつけた．
いやはや．
でそれを修正して 300 秒計算させると，
よーやくまともっぽい答えがではじめた．
収束まだまだなんだけど，
DIC = 1152
(なんとなくこの値が気になってしまう)．
-
もいっぺん Poisson 部にブロック間差を組み込んでみる．
335 秒．
これまたぜんぜん収束していない
(上もあわせてもうちょい長くとれば良さそう)．
DIC = 1243.7
-
じゃあもとにもどって小区画ごとに Poisson 部の random effects
があるとするのは?
これはやはり収束が劣悪．
DIC = 1215 か．
-
ブロック間差あり，
1000 回．
550 秒．
収束ぜんぜんダメ．
DIC = 1295
……
ブロック間差なし．
500 秒．
これまた収束ぜんぜんダメ．
DIC = 1419
……
-
とはいえ，
いろいろな計算結果に共通してあらわれる傾向，
というのはだんだんわかってきたな．
「やちパワー」の効力とその限界，
みたいなのが．
-
といった現状報告 +
まともな計算にはまだまだ時間かかるやもしれません，
と小山さんに報告．
まあ許容される進捗状況のようで．
-
yuji_nk さん
から Bayesian ZIP モデル情報おしえていただく．
ありがとうございます．
たくさん Bayesian 本かいている
Perter Congdon
の一冊
``Bayesian Models for Categorical Data''
というのがあって，
その 5.4 節 (p.167) あたりに Bayesian ZIP の解説あり，
とのこと．
対応する WinBUGS コードは
ftp サイト内
の
BMCD_Programs.zip
の Program 5.4.odc．
-
私も粕谷さん科研費でその本を買っていたので (しかし読んでない!)，
その zero-inflated Poisson (ZIP)
& zero-inflated Negative Binomial (ZINB)
を読んでみる
……
うーむ，
あまり Bayesian っぽくないふつうの ZI モデル解説のような気がする．
生態学まわりではあまり使われてない「ゼロぬき」
ポアソン分布を使ったモデルの説明とかもあり．
-
Congdon の WinBUGS コードの解読を試みると，
たしかにこれは ZIP の BUGS 的実装になってるよな
……
隠れ変数
T[i]
が 1 か 2 の値をとるようにして
(わざわざ dcat(...) で二項分布からのサンプル，か)，
ポアソン分布の平均をふたつ用意して
(そのうちひとつは mu[i, 1] <- 0.0 )，
y[i] ~ dpois(mu[i, T[i]])
とかやってるわけね．
-
私の方式は「隠れ変数の悪用」と称してるだけあって，
同じことを実現するのにもう少し邪悪に
(あるいはふまじめに)
やっている．
ブロック差・やちパワーなどに依存する
二項分布にしたがう隠れ状態
T[i]
は 0 か 1 をとる値で (T[i] ~ dbern(...))，
これを用いてポアソン分布の平均を
mu[i] <- T[i] * exp(線形予測子)
と表現している．
これでもちゃんと動作するし，
自分ではこちらのほうが簡潔だと思うんだけど
(もしかしたら Congdon BMCD 本が書かれたころの
version の WinBUGS ではこの方式の計算できなかったのかもしれない)．
-
森林総研の永光さんから Cc でメイルいただいたので，
7F 院生室に (放置していた) シウリザクラ作業の確認を
……
で，
わかったことは
先月 10 日
の「三者面談」以降，
字句どおりの意味において何ひとつとして進捗していなかった，
ということ．
さすがは野郎大学院生，
とやはり自滅型野郎大学院生のなれの果てである
私などは感銘をうけざるをえないのであった．
-
「けっきょくボクはこういうことをやったらいいんでしょうか?」
発言などまだ長き冬眠状態から覚醒しきっていない症状がみられたので，
観測データは取り上げてこちらで解析することに．
自分でも意味がある行為なのかどうかよくわからないんだけど，
図の修正などはやっておくように指示してみる．
-
「やちパワー」試行錯誤をつづけながら，
計算まち時間に
R
の中でそのシウリザクラ実験データを整理整頓．
まあ，
なんとかなりそう．
-
2100 研究室発．
2115 帰宅．
いつのまにか積雪．
晩飯．
[地環研前]
昼ごろには雪などまったくなかったんだけど，
いつのまにやらこういう状態に．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0700):
米麦 0.6 合．
ダイコン・ネギ・煮干・ワカメの味噌汁．
- 昼 (1310):
研究室お茶部屋．
食パン．
- 晩 (2230):
米麦 0.7 合．
豆腐．
ダイコン・ネギ・煮干・ワカメの味噌汁．
2006 年 11 月 23 日 (木)
-
0800 起床．
朝飯．
コーヒー．
あ，
今日は祝日か，
ということで洗濯．
怠業してる場合ではないんだが
……
と思いつつも洗濯機がまわると怠業もーどになってしまう私．
-
怠業ついでに運動だ．
1140 自宅発北大構内走．
曇天．
路上のあちこちに雪や水たまりがのこってて走りにくい．
そしていつもの往復区間である
低温研-獣医学部前の歩道は完全に雪に被覆されている．
1225 帰宅．
昼飯．
-
1310 自宅発．
曇．
1325 研究室着．
-
さーて，
昨晩かえりぎわに
「あとはもうどうにでもなれ，
とりあえず全種・全 stage で計算しろ」
と命じておいた MCMC 計算だが
……
小山さんから事後分布作図の
R
プログラムを拝領して図にしてみたところ，
やっぱり収束してないとわかった．
モデリングのいいかげんな状況での
なげやり計算はダメですな．
だめモデルの挙動の把握には役たつけど．
[事後分布の作図例]
まだぜんぜん信用できない計算結果です．
-
たしかにうまくいってない計算も多い
……
一種ずつ取りだして図にしてみると，
しかしここしばらく MCMC 計算の実験対象に使っている
植物コード
Solidago s (アキノキリンソウ実生)
はどのパラメーターもきれいに収束してやがるな．
標本数が多いせいか?
上の図でもミカヅキグサとかはひと山に見える．
ふーむ
……
-
お茶部屋会話．
「久保さん，
そこにおみやげのお菓子が」
「ありがとうございます．
しかしお茶部屋のお菓子を取るのは勇気がいるので，
あとで誰も見てないときにこっそりといただきます」
-
サロベツ ZIP 階層ベイズモデル，
無謀なモデル改造の一環としてブロック (line) ごとの
random effects を多変量正規分布
dmnorm(平均, 逆数分散共分散行列)
からださせるようにする．
今までは完全に独立した正規分布から
生成させるようにしていたんだけど，
収束の悪さはどうもこのあたりに何か原因があるような気がしてならない．
-
さて，
めんどうなのはこの新しい事前分布 (prior) たる多変量正規分布
の超事前分布 (hyper prior) だ．
単なる正規分布の場合だと，
その分散 (の逆数) の事前分布はガンマ分布を「なんとなく」
というか，
共役事前分布 (conjugate prior) だからという理由で使ってきた．
では，
多変量正規分布の分散共分散行列 (の逆数) の共役事前分布は何か?
それは
Wishart 分布
なんだよね
……
ようこそ，
アタマがヘンになりそうな多変量確率分布の世界へ．
-
R の
library(bayesm)
に入ってる rwishart()
や library(MCMCpack)
に入ってる rwish()，
すなわち Wishart 乱数はきだす関数でしばらく遊んでココロをおちつける．
で，
コードの改造にとりくんだわけだが
……
無情報 Wishart 分布にするためのパラメーターの与えかたがわからん．
調べてみてもナゾ．
-
えーと今回の多変量正規乱数の次元は 2
……
ということで，
なンか対角行列になっていればよい正定値行列 R を
R
の中で
R <- structure(c(1, 0, 0, 1), .Dim = c(2, 2))
として data
として R2WinBUGS
経由でわたす．
Wishart 分布の自由度は 2．
われながらでたらめだな
……
-
600 MCMC step 計算，
472 秒で終了．
収束はわるくない．
DIC = 1383．
問題の多変量正規乱数由来の事後分布もそこそこ値がゆれてる．
-
よくわからんので，
R <- structure(c(0.01, 0, 0, 0.01), .Dim = c(2, 2))
としてみる．
482 秒．
「逆数」分散共分散行列の値はむちゃくちゃでかい，
ということは分散とかすごく小さいとなるわけだが
……
そのあおりをうけて，
なのか DIC = 1293．
-
まだよくわからんので
R <- structure(c(100, 0, 0, 100), .Dim = c(2, 2))．
とうぜんこれは上とは逆で
「逆数」分散共分散行列の値はむちゃくちゃ小さい．
なぜか
deviance 小さく DIC = 1284．
すごくチェイン間差がある，
つまりぜんぜん収束してない．
-
なんとなく
R <- structure(c(0.1, 0, 0, 0.1), .Dim = c(2, 2))
というかもっと簡潔に
R <- diag(0.1, 2)
でよいのかしらん，
と勝手に決めて 1000 MCMC step 計算やらせてみる．
n.burnin = 700．
この問題に関しては 300 step のサンプリングでは足りないのかな?
-
786 秒．
DIC = 1242．
aP
とか収束あまりよくない．
-
R <- diag(0.01, 2)
にしてみるか
……
810 秒，
DIC = 1095．
収束はまだまだ．
-
お，
伊東さん苦闘の
「Mac OS X上で、WinBUGS+R2WinBUGSを使用する」
が
(とそのまわりの
経緯も)
公開されている
……
CrossOver と Darwine のコンビネイションわざ，
でしたか．
うーむ，
すごい．
おつかれさまです．
-
サンプル長かえてみる．
n.iter = 1300, n.burnin = 700, n.thin = 2．
……
計算時間 1030 秒，
DIC = 1294．
たいへんよく混交している．
-
多変量正規分布の共分散の事後分布の平均値は負だった
(かなり「ゆるい」けど)
……
で，
その解釈はどうなるわけ?
「植物の定着できる場所は限定されてる & いるところはたくさん」
VS 「どこでも定着できる & しかしどこでも低密度」
の区別がつかん，
ということかしらん．
そしてなぜ Gibbs sampler
は無相関の正規分布ペアから
「負の相関をもつ事後分布」
を生成できないんだろうか
……
いや，
生成してたのかな?
-
よくわからんところもあるけど，
ということで，
今晩の全種・全 stage 長時間計算は
n.iter = 1600, n.burnin = 1000, n.thin = 2，
昨日とのちがいは多変量正規分布の事前分布
& Wishart 分布の超事前分布ということで．
1900 研究室発．
謎のお供をする晩飯．
2100 帰宅．
-
データ解析こんさるメイル書きつづく．
R
FAQ．
「わくの外」に何かを描くには
par(xpd = NA)
が必要．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0820):
米麦 0.7 合．
海藻スープ．
- 昼 (1250):
タラコスパゲッティー．
- 晩 (1930):
北 14 西 1 のグロビュール．
スープカレー．
2006 年 11 月 24 日 (金)
-
0830 起床．
0840 自宅発．
曇．
ちょっと雪．
0855 研究室着．
-
昨晩めいじておいた計算，
Moliniopsis-f
(ヌマガヤ開花個体)
で停止しているな
(Trap 窓が出現)
……
集計してみるとたしかにここは 5
個体しかいなくて推定できない，
と．
> sapply(c("s", "j", "f"), function(st)
+ sapply(v.target.spc, function(sp)
+ sum(data[data$spc.stage == paste(sp, st, sep = "-"),]$Freq)))
s j f
Rhynchospora 232 367 287
Hypochaeris 299 180 39
Drosera 153 47 172
Solidago 172 165 13
Moliniopsis 169 132 5
Lobelia 12 153 8
Eriophorum 62 60 1
Carex 17 92 3
ということで，さっそく小山さんと相談．
また多変量正規乱数の活躍の場が増えてしまいそうな予感
……
とりあえず，
お茶部屋で朝飯．
コーヒー．
コードみなおしてるうちにセミナー時間．
1030
研究室セミナー，
本日は北村君の阿寒アカエゾマツ直径成長速度解析
(年輪データ)．
ガンマ分布の GLM つかった解析．
まあこんなもんか，
というかんじで
……
とりあえず同じ個体から何度も観測値とっている
縦断的 (longitudinal) なデータ構造なので，
random effects
考慮した混合モデル化が必要，
といった毎度のごとき指摘とか．
これでかなりすっきりするだろう．
このテの問題となると甲山さんコメントが
(よくも悪くも)
かなり「独特」なものになり，
つきあいの長い人間でなければなかなか理解できぬ甲山わーるどが展開する．
で，
甲山さん言うことをもし仮に全部マにうけたとすると，
すごく複雑な階層ベイズ + random field モデルになるな
……
とアタマの中でモデリングしながら拝聴してたんだけど，
ご本人はそのへんをどこまで理解しているのやら，
興味ぶかいところだ．
[いつのまにか雪]
かとー記念おふぃすからの光景
(東側)．
さてセミナー終了後は北大構内走によい時間なんだが
……
けっこう雪がふっている状況．
うんどう不足の日々をすごすうちに冬になってしまったのか．
計算はうまくいってないけど，
だからといって「自分を罰するために雪の中を走る」
という心境にはなれないしなあ．
-
成長段階 (
s, j, f の 3 段階)
をばらばらに解析するのではなく，
これらをまとめて，
かつ同時に stage 間の相違も推定するよう BUGS 計算まわりを改善する．
6 次元の多変量正規分布を事前分布としてつかわねばならぬ，
と．
-
Wishart 分布まわりのばぐにやられつつ，
なんとか改造終了．
お，
雪がやんでるな．
ということで，
北大構内走にでることに．
ついでにその間は長めの計算やらせてみる．
-
1345 研究室発北大構内走．
雪が路上にあるとやはり走りにくい．
そして突発的に雪がふったりやんだり．
1420 研究室もどる．
昼飯．
しかしお茶部屋内の院生密度が急激に高まる event が発生したので，
すばやく逃げ出す．
-
6 次元多変量正規分布モデル，
計算にえらく時間かかってるな．
40-50 分ぐらいで終了するかと思ったけど，
すでにその二倍以上が経過している．
まあ，
途中で止めるのもなンだし，
最後まで計算させてみるか
……
う，
きちんと確認できなかったんだけど 1000 step に二時間ついやしたような．
あと 600 step．
-
お茶部屋で宮田さんと
輪読会
の予習．
第 10 章の ESS 種子サイズだの移動分散の進化だの書いたの
Deborah Charlesworth か?
うーむ，
計算結果だけ示されてもねえ
……
-
アキノキリンソウ 3 stage まとめて ZIP 階層ベイズモデル，
3 時間 20 分を費して終了
……
しかし BUGS 言語で書いた model 定義に
一文字
の書きまちがいがあり，
遠路 12000 秒の MCMC 探索の旅路の果てに得られた何もかもは
一瞬にして鳥有に帰したのであった．
嗚呼．
気をとりなおして再計算を命じる．
-
メイリングリスト
stats
にながれた
R user 会 (12/8, 12/9)
の
abstract．
なんと Ripley 先生は 12/8 だけでなく 12/7 も講演するのか．
残念ながらこれには参加できないが．
それより自分の準備がですね
……
いやいやまだスケジュール破綻ではないけれど．
-
600 MCMC step ではまだまだ収束してませんなあ．
これでも 4440 秒ついやしたんだが．
計算そのものは破綻なく進捗するようなので，
8 植物種ぜんぶ，
合計にして十数時間かかりそうな計算を命じる．
-
計算まち時間は先日「とりあげた」シウリザクラ受粉実験データの解析．
「花序差」
を考慮した GLMM
による解析．
定番の
glmmML()，
cbind(結実した, しなかった)
表記法は使えるは stepAIC()
は使えるは (今日は使わなかったけど)，
ホントに便利になりましたねえ．
-
2004 年の自由集会
で紹介した
グループわけ
関数を使っていろいろな組みあわせ
(4 処理あるので 15 とーり)
のもとで
glmmML()
してみる
……
という部分はささっとできたんだけど，
見ばえのよいテイブルを作るところで苦闘というか，
かなり「汚い」 coding になってしまった．
set AIC autogamy geitonogamy control outcross SD
(naive prob.) (690) 0.08 0.17 0.33 0.51 --
1 (a)(c)(g)(o) 207.9 0.06 0.15 0.31 0.57 0.87
2 (a+g)(c)(o) 212.5 0.10 0.10 0.31 0.57 0.93
3 (a)(c+g)(o) 214.0 0.06 0.22 0.22 0.57 0.95
4 (a)(c+o)(g) 215.8 0.06 0.14 0.44 0.44 0.96
5 (a+g)(c+o) 219.6 0.10 0.10 0.44 0.44 1.01
6 (a+c+g)(o) 229.1 0.15 0.15 0.15 0.58 1.12
7 (a+c)(g)(o) 231.1 0.15 0.14 0.15 0.58 1.12
8 (a)(c+g+o) 233.9 0.06 0.32 0.32 0.32 1.17
9 (a)(c)(g+o) 235.9 0.06 0.33 0.31 0.33 1.16
10 (a+c)(g+o) 251.0 0.15 0.33 0.15 0.33 1.36
...
やはり 4 処理ごとに結実率が異なるとするグループわけ
(a)(c)(g)(o)
が最小になったか．
(naive prob.)
ってのは (結実合計)/(花数合計) という単純なわりざんち．
この AIC が劣悪である理由は
「花序間の overdispersion がとても大きい」
ということだ．
繁殖生態学のヒトたちって，
なぜかしら割算だい好き &
このあたりきちんと解析する気がないんだよねえ
……
よくあるのは花序ごとの割算値 (すなわち平均値)
のそのまた平均値の計算，
とか．
それって下の図のごとく花序サイズ
(○ の大きさ)
がけっこうちがうデータでやったら，
ホント意味不明な平均値ですよ
……
ということでこういうありがちな受粉実験データ解析なんかも
glmmML()
とかでやったほうが良いと思うんだが．
-
てなことを (苫小牧ではなくこちらの)
平尾君に (無理矢理) 見せて説明してるうちに，
上のややこしいテイブル生成の中にバグあること発見
(上に示しているのはすでに修正ずみ)
……
見せびらかして良かった．
ついでに
「上の図みたいに非アルファベット順 box plot横軸を生成するには?」
なるいつもの R FAQ うけたので，
バグとり御礼ということで，
factor(data$treatment, levels = c("a", "g", "c", "o"))
「レヴェル変換わざ」
の実演つき説明．
-
というあたりで本日は撤退．
2020 研究室発．
雪やんでる．
2040 帰宅．
体重 71.6kg!
ひさびさに 71kg 台の世界にかえってきましたなあ
……
今年の最悪の時期
(6 月ごろ?)
とか今より +3kg ぐらい太っていた．
晩飯の準備．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0940):
研究室お茶部屋．
食パン．
- 昼 (1440):
研究室お茶部屋．
食パン．
- 晩 (2210):
米麦 0.7 合．
ハクサイ・ネギ・ショウガ・豆腐・ワカメのシチュー．
2006 年 11 月 25 日 (土)
-
0720 起床．
朝飯．
コーヒー．
0910 自宅発．
晴．
0925 研究室着．
-
昨晩命じておいた計算は 15 時間ちかく経過したのにまだ終了せず．
えーと，
予想終了時間は
……
あと 30 分ぐらいか?
1 植物種に 1 時間 50 分 - 2 時間
(6400-6800 秒)
ぐらいを費しているんで．
-
この計算，
こんなに時間かかるようになったのは事前分布の設定がまずかったせいだろうね，
と今朝になってよーやく気づいた．
たとえば 6 個のパラメーターを多変量正規分布から生成させるときに，
(なんとなく R 的発想で，というか)
6 個まとめてえいやっと計算するのはこの場合はじつはかえってまずい．
というのも，
Wishart 分布を超事前分布とする
6 次元の分散共分散行列には 36 個の (自由に決められるのは 21 個か?)
パラメーターの Gibbs sampling が必要になるわけで
……
-
おそらくこのサロベツ ZIP 階層ベイズモデルの場合，
この 6 個のパラメーターは 2 個
×
3 回という Gibbs sampling で問題ないはず
……
もしそうであるなら，
これに必要とされる
2 次元の分散共分散を決めるのに必要なパラメーター数は 3 で，
7 倍のちがいがあるよな．
-
で，
このサロベツ計算に関する 2 個 ×
3 パラメーター化を
BUGS 言語で書くためにはおそらく 3 次元配列を使わんといかん，
と
……
このあたりの使えるデータ構造に制約があったりするところが
BUGS のヤなところだ．
-
うーむ，
最後の植物種
Carex
の計算が今朝の 0815 に始まって 110 分経過したんだけどまだ終了しない．
あ，
ヌマガヤ
(Moliniopsis)
とかは 8000 秒，
つまり 130 分以上かけて計算しているな．
Gibbs sampling の難しい植物種とそうでないのがいるわけ?
-
昨晩から始まった計算，
8000 秒弱を費して無事に終了
……
と思ったら，
ぜんぜん収束とかよくないし．
モデル定義とか書きなおして，
再計算せんといかんね．
-
再計算命じておいて，
ひさしぶりに映画みにいくことに．
1105 研究室発．
大学からもよりの (かなり小さな) 映画館
蠍座．
1130 から「うなぎ」，
今年なくなった今村昌平監督
(そういえば原作「闇にひらめく」の吉村昭も今年なくなっているな)，
役所広司・清水美砂主演
……
リアルそうに見せながら，
じつはわりと空想的・幻想的な展開のハナシ．
役者はうまいヒトばかりだと思うんですが，
それよりも舞台となり撮影も現場でなされた利根川ぞいにある佐原市
が何か不思議な水郷集落のように思えてしまうところがいいですねえ．
-
昼飯．
1520 研究室もどる．
計算のほうは
……
やれやれ，
初期値の与えかたにおけるいいかげんな改造のせいで，
二つめの植物種の計算途中でエラー吐いてとまっていた．
-
そしてめんどう計算担当 Dell 機の Wine が腐っている
……
いや，
またしても
$HOME/.wine/ ディレクトリがやられていて，
つまり「原因なにもわかりません」ということで
……
おそらくゐんどーづに独特な
($HOME/.wine 内におかれている)
れぢすとりファイルとやらが
wine 自身に壊されているんだろうね．
そんなものは手動で修正できんので，
$HOME/.wine
抹消してから，
手もとの ThinkPad から rsync.sh
でりもーと複写．
これでなおった．
-
R2WinBUGS
な R
コードをまた修正して再試行
……
1000 MCMC step を 1360 秒．
だいたい 6-7 倍の高速化か?
ただし混交はへぼへぼ．
もっと長く burn-in & sampling せんといかんわけ?
1000 MCMC step，
予想計算時間 2200 秒 (37 分間)，
か．
-
あららら，
上の計算試験はブタナ (
Hypochaeris
……
ってこんなタンポポみたいなやつだったのか)
で今回は Solidago だからサンプル数がけっこう異なる．
たぶん高速化も 2-3 倍程度だろう
……
-
1600 MCMC step に 84 分かかった，
ということで高速化は 2.2 倍程度．
これは 3 stage まとめての計算で，
1 stage ごとにみるとおよそ 28 分間
……
まあ，
こんなものなのかな．
しかし混交はぜんぜんダメ．
-
多変量正規分布の分散の事後分布がけっこう小さいので，
超事前分布にさかのぼってそれを少し「ゆるめる」よう設定して，
全種計算を命じておいて撤退．
1840 研究室発．
1855 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0800):
米麦 0.5 合．
ハクサイ・ネギ・ショウガ・豆腐・ワカメのシチュー．
- 昼 (1420):
我が家の近くにある定食屋
「卯和」
の昼定食 (注文は 14:30 まで) 750 円．
ベトナムふうカレー焼鳥が主菜のあとはふつうの和定食．
おいしかった．
- 晩 (2030):
米麦 0.7 合．
ハクサイ・ネギ・ショウガ・豆腐・ワカメのシチュー．
2006 年 11 月 26 日 (日)
-
0730 起床．
朝飯．
コーヒー．
怠業
……
この部屋，
ダンボール箱 (おもに Amazon 由来)
で足のふみ場もなくなってきたんで片づけなくては，
とのろのろと「室内移動用の通路開削」にとりくんでみたり．
何しろちらかりすぎてるんで，
ガスファンヒーターの温風経路が完全に閉ざされてる
(ので使用できない)．
お，
一ヶ月前に行方不明なってた腕時計が発掘された．
あ．
ボールペンも．
-
洗濯．
片づけつづく．
「倉庫番」
というパズルゲイムありますよね．
まさにああいう状況．
えーい，
われながら狭い部屋をせまくつかっているな．
あれこれひっくりがえしていくうちに，
またいろいろなものが発掘される．
そして室内の大気汚染は際限なく悪化していくので，
窓をあけて寒風をふきこませて息をつく．
-
とりあえず暖房の温風回廊は掘削できたので，
本日の作業は終了とする．
ばてた．
肺の中の空気を交換するべく
1055 自宅発北大構内走．
とうぜん部屋の窓はあけはなしたまま．
晴．
1150 帰宅．
体重 71.8kg．
-
昼飯の準備．
流しまわりに捨てるべきゴミが堆積していてとても不便だ．
昼飯．
1340 自宅発．
晴．
1355 研究室着．
-
うーむ，
ぼこぼこピーク，
すなわち，
ぜんぜん収束してない．
計算時間は 7-8 時間ぐらい．
-
やちパワー地図と見くらべると，
それなりに妥当な計算結果
(
Rhynchospora のやち嫌い疑惑とか)
ではあるんだが
……
こうも事後分布がでこぼこではねえ．
-
こういう問題で本来なら重要なはずの空間相関なんかも無視して
(これがよくないのか?)
単純化したモデルでなぜかくも計算が困難なのかよくわからない．
計算に莫大な時間が必要とされる理由も
……
-
とりあえずネット雑用．
-
室温 10°C
ぐらいか．
やはり (とくに全館スチームヒータ切られる休日は)
寒いので今しーづん初のかとーオフィスのガスストーヴ点火．
-
策がないので，
サロベツ ZIP 階層ベイズモデルの
「やちパワーの影響の stage ごとのちがい」
まわりを変更．
これまで多変量正規分布に生成させてたんだけど，
無情報な正規分布を事前分布にしてみることに．
とりあえずまた
Solidago
で試験運転．
願わくばさっさと計算終了してほしい．
-
そして R2WinBUGS
から WinBUGS を wine (非ゐんどーづ OS における
ゐんどーづアプリケイションソフトウェア実行環境)
経由で呼び出しているんだけど
……
どうもこの wine は
徹夜計算にヨワい?
疑惑が浮上してきた．
128% 意味・原因不明だが．
また wine 自身が
$HOME/.wine/
以下のれぢすとりファイルをぶち壊してくれたみたいなんで，
とりあえずまた rsync で復旧．
これで問題なく動く．
-
今日の計算まち時間はとある (短めの)
査読ぎょーむをぼちぼちと進行．
-
3915 秒で終了，
ということで昨日より 1000 秒ほど速くなった．
収束はまるでダメだけど．
-
また策がないまま，
こんどはくだんの無情報事前分布をやめて，
stage 間のばらつきをあらわす事前分布 (正規分布)
とその超事前分布 (ガンマ分布)
を使ったハンパなモデルに変更してみる．
空腹になってきたのでとりあえず 600 MCMC step に変更．
予測される計算時間としては 25 分間ぐらいか?
-
予想どーり 25 分間で終了．
混交はマシになった．
モデルとしてマシなのかどうかは疑問であるけれど．
長時間計算を命じて撤退
……
いや，
待て WinBUGS の挙動がヘンだ．
初期化でコケる．
超事前分布まかせにしてきた分散共分散行列に関して，
てきとーなる初期値を設定してみる．
今度は問題なく MCMC 計算が走りはじめた．
ひょっとしていままでうまくいってなかった原因ってコレ?
だとすると前と同じくまた手ぬき初期化でしくじっていたワケか．
-
ともあれ，
本日は撤退．
1940 研究室発．
1955 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0800):
バゲット．
ネギ・卵の炒めもの．
- 昼 (1240):
タラコスパゲッティー．
- 晩 (2130):
米麦 0.7 合．
ダイコン・ニンジン・ハクサイ・煮干の味噌汁．
2006 年 11 月 27 日 (月)
-
0800 起床．
ごみだし．
朝飯．
コーヒー．
0925 自宅発．
曇．
0940 研究室着．
-
昨晩命じておいた計算だけど
……
収束はややマシになったか?
というような．
しかしまだまだですなあ．
この簡単な問題について，
私がまだまだ理解できてない部分がある，
ということなんだろう．
-
A 棟 6F 実験室に連行されアヤしげ測定機器につながってる PC
の暴走処理
……
これって Windows95 か．
再起動する以外にテはないので，
再起動．
どうやらデータファイルは保存できてたそーで．
[VarioEL]
A-606 実験室の通称 C-N コーダー．
いきなり割算値
(窒素含有割合など)
を吐きだす邪悪な機械．
これまで何人かの院生がこれの扱いに苦闘してきたわけだが
……
私の偏見ではこの機械のおかげで何か現象が解明できた，
という事例はほとんど無いように思えるんだけど．
-
サロベツ ZIP 階層ベイズモデルのほうはよくわからないので，
事前分布から多変量正規分布を除去してみる．
どうなることやら．
-
計算時間 (600 MCMC step) は 1500 秒から 1200 秒ていどに短縮された
……
収束はあいかわらず良くない．
もうちょい長くサンプリングしてみるか?
-
査読報告ひとつ片づいた．
送信．
時刻はすでに 1155 か．
上記の長めサンプリングはいまだ終わらず．
イヤなことに次の査読もある．
締め切りは来週の水曜日なんだけど今から着手しないと破綻するだろう．
こちらはさらにややこしい内容．
来週といえば木曜日夜には東京出張に出発だな．
かなり危ないような気がしてきた．
-
3100 秒 (52 分) ほどかかって 1600 MCMC step 計算終了．
あいかわらずダメ．
line と stage を分割してみるか．
-
600 MCMC step に 1500 秒かかったけど，
やはりダメ．
策がないので，
漫然と長サンプリングを命じておいて
……
-
1300 より 2.5 時間ほど
IPPB 輪読会，
本日は宮田さん担当で Chapter 10 の生活史戦略の進化みたいなところ．
教科書としてはこうやって幕の内弁当的に説明せんといかんのだろうけど
……
今日も Silvertown & Charlesworth てぬき説明とかを
私などが黒板に図を描いたりしながら補足説明してみたんだけど，
こういう進化生物学がらみの個体群動態「理論」って
-
(あまり進化生物学に関係ないヒトたちばかりだけど)
植物生態学を研究している院生たちの
興味や学年にあまり関係なく知られていない
(勉強してない)
-
知らなくてもまあ
(今までの例から憶測すると)
修論・博論で何も困ることなさそう
というような気もする．
そりゃあ，
こういう「理論」にどっぷりひたった古典論文とかに対峙したときとか，
深く読んだりもできないだろうという気もするんだけど
……
うーむ，
現状を要約できないんだけど，
私が院生だったころとは何かが変わっているのかもしれず，
たとえば
(院生の現状とはちょっとまた別に)
こういう「植物集団の数の変化の理解」の最前線と
IPPB みたいな古典的教科書の乖離がより大きくなっている，
というか
……
やはりよくわからん．
-
その間に 1600 step MCMC 計算は 4100 秒を費して終了していた．
そしてあいかわらずトリプルチェイン間の乖離が大きい
(やちパワー関係のパラメーターとかで)．
やっぱり私がでっちあげた統計モデル
(ZIP + line ごとの random effects を階層ベイズモデルで表現)
が不適切なんだろうなあ
……
ということで，
R
に作らせたやちパワー地図などをココロしずかにながめてみる．
-
空間相関はとても気になるけど見えてないふりをするとして，
だ
……
やちパワーのすごく高いところ (赤) にあまり個体がいない，
というのを表現できてないよね．
これはやちぼうず直上でスゲがぼうぼうにはえてるから他の植物は
入り込めないわけだが．
じつは，
やっかいなことにハナシはそう単純ではなく，
ブタナとかは一部が入りこんでいる場合が少数あるんだよね．
モデル単純化を
まんまと阻害してくれた
この悪魔のごとき外来植物め
……
-
そして zero-inflated Poisson (ZIP)
ではなく zero-inflated negative binomial (ZINB)
モデルかなあ，
とか．
うーむ．
-
まあ，
やちパワー計算をひねくるとなると，
手持ちのプログラムではできず，
小山さんのマシンに入ってる R コード群が必要．
で，
小山さんは明日のやちセミナー準備でご多忙なので，
幸か不幸か
私はそれをじゃまできない状況
……
ということで，
これからしばらくはやち計算からはいったん離脱して，
ゆっくりとたてなおしをはかろう．
-
現状のしめきりあれこれを整理．
-
ややめんどう査読:
これは来週の水曜日 (12/6)
-
シウリザクラ原稿したうけ:
これは 12/1 (金)
に提出せよとの永光さん指示あり
-
12/9 (土) 東京
40 分間ひや汗発表:
そうとうにアブないような気もする
-
ホントに心底ど忘れしてたんだけど，
母子里林冠モデルの論文原稿，
さっさと文章修正してアップロードせんといかん
……
やちモデルで楽しく遊んでいてよい時間は，
東京出張がオワるまでは
じつはもうないという状況だったりするのかも．
他にも危険なニオイを発しているあれこれ
……
-
予測不可能な突発的な難題したうけ
-
ポスター発表することになった屋久島計算やりなおしとか
-
アリはどこ?
どこにいったの?
-
1800 研究室発．
どうすればよいのかよくわからなくなった
……
とりあえず散髪．
クラーク会館地下の理髪店．
1845 同発．
閉店まぎわはしあがりが速い，
というべきか．
1900 帰宅．
体重 71.6kg．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0840):
米麦 0.5 合．
ダイコン・ニンジン・ハクサイ・ワカメ・煮干の味噌汁．
- 昼 (1245):
研究室お茶部屋．
食パン．
- 晩 (1930):
うどん．
ダイコン・ニンジン・ハクサイ・ワカメ・煮干の味噌汁．
2006 年 11 月 28 日 (火)
[LaTeX な表]
R
の
xtable
package
使うとこういう表が生成される
……
問題は，
LaTeX とか使わぬヒトたちとの共同作業時ですなぁ
……
というか表だけでなく数式の変換もたいへんそう．
うう．
-
で，
他の図のつくりなおしとか終わったのかしらん，
と 7F に様子を見にいったんだが
……
本日は不在のようで．
ということで，
ここまで書いた文章のみなおしとか．
-
みなおしやったことにして，
本日のシウリザクラしたうけ終了．
-
めんどうな予感のする査読の論文原稿をぱらぱらとながめる．
1935 研究室発．
1950 帰宅．
晩飯．
-
ひきつづき，
ぱらぱらとながめつつあれこれと考えてみる．
やはりけっこうめんどうだ．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0750):
バゲット．
- 昼 (1330):
研究室お茶部屋．
北大生協の巻きずし．
258 円．
- 晩 (2130):
米麦 0.7 合．
ダイコン・ニンジン・ハクサイ・ワカメ・煮干の味噌汁．
2006 年 11 月 29 日 (水)
-
0720 起床．
朝飯．
コーヒー．
0830 自宅発．
晴．
昨晩はけっこう雪がふったもよう．
0845 研究室着．
-
ネット雑用．
-
ついでにというか，
なけなしのきあいを発動して，
香川の
伊藤さん
にアリデータに関するご相談メイルを書いて送信してしまう．
はてさて，
どうなることやら．
-
めんどう査読のつづき．
何がめんどうかっていうと
「簡単に解決できる問題をわざと難しくして解いてみせる」
というハナシなんだよね
……
たとえて言えば，
札幌から東京に出張する状況で
「移動時間を短くすることは重要」
とか言ってるのに飛行機つかわずに実現しうる最短到達時間の計算
→
「鉄道・道路の人身事故頻度がこれこれの条件のもとでは
複数の最短解が同時に存在可能です!」，
てなことをやってるよーな研究かな
……
複数あったからといって何がうれしいかとか，
そもそも
なぜ飛行機を使ってはいけないのかといった理由は示されていない．
-
私はふだんの査読では解決される問題の価値についてはあれこれ言わずに，
解決方法がきちんとしているかどうかだけを
入念にチェックすることにしてるんだが
(いいかげんな論文原稿はだいたいにおいて，
こちらがまずコワれているから)
……
今回のはちょっと，
という気もしてきた．
計算だけはちゃんとしてるようには見えるんだよねえ．
-
論文原稿を引用文献などながめつつ読みススめていると，
何やら「助教」がらみのメイルが事務から送りつけられてくる．
その抜粋．
授業が増えそうなのはしょうがないとして，
教授会とかに時間とられたり雑用
(「当研究院の運営に責任を持って協力」)
をおしつけられたりするのはイヤだなぁ．
そしてまた書類かきか
……
-
すでに正午．
寒冷前線通過中でみぞれがふっている．
本日は北大構内走は無理そう．
あ．
雪になった．
-
たんたんと査読つづけている．
文章はすごくわかりやすい
……
おっと，
計算まちがいというか定式化のまちがいみつけた．
本人は単なる手ぬきのつもりかもしれんけど，
結果はかなり変わるはず．
とりあえず，
こういうわかりやすいまちがいは指摘しやすいよね．
ついでに説明てぬきも発見．
しかしこれは「言い換え」で何とでもなるな．
しかしながらこっちは
……
もとネタ論文を読んでみたんだけど，
そんなおかしな仮定はもうけられないと思うんですが．
-
なんとなく査読報告かきのねたはそろってきたような気がする．
休憩して昼飯．
-
査読報告かきにとりくんでみる．
また時間かかりそう．
しかし今日中に終わらせたいんだけど
……
-
なかなかススまず，
お茶部屋に逃げる．
report 書きがなかなか進捗しない大学院生と現状に関する不満を
相互に報告する．
-
じりじりと書き進めていると，
あやしげ統計モデリングメイルを受けとったので
ゴミ箱に捨てようかと思ったけどとりあえず放置
……
またお茶部屋に逃げたりすると，
R
作図で苦闘している院生がいたので
現実逃避的「R で図のつくりかた」講習会．
-
席を離れているあいだに親切なだれかが作文をススめてくれている
……
はずもなく，
また作文のつづき．
-
また少し書きすすめてから，
またお茶部屋に逃げると report 苦闘の院生から酢をススめられたので，
二人してそれを飲む．
-
……
1900 ごろ，
なんとかひととーり書けた．
みなおし．
ひととーりみなおした．
プリントアウト．
プリンターは院生部屋にある．
そのままプリンター近傍の院生のところで
R 作図講習．
30 分ほどでそれを終えてお茶部屋にふらふらとさまよいこむと，
そこで展開されていた
お茶部屋雑談ストリーム奔流の中でくるくる回されて呆然としてるうちに
「ぶらうざー係」
命じられた．
これはお茶部屋共用 Mac の web browser を起動し，
院生たちの
「アレってどうだっけ?」
「反町くんの最近の仕事って
……」
「ほら久保さん，
アレ検索しなさいよ，
アレ」
「そっちじゃなくてその上!
……
こりゃダメだ，
もっとちゃんと説明してるところ!」
といった著しく抽象化されたコマンドを
前後の (ぱられるちぇいんな) 文脈にそって適切に解釈して
計算機に理解できるカタチで伝達してネットの海から情報を拾いあげる，
といった高度に知的な能力を要求されるお茶部屋社会最下層階級の仕事である．
-
で，
お茶部屋を脱出して査読報告の最後のみなおし．
2040 送信，
これにて一件らくちゃく〜，
と帰宅しようとするとゐんどーづに呪われた院生からおはらいの依頼
……
TCP/IP ネットワークぷろとこるが使えません，
と．
可能な原因をひとつずつ消去していくと
……
どうやら Windows Update で OS 書き換えられたんだけど，
げいつ OS 必須の「再起動」やってなかったことが原因とわかった
……
いやはやー．
たぶん Vista とやらになってもこの調子なんだろね．
-
2130 研究室発．
路面白凍結．
気温氷点下．
買いもの．
2155 帰宅．
晩飯．
-
一昨日の
「現状のしめきりあれこれを整理」
の現状
-
ややめんどう査読:
これは来週の水曜日 (12/6)
-
シウリザクラ原稿したうけ:
これは 12/1 (金)
に提出せよとの永光さん指示あり
-
うーむ
……
花粉実験まわりの作文は
とりあえず昨日なんとかしてみたわけだが
……
他の図とか，
いったいどうなってるんだろうね
-
12/9 (土) 東京
40 分間ひや汗発表:
そうとうにアブないような気もする
-
ホントに心底ど忘れしてたんだけど，
母子里林冠モデルの論文原稿，
さっさと文章修正してアップロードせんといかん
-
[今日の運動]
-
[今日の食卓]
- 朝 (0810):
米麦 0.5 合．
ネギ卵炒飯．
海藻スープ．
- 昼 (1340):
研究室お茶部屋．
食パン．
- 晩 (2220):
米麦 0.7 合．
大家さんにめぐんでもらったロールキャベツ．
2006 年 11 月 30 日 (木)
-
0705 起床．
朝飯．
コーヒー．
0845 自宅発．
晴．
すでに毎日の最低気温は氷点下になっている．
0900 研究室着．
-
A804 室廊下の天井で何やら工事．
けっこう騒音が．
-
シウリザクラ図
……
なるほど，
「できてません」
か．
まあ明日あたりからの展開次第だな．
今日あわてて強制執行することはない．
-
ということで，
来週の「アブない」東京出張の準備に着手してみる．
具体的には 40 分間ほど R
と
MCMC 計算
について話さないといけない．
いろいろと難しいところがあるんだけど
……
たとえば聴衆の構成．
おそらく以下のように極端なヒトたちからなっているだろう:
-
このあたりむちゃくちゃに詳しい
(もちろん私などより遥かによく理解している)
ヒトたち
-
「MCMC 計算って何?」
なヒトたち
どっちが多数派なのか判然としないのが
この集会のオソろしいところなんだが
……
まあ，
2. 方面にあわせて準備するのが妥当であろう．
-
12/9 (土) の日本語版集会で他に発表される
ヒトたち
って，
generic というかどういう聴き手にもウケそうなハナシだよな
……
私のだけが中途半端な気がする．
-
うーむ，
終了点 → 出発点方向で考えてみようか．
全体 40 分間の中の 25-30 分めあたりから
最後のほうの 10-15 分間ぐらいを使って何を話していればよいのかといえば，
おそらく
「まあやっぱり WinBUGS と
R2WinBUGS
が今のところは便利なんで，
しばらくはコレをこう使ってみますかね」
といった内容だろう．
-
とするとその前は「有力な対抗馬かつ R と縁の深い」
JAGS
について紹介しとかんといかんだろうし
(何しろ「R の整備と利用」研究会なのだ!)，
OpenBUGS についても言及せんといかんだろう．
OpenBUGS といえば，
lindley さんの
R と WinBUGS
における OpenBUGS まわりの解明はすごいよね
(R2WinBUGS 改造わざとかもスゴい)．
-
で，
こういった *GS なソフトウェアを導入せねばならん理由をあげねばならず，
それは
「R
単体では MCMC 計算やるのはなかなかしんどいから」
と言いつつも
library(MCMCpack)
とかにも言及せざるをえなくなる，
と．
すると今年 3 月発行の
Rnews
Volume 6/1 (べいづ特集号) とかにもふれる必要があり，
WinBUGS が優占種になってしまってる世界で
どうして OpenBUGS や JAGS が必要かといったハナシもありかもね．
-
ここまでのハナシで 40 分間のうち 30 分ぐらいは使っているかもしれん．
というのも，
R & BUGS のコードをみせたり，
実演とかもやらんといかんだろうから．
-
とすると，
「残された最初の 10 分間」
(10 分間もないかも?)
では何をしゃべるべきか?
というと，
まあ「MCMC 計算って何?」という聴き手もそれなりにいると想定して，
「MCMC 計算はベイズ推定に必要」
「ベイズ推定とは，
階層ベイズモデルとは，
GLMM
(経験ベイズ法)
とのちがいは」
といったハナシになるだろう
……
やっぱり今回も
「簡単な例題を解いていきつつ手法を紹介する」
というわくぐみでやるか?
-
まあ一月末の授業のうち一回はそういうハナシだろうし，
(もし採択されれば)
三月の生態学会松山大会の自由集会もそういうハナシだろうから，
今のうちから準備しとくのは悪くない．
授業・自由集会と今回のちがいは何か，
と考えると
……
今回のが一番
「ソフトウェア寄り，
考えかた省略」
というかんじだな．
M1 あいての授業とかでは
MCMCpack
なんて絶対に紹介しないよ (ホントか?)．
-
使う例題とかまた
架空植物の結実
とかだろうなあ
……
おっと，
この生態学会誌の解説記事，
新しくなった glmmML()
に対応させなくては．
-
食材ローテイションにしくじり，
いったん帰宅して昼飯となる．
北大構内走向けの好天なのにもったいない．
1230 帰宅．
ついでに洗濯．
昼飯の準備．
昼飯．
1340 研究室もどる．
いつのまにか雪．
-
和文誌の
GLMM 解説記事 web 補足
に
glmmML(cbind(Yes, No) ~ ...)
記法に関する追記，
と．
-
窓の外はすごい雪．
冬型気圧配置．
-
文部科学省共済組合北海道大学支部からナゾの紙切れ．
医療費の合計，
とかで 6 月に通院してた歯医者に合計 23520 円だそーで．
-
めんどくさそうな作業をまず片づけるか，
ということで，
シウリザクラ下うけ作文の私のここまでの分担部分の見なおし
&
TeX2Word
をもちいたわーど化．
-
……
の作業には着手できず，
松田さんと駒ヶ岳データを
R
で楽しくデータ解析．
かなりトリッキーな
data.frame()
変換が必要とされた．
-
もちろん
R
があれば何とでもなるわけだが
……
いかに R が「データまちがい探し」においても
きわめて強力なツールであるとはいえ，
事後的にまちがいを search & destroy
するのがいつも正しい方策であるのか，
と．
-
代替案としては入力段階でチェックすることで，
より具体的には
-
ゑくせる上で入力時まちがいのチェック:
VBA によるプログラミングが必要だろう
(かなり面倒そう)
-
web interface とかをもつ何かリレイショナルデータベイス:
チェック部分は Perl とか?
-
R そのものを入力ソフトウェアとする
じつは 3 が案外妥当なのでは，
という気がしてきた．
「データ入力はゑくせるでなきゃダメだ!」
という院生たちのゑくせる信仰をくつがえすような
使いやすい interface (もちろん非ゑくせる的な)
を設計するのはたいへんそうだけど．
-
1630 より
Trendy セミナー，
本日は 崎尾均さん
(埼玉県農林総合研究センター)．
秩父山地の渓畔林の林冠木である
シオジ，サワグルミ，カツラの特徴を抽出するような調査・実験いろいろ．
どのように渓谷内の撹乱に依存しているか，
三者三様のやりくちがあるらしい，
ということが示されていく．
-
いろいろ紹介された中で，
わりと細かい問題が気になった．
シオジの結実はあるていど広いの個体群内で同調しており，
3-4 年に一度ぐらいは落下種子量ゼロとなる年がある．
種子生産の個体群内同調は
井鷺さん
のいわゆるタンクモデルで説明できるだろう．
-
それはよいのだけど
……
ここで奇妙なのは，
シオジは雌雄異種であり，
開花の同調はメスだけでなくオスででも見られるんだよね．
オスの個体では種子を作らないんだけど，
花を作り花粉をばらまくだけ．
そして井鷺さんのモデルではオスの花の同調開花は説明できないんだよねえ
……
-
崎尾さんからうかがったあれこれを列挙しておこう．
気象要因ではこの同調は説明できなさそう
(key となる気象要因は見つけられなかった)．
シオジの性比はだいたい 1:1 (おそらく成木の)．
サイズなどによる性転換なし．
花粉はたぶん風で運ばれているのだろう，
と．
-
セミナー終了後，
また駒ヶ岳データにとりくむ．
なんとかなったので
1930 研究室発．
1945 帰宅．
体重 71.4kg．
晩飯．
-
なぜか今宵は統計メイルがたくさんきてるので (masting 日?)，
いろいろ返信してみる．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0750):
米麦 0.6 合．
味噌汁．
- 昼 (1310):
米麦 0.6 合．
鶏レバ・ショウガ炊きこみ飯．
たまに鶏レバとか買ったりすると，
こう調理のタイミングを逸してですねえ
……
しかし考えてみれば，
昨晩のうちに加熱調理して煮物にでもしてしまえばよかった．
- 晩 (2030):
米麦 0.7 合．
鶏レバ・ショウガ炊きこみ飯．
コマツナあえもの．
