ぎょーむ日誌 2005-12-(21-31)
2005 年 12 月 21 日 (水)
-
0840 起床．
朝飯．
コーヒー．
0945 自宅発．
晴．
1000 研究室着．
-
昨日の北大ネット途絶に関して
……
管理者連中はまた「なかったこと」にしてるよーで．
あるいはホントに途絶したことに気づいてなかったとか?
-
苫小牧直径成長問題，
樹木個体 random effects の積分のやりかたを変える，
とゆー昨日の改造プログラムの試運転．
最尤推定はうまくできるけど，
事後分布の作図がうまくいかない．
試行錯誤しつつばぐとり．
正午前にまっとうそうな結果が得られるようになった．
-
こんなかんぢで．
二ヶ月前の 10/27
ごろに計算してたのとほぼ同じなんだけど，
事前分布の「下の端」まで計算するようになった，
というのが相違点である．
-
1130 北大ネットまたしても原因不明の通信途絶．
少なくともこの建物 (地環研 A 棟)
は完全に孤立している．
A 棟 hub が「ツマってる」のか，
ということで 3F まで行ってまた無断強制再起動を試みたんだが
……
ぜんぜん効果ナシ．
ということで北大情報基盤センターとやらが管理している
地環研ルーターかそれより上流
(ただし北大内)
のどこかが壊れているもよう．
1157 ごろ一瞬だけ外界につながったけどすぐに沈黙．
-
昼飯．
1300 すぎてもネットはネット復旧せず．
-
かとーさんが i-mode などで調査し考察したところ，
どうも A 棟だけがネットから分断されていて，
その原因は
(以前にも同じ事例があったんだが)
A 棟内の
virus 呪われゐんどーづ機暴走なんかにあるんでは，
とのこと．
そこでまた棟 hub を無断再起動してみたんだけど，
状況は変わらず．
-
しょうがないんで内線 4494 でネットワーク掛に電話．
状況説明．
やはり先方は管理者のくせにこの通信途絶に気づいてなかった．
えーと，
ここ 8F の hub 番号は
T2208，
使用ポート番号の例としては 02
など．
調査するとのこと．
-
しばらくしてから脱力返答．
通信不能なので状況わからず．
ぎょーしゃさんとやらを現場に派遣する．
時間かかるけどそれでいいか?
……
これって「だめだ」といったらもっと迅速に復旧してくれるわけ?
-
ぼーっとしてるのも何なのでかとーさんと調査つづける．
どうも A 棟 hub と 8F hub のどちらかに問題あるらしい．
8F 闇ネット内には問題なし，
と実験によって確認．
となると，
当方には打つテなし，
か．
-
1440 きゃたつかかえた「ぎょーしゃさん」とやらが登場．
北海道ゼロックスから二名きた．
ThinkPad X40 の USB から 8F hub の RS232C シリアルポートに直結
……
するまでもなく uplink ぎれ，
とわかる
(私はそのあたり注意してなかった)．
で，
一人が 3F hub まで行ったところ，
8F につながる (12 番ポートの) ケイブルがぬけかけていた，
とのこと．
さしなおしたら復旧．
誰がやったんだ?
-
ともあれ北大ネット途絶時の点検手順を更新．
-
かとーオフィス内の断線を調べる
-
A 棟 8F 内の通信状況を調べる
-
A 棟 8F の uplink を調べる
-
3F hub の 12 番ポート
のケイブルを点検
-
3F hub の電源再投入
-
それでも動かなければ
情報基盤センター 4494 に電話する
というところか．
-
なぜかしら間瀬さんが
RjpWiki
内に
Gmail 宣伝
ペイジを作っておられる．
先日
このへんで書いた ``View as HTML'' 機能紹介しました，
とわざわざ連絡をくださったので
……
-
苫小牧樹木直径成長問題，
本日えられた結果などアップロードしてみる．
「個体差」 (random effects)
をあらわす事前分布のいれかた
(積分のやりかた)
を変えると推定結果がびみょーに変わった．
全体の最大化対数尤度も改善しているので，
喜ぶべきことなのかもしれないが
……
まあ，
ひじょーにびみょーな差位であるし
……
-
イタヤカエデ (AM) なんかにおいて，
大サイズで成長が「へたれる」ようになった，
といったところが見てくれに関して変わったところ．
これは「私は集団内のすごくへっぽこ樹木個体です」と自称する連中の
「個性」を (雑ながら) 端まで計算できるようになったため
……
と考えたいのだが，
あまりにも苫小牧ボスの意向にそった推定結果なので，
改造点を中心に計算プログラムの再点検にとりくまねばなるまい．
というのも，
私が「無意識」のうちにボスの意向に迎合する計算結果を出すような
コード書きをやってた可能性が
……
-
変数名改善などもかねて，
一行ずつ点検．
変数名を変えたので試験運転．
どちらも問題ナシ，
のように見える．
そもそも，
すごく単純な改造だし．
-
ということで，
じつは大サイズではやっぱり落ちていたのか．
成長速度．
-
こうなってくると，
例の hyperspecies / species / individual
三階層ベイズ推定とかやったらどうなるんだろう，
と好奇心わくわけだが，
これはすけぢゅーる破綻という墓穴への直通経路になっているんで，
今回はあきらめることに．
-
こちらのモデリングに関してはボス指示待ちだな．
うかつに動けん．
-
1910 研究室発．
買いもの．
1935 帰宅．
晩飯．
-
直径成長モデル，
さっそく質問いただく．
てきとーに回答しようとして
……
ふと手が止まる．
そういや，
random effects の積分が少し手ぬきぎみだったかも．
予想としてはほとんど影響ないはずなんだが，
だんだんコワくなってきた．
ということで，
とりあえずまずはそちらを点検．
dx.re <- 0.05
x.re <- seq(dx.re * 0.5, 6.0, dx.re)
sum(dgamma(
x.re,
scale = scale,
shape = 1.0 / scale # because mean == 1
)) * dx.re
[いいかげんな積分?]
いいかげんに離散化したガンマ分布を積分した値．
scale という分散を定めるパラメーター
がでかいところでは 1 からのずれが大きくなる
(分散 > 平均
2 な部分)．
まあ，
不幸中のさいわいでこんなに「個体差」分散が大きいのは
オオモミジ (AP) だけ，
なんだけど
……
まあ，
きちんと補正するようにしよう．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0900):
米麦 0.7 合．
ダイコン・ニンジン・ネギ・ブナシメジ・タラすりみの味噌汁．
- 昼 (1240):
研究室お茶部屋．
ツナマヨパン．
- 晩 (2030):
米麦 0.7 合．
ダイコン・ネギ・チカの味噌汁．
ホウレンソウ卵とじ．
2005 年 12 月 22 日 (木)
-
0730 起床．
本日は冬至か．
朝飯．
コーヒー．
-
朝から
……
というか朝ならではのみょーな問題にとらっぷされる．
屋久島階層ベイズモデルのほう．
パラメーターの事後分布なんかに関しては
このへんで図示してるとーりなんだけど，
「ある現象の発生する事後確率」
はどうなるんだろう，
といったわけのわからぬもの．
観察された現象の尤度の事後分布は計算されている．
しかし，
「観察されなかった」
現象の尤度は?
各 step ごとに素朴に尤度計算やりなおせばいいだけ，
なんだろうか．
-
考えてもよくわからんので
1010 自宅発．
晴．
まあ，
修論したうけぎょーむという観点からは，
この問題はしばらく放置してさしつかえないんだが，
と自分を納得させようとしてみる．
1025 研究室着．
-
いきなりネット途絶?
しかしこれはなぜかしら，
かとーオフィス内の switiching hub が
「つまっていた」
ためとわかった．
-
苫小牧樹木直径成長モデリング．
昨晩確認していた積分の補正とか．
あらら，
事後分布計算まわりに少しバグがあるような気がする．
いろいろ修繕して計算やりなおし．
-
う．
推定結果，
けっこうかわってしまった
(しかも冷汗ながれる方向に)．
再点検しなくては．
-
かなりうろたえつつメイルやりとり．
-
R
の
data.frame()
の仕様が 2.2.0 (?)
ぐらいから変化していたらしい，
とようやく気づいた．
以前の仕様に依存してプログラム作ってたので，
これは改造する必要ある．
-
以前の
data.frame()
に label つき vector を列として追加すると，
label までいっしょに格納されていた．
しかしながら現在の version は label がはぎとられ，
値だけが保存される．
これは data.frame のサイズを小さくするため，
というような説明を以前どこかで読んだような気がする．
-
人類の歴史上おそらく最も懇切丁寧なる R 本，
データ解析環境「R」，
送っていただいた．
ありがとうございます．
[データ解析環境「R」]
なぜか工学社．
著者たちによる
サポート頁
もすでに完備．
著者の一人・舟尾さんに教えていただいたところによると
R-Tips
が (あんなに丁寧なのに) まだ難しい，
という意見があって，
このさらに易しい本が作られたとのこと．
-
1600-1950，
塩寺さんデータすなわちインドネシア山地熱帯林葉っぱデータ解析のこんさる．
屋久島につづいてまたしても葉っぱですよ．
けっきょく安易な回避は見つからず，
屋久島問題と同じ方向 (葉の死亡確率は何で決まるか
→
推定された葉齢依存死亡確率から葉寿命分布の計算)
になってしまった．
ただし多正面作戦は私にとってまったく自滅的なので，
葉寿命ぎょーかいにたむろしてるヒトたちのやっていること
(の根底にある暗黙の仮定もろもろ)
が，
あたかも正しいものであるかのようなふりをして，
個体差もなければ種差もナシ，
という計算方法をまずは説明してみる．
-
常緑樹の葉っぱ寿命計算 (というか葉寿命なるものの定義ですな)
ってのはけっこうぼろぼろというか個体群生態学的に見ると
「なんぢゃこれは」
ものが多い．
菊沢さん「葉寿命」本の (季節性よくわからぬ)
常緑植物の葉っぱデータにおける葉寿命計算の説明をみて，
ぱっと理解できるヒトは相当にアタマおかしいと認定されてもよい．
長年放置されてた生ごみバケツのフタあけて顔面をつっこんだようなかんぢだ．
-
2005 研究室発．
2030 帰宅．
晩飯．
苫小牧ボスから計算指令が
……
-
[今日の運動]
-
[今日の食卓]
- 朝 (0820):
米麦 0.7 合．
ダイコン・ネギ・チカの味噌汁．
- 昼 (1250):
研究室お茶部屋．
米麦 0.5 合．
ホウレンソウ卵とじ．
- 晩 (2130):
米麦 0.7 合．
ハクサイ・ネギ・豚ハツの炒めもの．
ダイコン・ネギ・チカの味噌汁．
2005 年 12 月 23 日 (金)
-
0900 起床．
こんさるばて，
か?
そして本日は祝日である，
としばらくしてから気づいた．
朝飯．
コーヒー．
1030 自宅発．
晴．
1045 研究室着．
-
メイリングリスト
R-help
への質問投稿
「
ls() で見えない関数を作る方法は?」
に対する常連 Martin Maechler 氏の
回答:
-
package にいれる
(namespace の利用)
-
.abc()
というようにドットつき関数名にする
どちらもおもしろい．
-
苫小牧樹木直径成長モデル，
glm(..., family = quasipoisson)
と GLMM のこんびねいしょんワザで段階的に
気象フィルターを順位づけ・選抜する方式に切り替えたので
(昨年 11 月末にはこの方式を試していたんだけど，
おくらいりになっていた)，
気象フィルターの候補を増やすということに
……
どうなることやら．
ボスを立腹させる結果にならなければよいのだが．
-
これまで
前年気象値 54
×
前年気象値 27 = 1458 とーりで探索してたんだけど，
今回は前年気象値 256
×
前年気象値 64 = 16384 とーりで試みてみる．
11 倍に増えたわけだが
……
イタヤカエデの計算時間は 6 分から 12 分に増大，
か．
-
昼飯．
-
苫小牧樹木直径成長モデル，
推定計算の拡大版ぢりぢりと進行中．
ひとつの気象値セットの計算時間は 47-49 分間ぐらいか
(昨日までの試算版だと 6-7 分)．
ここでいう気象値セットとは，
「雨ふりがおよぶ効果」
パラメーターを固定した 16384 セットの計算である．
この雨パラメーターは全樹種共通にしていて，
これまた何とおりか試してみなければならない．
-
計算は計算機にまかせて，
人間は他の仕事をススめればよいはず．
ところがススまない．
不思議だ．
-
1750 本日の計算，
ひととーり終了．
今回はボス追従的な結果になった．
計算者のココロの状態が計算に反映されてしまうのだろうか．
-
といったえせ深層心理学的なハナシはおくとして
……
計算結果は幹部はここ一年ぐらいほとんど変わっていない，
と言ってもよいのではないか?
ただし「推定結果の末端」とでもいうべき部分は「ふるえている」ように見え，
おそらく統計学的な意味においても
定まらない状態にあるのだと考えるべきかもしれない
……
と，
アタマの中が徐々にべいぢあんな方向に偏向しつつあるわけだが．
-
このようにふらつく原因として
-
線形予測子だの link 関数だのは
生物学の何かの原理原則
(そんなもんどこにもありませんね)
から導出されたものではなく，
まあ解析者たる当方に都合のよい近似にすぎないので
「つねにあてはまりが悪い」
とでもいうような状態にある
-
さらに通年毎日の気象データなる莫大な数値の羅列から，
生物の挙動を説明できるような
「意味のある情報」
をすくいあげることの難しさ
後者に関して補足するなら
……
たとえば気温データ 10 年分と称するじぐざぐの波形どもを渡されたときに，
樹種 X の個体たち 10 年間の成長量の年変動を「説明」できる
気象フィルターセットを構成できるか，
という問題を考えてみよ．
-
回答は簡単で，
それはいくらでも構成できる，
ということになる．
このときに「○○はするな!」
といった制約条件を加えていくことで，
良さそうな気象フィルターセットの候補の個数を多少は減らせる．
というかこれは必要なことで，
もしホントに気象フィルターを自由に選ばせてみたら，
人間などが考えだすよりよほど
「独創的な解」
を最尤推定数値解法プログラムが探しあててくれることにあなたは驚くだろう．
-
しかしながら，
この制約を加えるという行為は計算するまえに
その結果を決めてしまうようなもので別の意味で危険きわまりない．
で，
制約はゆるめにしなければならないのだけど，
そうするとまったく異なる気象フィルターセットたちが，
ほとんど同じような
(むろんびみょーに異なるけど)
年変動パターンを生成しうるらしい
……
というのがよーやくにして私などのアタマに浸透しつつある，
この現象論的データ解析の実相の近似的理解である．
-
強引な制約をつけられるほどには，
われわれは植物内部の挙動を理解できていない．
つまりこの強引な制約とは，
気象データセットなるものを「絶対値」の集まりとして解釈できるような，
定量的な植物生理学の知識である．
幸か不幸か，
そんなものは誰も知らない．
-
べいづ推定でいえば，
すごくカタチの制約された分散の小さい事前分布みたいなもんかな．
-
うーむ，
とアタマをひねりつつ撤退．
1930 研究室発．
2030 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0940):
米麦 0.6 合．
ハクサイ・ネギ・豚ハツの炒めもの．
- 昼 (1320):
研究室お茶部屋．
米麦 0.6 合．
ハクサイ・ネギ・豚ハツの炒めもの．
- 晩 (2100):
米麦 0.8 合．
ネギ卵炒飯．
麻婆豆腐．
2005 年 12 月 24 日 (土)
-
0800 起床．
朝飯．
コーヒー．
たまってた洗濯．
怠業．
-
昨晩，
北大ちかくの古本屋で
Greg Eagan
(SF 作家)
の短編集二冊を買った．
今日は「祈りの海」のほうにはまりこんでしまった．
-
アイデアねたとしては情報科学と生物学方面が多い．
たとえば
「ミトコンドリア・イヴ」
はミトコンドリアのプラスミドの塩基配列だけではなく
EPR 相関
が観測できるとすればというホラ話から，
それがひきおこす社会のどたばたにつなげている．
怪獣つくったりナノテク汚染やったりができないデータ解析者にすぎぬ
われわれもまた，
この世界を混乱のどんぞこにたたきこむ，
あのあこがれの
「マッドサイエンティスト」
になれるんだ，
と勇気づけられる一篇である．
-
1650 自宅発．
とっくに日没．
雪．
大学に向かう
……
が「手ぶら」で登校したら大学院生に怒られそうな気もしてきたので，
北 12 生協でお菓子など少々買ってから
1720 研究室着．
-
とりあえず，
ぢりぢりとメイル書きとか，
計算プログラム & データなどのアップロード．
-
生態学会さーばー雑用．
ユーザー追加する場合，
あれこれと手作業が必要な場合がある．
symbolik link 二箇所を設定しなおすだけだが．
-
屋久島葉っぱ寿命モデリング復帰にむけてアタマをうごかす
……
やはり，
survival analysis のワザは流用できん，
か．
-
2150 研究室発．
2230 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0900):
バターロール．
- 昼 (1320):
スパゲッティー．
麻婆豆腐．
- 晩 (2330):
蕎麦．
エノキダケ卵とじ．
2005 年 12 月 25 日 (日)
-
0850 起床．
朝飯．
コーヒー．
怠業．
-
Ruby 開発者の
まつもとさんも
ThinkPad X31
に呪われてるよーで
(症状が当方のそれに似ているような気がする)
……
私は製造元製の AC アダプターやめて
サードパーティーのものにしてからこういう
「ふりーづ」
呪われからまぬがれているような気がする．
なんか X31 って欠陥製品のように思えてきた
(今は順調に動いているけれど)．
-
1245 自宅発北大構内走．
晴．
雪上走してるヒト，
ちらほらと．
かなり歩く．
1345 帰宅．
体重 73.4kg．
昼飯．
-
屋久島モデル，
シュート伸長休眠のほうをちょっといじる．
-
1535 自宅発．
曇．
1550 研究室着．
-
先ほど手を加えた屋久島シュート伸長休眠の
階層ベイズモデルの MCMC 計算プログラムを動かしてみる．
Dell 機上で 21 分間．
今まではある step 内で個体 → 樹種 → hyperspc
という順にメトロポリスしてたんだけど，
これを逆にしてみた．
やはりというか，
あたりまえのことなんだけど，
計算結果は何も変わらなかった．
-
何がやりたかったかというと，
ひとつの MCMC step として出力されるパラメーター値が
こっちのほうがわかりやすいんではないかな，
というよーな
(上の図ではそのあたりがまったく見えてないし，
後述するようにわりとどうでもいいハナシなんだけど)
……
もちろん MCMC 計算の主眼は「同時分布」
(joint distribution)
からの標本らしきものを生成していくことにあるので，
どのパラメーターが先であとだというハナシにはあまり意味ないわけで．
-
しかしながら，
この問題のように hyperspc の値は各個体のパラメーターをふりまわすけど，
個々の個体の値の増減はそれほど hyperspc には影響しない，
という非対称性があるときには，
あたかも hyperspc → 樹種 → 個体という順で
パラメーター生成させて，
これを MCMC の 1 step とくくるほうが
……
いや，
まあどっちでも同じかな．
とはいえ，
現在のやりかたのほうが，
「『下』は『上』にふりまわされる」
感が少しばかり強調されてるような．
しょせんはキモチの問題にすぎないのだけど．
-
で，
何が気になってたのかとゆーと
……
上の図で見えてるのは周辺分布 (marginal distribution)
なんだけど，
これって同時分布のよい「断面図」になってんのかしらん，
というあたり．
今回の計算結果で作図してみると
……
まあ，
「『下』は『上』にふりまわされる」
んだろうけど，
思ったより hyperspc - spc 間で相関はヨワくてひと安心，
と．
-
なるほど，
こういう階層モデルって，
負の相関があるのがたぶん正しい推定結果なんだろう．
正の相関だったらむちゃくちゃ，
というか．
-
ついでに樹種-個体間でも同様の作図してみたんだが
……
まあ，
これまた当然ながらというか，
より「無相関」な方向で．
「つよい個体差」が言えるほど個体データがないためだろう．
-
お茶部屋で修論さっさと片づけてスキーに行きたい，
というハナシを拝聴する．
そう，
修論はさっさと
……
-
屋久島葉寿命の推定プログラム，
すすまん．
-
理由のひとつは「年枝」のデータ構造を決めてないから，
というあたりのようで．
全体 / 樹種 / 個体 / 年枝 (先端) / 年枝 (次) / 年枝 (...) / ...
とゆー構造なんだよね．
どうしたもんかな．
ばか正直にこのとおりにやるべきかな．
-
年枝がもつべきデータは
……
-
見かけの年齢
(観測)
-
ホントの年齢
(先端方向の休眠・二度伸びに左右される)
-
葉数
(観測)
-
葉痕数
(観測)
-
先端方向の節における休眠回数
-
根元方向の節における二度伸び回数
-
根元方向の年枝への reference
-
先端方向の年枝への reference
(いらんか?)
さらに欠測末端年枝なるものも必要で
(継承するか?)，
上の属性に加えて
が必要とされる．
うへー，
めんどくさい．
「尤度」
は持つ必要あるのか?
計算時間をケチるなら持たせるべきなんだが
……
いや，
あほらしいからやめておこう．
-
で，
Segment
なるクラスを独立させて定義していくわけだが
……
けっこうあれこれと動作を定義してやらねばならぬ．
-
かなり寒くなってきたので撤退．
2100 研究室発．
大雪．
2120 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0900):
バターロール．
- 昼 (1420):
米麦 0.8 合．
納豆．
エノキダケ卵とじ．
- 晩 (2230):
米麦 0.8 合．
ハクサイ・ニンジン・豆腐のカレー．
2005 年 12 月 26 日 (月)
-
0800 起床．
朝飯．
コーヒー．
屋久島シュート伸びちぢみ問題んアタマが拘束される．
1010 自宅発．
曇．
歩いてるうちに
「休眠だのなんだのは『隠れ状態』なんだから
Gibbs sampler とかわりと導入しやすいかな?」
とか
「休眠と二度伸びがあるから二次元幾何分布か?
いや，
休眠と二度伸びは排他的だから一次元が二本，
か」
などなど考えがまとまってくる．
1025 研究室着．
-
うーむ，
苫小牧ボスからのご指示．
これは作図関数まわりを再整備せんといかんのかしらん．
どうしたもんかなぁ．
ちょっと以前の清書作図関数どもを調べてみよう．
一件落着への道はいまなお遠し．
-
で，
この苫小牧モデルまわりの作図関数は設計のできがよろしくなくて
……
と改造に苦闘する．
自業自得というやつか．
-
だんだん抜本的に作りなおしてみたくなってきた．
昼飯．
-
1300 から 1.5 時間ほど，
研究室セミナー補足篇．
先日の M2 おいこみから逃れた石橋君でミズキのハナシ．
ミズキとかヤマボウシでは横枝が「展伸」するのが特徴的なんだけど，
これがミズキ全体のカタチとどう関係するかよくわからない
(別に展伸などやらなくても横枝は当然ながら伸びるものなんで)．
そして枝の垂直配置を決定しているらしい幹の挙動がよくわからない．
-
苫小牧直径モデル EPS 作図関数の整備．
われながらヒドいプログラムというか，
R
の
env = .GlobalEnv
汚染されまくり，
とでもいうようなアタマ痛い手口など使っていて，
たいへんわかりにくく改造しづらい．
-
この
.GlobalEnv
呪われのせいで，
計算結果があわなくてアセる．
修復に手間どる
……
-
すっかり夜になってからよーやく終わった，
とアップロード
……
したんだけど小バグが残ってて，
翌日また冷汗ながすことに．
-
屋久島問題もちゃんと考えてますということをアピールするために，
牛原さんのところで作図 R プログラミングおてつだい．
2030 研究室発．
2045 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0830):
米麦 0.7 合．
ハクサイ・ニンジン・豆腐のカレー．
- 昼 (1240):
研究室お茶部屋．
米麦 0.5 合．
ハクサイ・ニンジン・豆腐のカレー．
- 晩 (2200):
米麦 0.7 合．
ハクサイ・ニンジン・豆腐のカレー．
2005 年 12 月 27 日 (火)
-
0800 起床．
朝飯．
コーヒー．
朝からまた屋久島シュート問題に拘束される．
0945 自宅発．
晴．
1000 研究室着．
-
いきなり 1300 ごろまで，
(A 棟 5F の)
野田さんと群集モデル選択ハナシ．
野田さん学生も質問にやってきて，
師弟だぶる攻撃状態．
確率分割ですかね．
確率分割．
-
じつは群集多様性にしろヘテロ接合度にしろ屋久島多樹種系にしろ，
集団の階層構造がある問題は，
どう考えるべきなのかよく整理できてないところがある
(私のアタマの中で)．
-
昼飯．
-
苫小牧ボスからバグ指摘いただき，
自分のうかつさに脱力する．
教訓．
ファイル名はできごころで変えない．
ゐんどーづ上で動作させるディレクトリシステムに
Unix の symbolic link をまぜない．
-
さらに本日はどういった日なのか
……
苫小牧から村上さんが見える．
「年末に村上さん登場」
というと凶事の別名として使えるぐらいなんだが，
本日は札幌での雑用ついでにこちらに寄られた，
とのこと．
-
亀山さん・工藤さんとツガザクラ雑談少々．
「分布のはし」
においてアオノツガザクラはクローンばかりだった
……
というハナシは種子繁殖ってそれほど重要じゃないかも，
といったことにはならないのかな．
しかしこの植物の場合，
栄養繁殖 (栄養成長)
はかなりゆっくりしてるよーで．
1 年 1cm?
ホントだろうか．
-
北大生協にいって RDRAM うけとる．
512MB 34800 円．
高い．
院生が使ってる研究室 PC 増強用．
研究室にもどってつけかえてみる
(64MB 二本ぬきとる)．
起動してみると合計 768 MB．
問題ナシ，
だったか
(RIMM の PC1066 を PC800 で動かしてる)．
-
……
とどたばたしてて，
なかなかアタマが屋久島問題にきりかわらづ．
時刻はすでに 1730．
ぢつは昨晩なかなか眠れなかったんでそもそもアタマがホワイトアウトぎみで
……
寝ようとしたら，
データ解析がアタマの中で暴走して，
ですね．
-
ぢりぢりと屋久島照葉樹葉寿命の階層ベイズモデル化をススめる
……
-
しかし前途遼遠なまま撤退．
1955 研究室発．
2010 帰宅．
晩飯．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0820):
米麦 0.7 合．
ハクサイ・ニンジン・豆腐のカレー．
- 昼 (1310):
研究室お茶部屋．
米麦 0.7 合．
レトルトパウチド親子丼．
- 晩 (2100):
米麦 0.8 合．
コマツナあえもの．
ネギ・ショウガ・鶏レバ炒めもの．
2005 年 12 月 28 日 (水)
-
0850 起床．
10 時間ぐらい寝た．
おかげでよーやくアタマが修復．
朝飯．
コーヒー．
1010 自宅発．
すげー大雪．
1025 研究室着．
-
屋久島照葉樹の葉寿命モデリング．
年枝を
Segment
としてたのを YearShoot
と改名してコード書きをつづけているんだが
……
ここまでやらなくてもいいかも，
という気がしてきた．
-
なぜこのへんが複雑になりがちかというと
……
-
末端問題:
いちばん根元方向の年枝がどうなっているのかわからない
場合がある
(欠測)
-
休眠・二度伸び問題:
休眠 (根元方向の年枝齢が増加)
・二度伸び (根元方向の年枝齢が減少)
が各芽鱗痕で何回生じたのかわからない
これらをさらに「圧縮」すると
-
いちばん根元方向の年枝の葉痕数がわからない
-
各年枝の
(みかけの齢ではなく)
「ホントの齢」
がわからない
となる．
逆にいえばこの二点以外は観測値として与えられており，
推定計算の最中に「変化」することはない．
-
リスト構造になってる年枝オブジェクト，
というのはある意味かんがえやすいデータ構造なんだけど，
上のような計算のためだけにわざわざ再帰を使いまくる必要あるか?
ということで設計をみなおし中．
-
これはデータを読ませるところから再検討したほうがいいのかな．
いやはや．
-
このあたりで苦闘させられる遠因のひとつは，
いままで葉寿命だのシュート構造だの
研究してきた連中がいかに何も考えずに調査してたか，
というあたりにあるような気がしてならない．
-
よけーな仮定をできるだけ少なくして
-
シュート伸長の休眠だの二度伸びだのの影響を考慮して
葉寿命を推定するにはどうしたらいいか，
といったすごく基本的な問題について，
いままでマトモに考えられていない．
ということで，
私などがここで
「欠測データの処理の方法」
とかいまさらながら検討させられてるんではないかな．
雑学的知識ばかり重視するくせに基本的計算の方法は軽視する連中，
というか．
-
で，
データ読みこみ部分をみなおすと
……
いつものごとく，
こういう階層構造のあるデータをゑくせるで入力するのはよせ，
と言いたくなる．
「集団の中の個体の中の枝の中の年枝の中の葉っぱ」
とゆー階層構造が認識されてないから，
ゑくせるになってしまうんだろうか．
-
昼飯．
外はあいからわず大雪．
-
来月の東京出張は甲山さんからの「依頼出張」になる
……
といった意味不明ぎみの趣旨の自動発信メイルを北大旅費システムからうけとる．
北大の旅費システムにもぐりこんでみたんだけど
(事務処理劣悪な北大にしてはオンライン化されただけでも偉いと評価すべきなんだろうけど，
ぶらうざー依存性などがひどくていまいち感も濃厚)，
状況がよくわからない．
結局，
雪野さんのところに伺うと
……
全てはここに集約されていて，
雪野さんが甲山さん名で依頼出張の手続きをだし，
つぎに久保名でその依頼に対応しておられた．
なンかすごいです．
-
来月始めに「東京三菱銀行」が「三菱東京 UFJ 銀行」になるので，
その変更届を会計掛に提出する．
-
屋久島モデリング，
データよみこみまわりの再整備，
葉数・葉齢分布のデータ構造の簡略化，
階層ベイズモデルの線形予測子の構成，
尤度計算
……
という順で作業をススめていく．
-
さて葉っぱ生き死に尤度計算でいつも注意せんといかんところがある．
ある齢 a の年枝に葉っぱと葉痕があるとする．
葉っぱに関しては齢 a まで生存している，
と言える．
このようにある葉が齢 a まで生残する確率を p(a) とおいてみよう．
いっぽうで，
葉痕は「いつ死んだかわからないけど，
齢 a までには死んだ」という情報だ．
こういう葉痕一個が観察される確率は?
これは簡単で，
1 - p(a) でよいのである．
-
1630 ひととーり書けたように思うので試験運転開始．
現時点の計算プログラムは，
単純に (樹種間共通) + (樹種差) + (個体差)
だけを計算させる階層ベイズモデル版である．
休眠や二度伸びといった常緑樹の呪われ属性は
まだ考慮してない．
こいつらの効果の設計・実装はたいへんそうなんであとまわしだ．
-
あ，
output まわりの関数準備するの忘れてた．
output まわりもめんどくさい．
-
しかし，
やってみたら上部はともかく下部構造は意外と前のやつと重複していた．
Perl のアヤしげなる
@ISA 継承
でもってインターフェイスと実装を引き継ぐべきか?
……
いや，
びみょーにずれてるからダメだな．
-
実装の継承だけやりたい
……
と書いて，
ふと気づいた．
ケガれ言語 Perl の場合，
それは単に「ふつーの関数」
とかに
$self
(べつに
$this
と書いてもよいが)
を渡してしまっていじくってもらえばよい，
ってこと?
いかにもケガれた解決策だな．
なにしろ クラス::メソッド
だってどこからでも呼べるわけで．
呪われ言語 C++ ふうの用語で説明すると，
Perl においては「メンバー変数」がことごとく
``public 属性''
なんで，
それが可能，
ということ．
-
さすがにこれはやりたくない．
-
お茶部屋で休憩．
また院生たちの家庭事情に詳しくなってしまった．
女性大学院生たちの挙動を理解する上で
「パパ」
たちは重要なのかも．
どうなんだろう．
-
事前分布の分散パラメーターは別に出力すべきかな
……
まあ，
しばらく放置．
-
いろいろ手直ししてよーやく MCMC 疾走するようになった．
時刻は 1930．
といっても，
ぢりぢりとしか進捗しないが．
かなりサンプリングの間隔をあけないと
Hyperspecies の階層ではパラメーターが動かない．
いまのところ 20 MCMC step とばしでやっているんだけど，
毎回変わるパラメーターはない．
-
何やら M2 なヒトたちが
「修論おいこみしーづんをがんばろう夕食会」
にでるというので，
修論したうけの私などもそういうトコロに参加可能でしょうか，
とお尋ねしたところ一瞬で却下された．
嗚呼，
身分のいやしいしたうけの悲哀．
-
メトロポリス法で MCMC 動かしてると，
ときどき
「葉っぱの生残確率 10-100」
となってしまうようなトンでもないパラメーターを
「これどうでしょう」
と「提案」してくる．
そういうのは尤度が悪くなりすぎるので，
とうぜん「却下」される．
しかしながら，
ヒドい場合には計算機の数値精度を超えるような値になることもあり
(なにしろ無情報事前分布なんで)
……
その場合，
「生残確率ゼロ」
にされてしまうこともある．
こうなると対数尤度が評価できなくなるんで，
ごまかし補正をいれてみる．
-
この小細工は最終的に，
確率 200 万分の 1 の「上下げた」をはいた logistic 関数，
という奇妙な関数型にまで発展した．
私の憶測なんだけど，
この方式は世界中の MCMC 使いたちのあいだでで
こっそりと愛用されてるにちがいない
……
いや，
ぎぶすサンプラーをちゃんと書けば，
かかる面倒は不要なのか?
-
パラメーターひとつひとつのメトロポリスサンプリングごとに
すべての個体のすべての枝のすべての年枝のすべての葉っぱでの
尤度を評価しなければならないので計算がおそい．
-
いまいちうまく計算できてないような気がするまま撤退．
2105 研究室発．
2120 帰宅．
晩飯．
-
試験運転，
真夜中すぎまでつづく．
というか，
計算機には徹夜計算を命じて寝る．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0920):
バターロール．
- 昼 (1300):
研究室お茶部屋．
米麦 0.5 合．
コマツナあえもの．
ネギ・ショウガ・鶏レバ炒めもの．
- 晩 (2150):
米麦 0.8 合．
コマツナあえもの．
ネギ・ショウガ・鶏レバ炒めもの．
2005 年 12 月 29 日 (木)
-
0800 起床．
朝飯．
コーヒー．
-
昨晩からやらせている計算，
8 時間ほど経過したけど 900 step (MCMC step 数でいうと 4500)
しか進捗してない．
いったん停止して，
「とりあえず作図」
してみる．
-
うーむ，
妥当そうに見える部分と，
ちょっと「?」な部分がある．
-
定数部分はまあ (いまのところは) どうでもいい．
齢依存のところは妥当にみえる．
こいつだけ「個体差」が重要になっている．
これはおそらく「休眠」とかで説明できるにちがいない
(現時点の推定計算では休眠の効果は考慮してないんで)．
-
明るさ依存もまあこんなところか．
おそらく休眠とか入れると明るさ依存性がわかりやすくなるかも．
暗いところで眠りこんでしまう連中なんで．
-
わかりにくいのは窒素依存性の明瞭さ，
だ．
私はここでいう「窒素 (濃度)」
(当年葉内の重量ぱーせんと)
などという測定値は葉っぱ死亡確率の説明に何の役にもたたないだろう，
と予測してたんだけど
……
この樹種差を超えたそろいかたはどういう意味なんだろう?
個体差がめちゃくちゃになってるのは，
測定値が信用できないことのあらわれ，
という気もするけど．
-
しかもおもしろいことに，
「窒素濃度が高いほど長生き」
という方向で．
私はそもそも葉っぱぎょーかいのヒトではないんで，
いんちき解析からひねくりだされた
「窒素濃度が高いほど短命」
といった俗説
(ぎょーかいの中のヒトたちによると
「だれもが認める真実」)
なんぞは気にしてないんだけど
……
しかし「長生き」な方向に樹種間でそろってしまうとは，
ねえ．
何があったんだろうか?
単純な計算まちがいとかだったら気楽なんだけど．
-
とりあえず
1100 自宅発．
曇．
1115 研究室着．
-
もとデータにたちかえって検討．
-
よくわからんな．
樹種によっては，
明るい暗いと窒素のあいだに交絡ありそうな気はする．
それから牛原さん定義するところの span
は寿命の長短と必ずしも関係ないからなあ
(むしろ個体差，
というべきかも)
……
-
さらに以前つくってみた葉寿命分布図ながめてると，
たしかに「明るさと呼んでるものは個体差」なるあつかいかも生じうるかも，
という気はするな．
個体差ナシモデルのほうが，
やや強引なあてはめだった
(「明るさに依存してる」は言いすぎだった)，
ということか?
-
窒素に関してはよくわからんので試験的に，
窒素を
unused
に指定して推定計算やらせてみる．
-
と同時に，
末端欠測問題なんかもうまくあつかえていないような気がするので，
このあたりの改善に着手．
-
いろいろ調べてるうちに，
計算プログラムのばぐは見つかった．
なんとまあ個体差の prior から樹種差を生成させていた．
どちらも平均ゼロのガウス分布とはいえ，
いやはやなまちがいだ．
しかしこれを修正しても結果はあまり変わらないと思うんだよね
……
-
じつはこのばぐは以前のシュート伸長休眠・二度伸びの推定計算にも
影響およぼしているんで，
まずはそちらから計算やりなおし
……
こちらは計算が多少は速い．
そして様子をみていると，
幸か不幸か，
あまり変わらないように思える．
ただし Hyperspecies の階層にはおいては，
ということだけど．
-
ということで，
葉寿命推定計算のほう，
まだ何かばぐが残存しているにちがいない．
-
昼飯．
院生たちはまだいる．
-
コードみなおし & ちょい試験運転，
つづく．
-
ばぐは見つからないけど，
継承 (inheritance)
濫用したりすとらはおそるべき勢いにて進行中．
つまりばぐも継承する，
と．
しかしながら，
ばぐがあちこちに「こぴー」されるよりは 100 倍マシ．
-
しかし，
この教科書的とでもいうほど素直な「差分プログラミング」，
親クラスの挙動がアタマに入ってないと
「えーっとこいつのコンストラクターって何やってたっけ?
……
あ，
親クラスはじーさんクラスのコンストラクター呼んでる」
といった迷宮放浪になってしまいそうな．
-
むろんどんどん先祖たどりをやれば，
すぐに終点にいきつくわけだが．
たとえば
TreeL
だと
OwnerParameter
<- OwnerFixedEffects
<- TreeP
<- TreeL
で終了．
これが一番長い．
-
今日から全館暖房は止まってるんで寒い．
ガスストーヴつける．
-
シュート伸長の休眠・二度伸びの再計算，
105 分で終了
(100 step とばしの 1500 sampling)．
幸か不幸かほとんどかわらなかった
(以前の図)．
hyperprior は (そこそこの値であるなら)
なんであってもそれほど影響が大きくない，
ということなのか．
-
とゆーことで，
葉寿命推定計算のばぐらしきものはいまだ特定できづ．
アタマをひねってても進展しないので，
A 801 室の Dell 機では試験運転を実施し，
手もとの ThinkPad ではコードみなおしのつづきを．
尤度計算まわりもアヤしいかもな．
-
ばぐわからん．
hyperprior ばぐ修復の功名なのか，
試験運転中の葉っぱ生き死にの推定 MCMC 計算は以前よりマシなかんぢ
の値を吐きだしつつある．
明るい・暗いの効果がかなりはっきりしている
(明るいと短命)．
しかしデータとの対応がよくわからぬ窒素長命化効果とやらも
(今朝えられたものよりマシとはいえ)
まだ残存してるんだよなぁ
……
明るさとの交絡 (confounding) が問題なのかな．
-
重量あたりの窒素量だの面積あたりの窒素量だの，
なンか気分の悪くなる割り算値ばかりあつかってるわけだが
……
とりあえず図にしてみる．
ぜんぜん違ってくるわけで．
呪われ割り算値．
-
一番単純なハナシとしては
……
窒素量なるものは単なるノイズ，
とでも言ってしまうことか．
こういう推定プログラムはゴミのような情報に
存在するゴミ相関を悪用したりもできるんで．
-
コードながめててもばぐらしきものがみあたらない．
あるいはデータの問題なのかもしれないけど，
よくわからない．
MCMC 計算そのものがまずいのかもしれない
……
どこかにある尤度の小岩峰で「引っかかっている」
(「登ったはいいけど降りられなくなっている」?)
とか．
-
で，
このあたりを解明するためには，
計算そのものをよく観察し，
挙動を把握する実験を試みる必要あるんだけど，
計算がおそくてなかなか進捗しない．
これは，
昨日もかいたとーり
「パラメーターひとつひとつのメトロポリスサンプリングごとに
すべての個体のすべての枝のすべての年枝のすべての葉っぱでの
尤度を評価しなければならない」
てなことを Perl なんかでやらせているためである．
-
とゆーことで，
SWIG
で尤度計算の部分だけを C コードに置き換えること検討中．
これは Perl その他から C/C++ 関数を
「お手軽に」
よびだせるしくみ．
いままで使ったことないんだが．
Vine Linux だと
swig-1.3.21
の RPM pacakge が準備されてるな．
-
1850 研究室発．
1900 帰宅．
晩飯．
-
SWIG
の勉強ぢりぢりと．
便利なのか不便なのかよくわからんシステムだな．
「配列を渡す」
というところはかなりとりっきーだ．
このあたり，
また明日以降に解説の予定．
-
[今日の運動]
-
[今日の食卓]
- 朝 (0850):
米麦 0.7 合．
コマツナあえもの．
納豆．
ダイコン・ニンジン・ネギ・ナメコ・タラすりみ味噌汁．
- 昼 (1340):
研究室お茶部屋．
米麦 0.5 合．
ダイコン・ニンジン・ネギ・ナメコ・タラすりみ味噌汁．
- 晩 (1930):
米麦 0.8 合．
ダイコン・ニンジン・ネギ・ナメコ・タラすりみ味噌汁．
2005 年 12 月 30 日 (金)
-
0740 起床．
朝飯．
コーヒー．
-
屋久島常緑樹葉生残過程モデリング，
昨日のつづき．
Perl による数値計算はおそいので，
葉っぱ生き死にの尤度だけを C の関数でやらせよう，
という方向にむかいつつある．
Perl から C 関数の呼出しに
SWIG
(simplified wrapper and interface generator)
使ったらいいのでは，
と考えて勉強中．
-
たとえばこの
Tutorial
にある例にならって
calc.c
なるものを作ってみる．
double z = 3.0; /* ぐろーばる変数 */
double calc(double x)
{
int i;
double y;
i = 2;
y = i * z + x;
return y;
}
これに対応する SWIG interface である
calc.i
をこう定義する．
%module calc
extern double z;
extern double calc(double x);
こいつらをコンパイルする Makefile
を書いてみる．
TARGET=calc
CORE=/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/
all:
swig -perl5 -module ${TARGET} ${TARGET}.i
gcc -fpic -c ${TARGET}.c -I${CORE}
gcc -fpic -c ${TARGET}_wrap.c -I${CORE}
gcc -shared ${TARGET}.o ${TARGET}_wrap.o -o ${TARGET}.so
しかし，
make
するとエラーがでる．
gcc -c calc_wrap.c -I/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/
In file included from /usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/op.h:484,
from /usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/perl.h:2344,
from calc_wrap.c:291:
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/reentr.h:611: error: field `_crypt_struct' has incomplete type
In file included from /usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/perl.h:3552,
from calc_wrap.c:291:
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:199: error: 文法エラー before "off64_t"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:201: error: 文法エラー before "Perl_do_sysseek"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:201: error: 文法エラー before "off64_t"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:201: 警告: data definition has no type or storage class
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:202: error: 文法エラー before "Perl_do_tell"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:202: 警告: data definition has no type or storage class
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:1308: error: 文法エラー before "Perl_PerlIO_tell"
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:1308: 警告: data definition has no type or storage class
/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/proto.h:1309: error: 文法エラー before "off64_t"
ということで，
計算機の訴えを注意ぶかく検討し，
自分の置かれた環境にあわせた対応が必要になる．
いろいろ調べてみてわかったんだけど，
いちばん安直な解決策としては，
各自の環境において
perl -e 'use Config; print $Config{ccflags};'
が生成する Perl まわり用コンパイラーオプションを
そのまま流用せよ，
というもの．
Makefile
をこう書き直すと無事にコンパイルできる．
calc.pm
や
calc.so
が生成される
(dir)．
TARGET=calc
CORE=/usr/lib/perl5/5.8.2/i686-linux-thread-multi/CORE/
OPT=-D_REENTRANT -D_GNU_SOURCE -DTHREADS_HAVE_PIDS \
-fno-strict-aliasing -D_LARGEFILE_SOURCE \
-D_FILE_OFFSET_BITS=64 -I/usr/include/gdbm
all:
swig -perl5 -module ${TARGET} ${TARGET}.c
gcc -fpic -c ${TARGET}.c -I${CORE} ${OPT}
gcc -fpic -c ${TARGET}_wrap.c -I${CORE} ${OPT}
gcc -shared ${TARGET}.o ${TARGET}_wrap.o -o ${TARGET}.so
このようにして自作した
calc package
を呼びだす Perl スクリプト例はかくのごとし．
まあ，
さすがに簡単さをウリにするだけはある．
use calc;
print &calc::calc(3.0)."\n"; # 9.0 と表示される
$calc::z = 9.0; # ぐろーばる変数をいじる
print &calc::calc(3.0)."\n"; # 21.0 と表示される
面倒なのは calc::calc()
にスカラー以外のオブジェクトを渡したい場合だ．
たとえば配列を引数にしたい，
という問題を考えよう．
SWIG の機能をもう少し使いこなす必要がある．
いろいろな解決策がありそうなんだけど，
ここではヘルパー関数による方法を試みてみよう．
SWIG 本家サイト
にある解説文書
Perl Extension Building with SWIG
の ``3.4 Helper Functions''
まわりそのままである．
calc.c
をこう変えてみる．
double z = 3.0;
double calc(double x, int imax, double* p)
{
int i;
double y;
y = z + x;
if (imax > 0) {
for (i = 0; i <= imax; i++) {
y += p[i];
}
}
return y;
}
対応する SWIG interface
calc.i
に double 配列の
(いわば)
constructor / destructor / setter / getter
を定義するわけだ．
/* calc.i */
%module calc
%inline %{
double* double_array(int size) {
return (double *) malloc(sizeof(double) * size);
}
void double_destroy(double *a) {
free(a);
}
void double_set(double *a, int i, double val) {
a[i] = val;
}
double double_get(double *a, int i) {
return a[i];
}
%}
extern double z;
extern double calc(double x, int imax, double* p);
Makefile
は上記と同様でよい．
こいつの Perl ドライヴァーはこういうかんぢで．
use calc;
sub create_array
{
my $len = scalar(@_);
my $ia = calc::double_array($len);
for my $i (0 .. ($len - 1)) {
my $val = shift;
calc::double_set($ia, $i, $val);
}
return $ia;
}
$calc::z = 9.0;
my $imax = 3;
my $p = create_array(1, 2, 3, 4);
print calc::calc(3.0, $imax, $p)."\n"; # 22 と表示される
calc::double_destroy($p)
ついでに $p
はどういう reference になってるんだ，
と表示させてみると _p_double=SCALAR(0x804c464)
となってた．
C ポインターへの reference?
とりあえず配列を渡すことができれば，
葉っぱ生き死に尤度は計算できそうだ．
SWIG はなかなかこったつくりになっていて，
さらに複雑なデータ構造のやりとりもできるのもウリらしい
(けど今回は使わない)．
1040 自宅発．
曇．
1055 研究室着．
昨晩命じておいた葉っぱ生き死に階層ベイズモデルの MCMC 計算，
10 時間 40 分かかって終了していた
(これはまだ all Perl)．
100 MCMC step とばしの
1500 sample の最初の 500 を捨てて作図．
収束はすごくおそい．
500 ではだめかも．
-
このモデルでは休眠などの効果をまだ考慮してないんで，
葉齢は「みかけの葉齢」である．
で，
こいつは「みかけ」であるためなのか，
ばらついている，
と．
-
明るさ依存性は「明るいと短命」な方向だけど，
よくわからん，
と．
樹種差・個体差がそろってしまっているのはなんか意味あるんだろな．
-
(葉面積あたりのとゆー)
窒素割り算値依存性は
……
妥当というかナゾというか．
全体としては効果ゼロ，
そのばらつきはおそろしく小さい．
そして樹種差・個体差が激しい
……
ように見える．
樹種差・個体差はどちらかといえば，
「窒素長命化」
な方向で．
なンかの交絡なのか?
-
よくわからん．
そもそもこの計算だってアテになるのかアヤしいもんだ．
とはいえ改善案もないので，
窒素濃度を重量ベイスにして計算やりなおし命じる．
無策だ．
-
SWIG まわりのことなどうだうだ書く．
塩寺さんに熱帯林葉っぱデータの準備できたと言われてあせる．
しばらくお待ちください，
とゆーことで．
-
昼飯．
今日は院生密度がさすがに低い．
-
1/2 にいきなり京都往復，
というよーなハナシが．
うーむ．
とりあえずいろいろ調べてみる．
-
のんびりと修論にとりくむ石橋君のミズキ
モデリングこんさる．
個体差が見えにくい．
-
歳末も遺伝子をよむ亀山さんとツガザクラ雑談．
この植物，
クローナル成長が重要なのか重要でないのか，
まったくわからなくなった．
場所とりは何で決まるか
……
さらにこれと自殖の可否との関係もよくわからん．
たとえば，
単一クローンにみえるのが全部自殖種子由来，
とかありえんのかしらん
(親がヘテロなんだけどどちらの染色体も劣性致死になってる，
とか)．
自殖率と近交係数から近交弱勢が決まる，
と教えてもらったんだけど，
そうなると近交係数の空間構造依存性が気になってきたり
……
このあたり，
あまり重要とは考えられてないよーな．
-
SWIG
による屋久島モデリングの尤度計算高速化，
つまり Perl を部分的に C に置き換える，
という試み．
やってみたら，
100 MCMC step 計算に要する時間が
56.2 秒
から
36.3 秒
に短縮された
……
うーむ，
改善と言えるかどうかびみょーな気も．
もうちょい速くなってほしい．
-
尤度計算だけでなくメトロポリス法とかでもてこずっているのかね
……
よくわからんな．
ぎぶすさんぷらーだと速くなる，
とか?
どうなんだろう．
-
ぎぶすさんぷらー化するのはたいへんなので，
安易なトコロから
……
正規乱数生成とか事前分布まわりの尤度計算も C に置き換えてみる．
100 MCMC step 計算に要する時間が
36.3 秒
から
35.3 秒
に短縮された
……
うーむ，
努力がむくわれん．
Perl もそこそこに速い，
ということかな
(いちおーコンパイルしてるそーで)．
-
まあ，
これ以上速くしようとするなら，
コード全部を呪われ言語 C++ で書き直すしかあるまい．
その時間的よゆーがなく，
無理っぽいような気がする．
-
ぢたばたしても改善されないので，
撤退準備
……
すると，
なぜか晦日も
R
作図にはげむ西村さんから
「
par(mfrow = c(縦, 横))
で図を並べるとつぶれてしまう」
質問．
これは FAQ というか
……
(grid を用いない場合に)
各わく内でつぶさないポイントは
par(mar = c(下, 左, 上, 右), mgp = c(軸名, 軸数字, 軸線))
の調節 (特に 上, 右)，
そして出力 device のサイズ調節である．
-
2210 研究室発．
買いもの．
食料の多くを輸入にたよっている北海道では，
冬の荒天による輸送とどこおりに敏感である
……
野菜の値段とか，
が．
たとえば，
ホウレンソウ一把 238 円．
私は 1-2 食ぶん 150 円を超える野菜はなかなか買えない．
2230 帰宅．
晩飯．
-
窒素の重量ぱーせんととやらを使ったモデルの計算結果がでていた．
501 から 1500 番のサンプルを使って様子をみてみると，
昨日朝に「とりあえず」作図したのとあまり変わらないような気がする．
-
おそらく窒素まわりが明るさと交絡してるンでは，
という予感．
シュート伸長の休眠とかをいれると，
何か変わるんではないかな．
現時点の予想としては，
「明るさ」の効果が明確になり (暗いところでよく休眠するので)，
窒素うんぬんはそれによって抹殺される
……
となってればハナシは簡単なんだけどな．
-
とはいえ，
このあたりの設計・実装にはもうちょい時間かかりそう．
明日ぐらいには何とかなるか?
ともあれ，
現在うごかしている MCMC 計算プログラムの信頼性には，
まだまだ疑うべき点が残されているように思うので，
昨晩から本日にかけての計算の検算 (異なる初期状態からの MCMC 計算)
を A 棟 8F の Dell 機に命じておく．
-
[今日の運動]
-
うーむ，
このうんどう不足状態はどうやったら解消されるのか
……
ぢみちに腕立伏せとか?
やらないよりはよほどマシだろうな
(肩こりにも効きそうだし)
……
-
[今日の食卓]
- 朝 (0820):
米麦 0.7 合．
ダイコン・ニンジン・ネギ・ナメコ・タラすりみ味噌汁．
- 昼 (1330):
研究室お茶部屋．
こんびに握り飯二個．
- 晩 (2250):
米麦 0.7 合．
ダイコン・ニンジン・ネギ・ナメコ・タラすりみ味噌汁．
2005 年 12 月 31 日 (土)
-
0800 起床．
朝飯．
コーヒー．
洗濯．
怠業
……
-
……
ぎみなんだけど，
よゆーがないのでネットから北大にもぐりこんだり，
昨日の計算結果とかを調べ直したり．
-
通販で発注してたパンツ (ズボン) 三本とどく．
通販といってもゆにくろですが．
いやはや，
数すくない手もちがことごとく穴あき状態なんで．
スソとかに穴あいてるのは機能上はさほど不便はないんだけど，
ヒトがみると
「何か特別な主義主張にもとづいて穴あきのままなのか?」
と誤解されがちなんで
(あるいは大学院生たちから
「ちゃんとしてくださいよっ」
と譴責されるんで)，
今回の発注におよんだ次第．
-
素材は綿でも悪くはないけど，
ひとつぐらいはウールにしとけばよかった
……
あ，
いま思いだしたんだけど，
ウールパンツというとついつい
「野外でも使える良いものを」
とゆー厳しい規準で選ぼうとしてしまい，
ゆにくろにはそれがなかったんだよね．
-
三本の中の一番の安物であるポリエステルものをはいて
1250 自宅発．
晴れてるけど雪ふり．
化繊でも意外と感触は悪くないかな．
ちょい冷えるような気がするけど．
1305 研究室着．
-
さて，
屋久島常緑樹の葉っぱ生き死に問題，
おそらく最後に残された難問である
(最後であってほしい，
最後ぢゃなきゃイヤだ)
シュート伸長の休眠 &
二度伸びの効果の具体的なモデリングについて朝から検討している．
-
なんでこんなものなんぞを考慮するのかといえば，
ぢつは葉齢なる観測値は「芽鱗痕」を利用して憶測されるからだ．
ある枝において芽鱗痕が毎年毎年ひとつづつ増えていくなら，
何も問題なく葉齢は特定できる．
しかしながら，
年によって芽鱗痕がまったく増えなかったり，
ふたつも増えたりすると「みかけの葉齢」と「ホントの葉齢」
は食いちがってくるのである．
-
これまでの常緑樹の葉齢研究では，
こういう「みかけ」と「ホント」のずれみたいな
都合の悪いことは「なかった」ことにされている
(もちろん根性だしてマークつけ
& 5 年以上もの追跡やるとゆー戦略もあるわけだが)
……
そのくせ
(別の文脈では)
シュート伸長の休眠があーだこーだと言いたがるのは不思議なコトなんだけど．
-
伸長の休眠と二度伸びに関する hyperspecies - species - tree
の階層ある階層ベイズモデルの推定はすでに終了している．
おそらく一番首尾一貫してんのは，
こういう休眠推定と葉死亡確率推定を同時にやることだろう．
しかしながら，
これはたいへんだと思うんで，
今回は見送り．
-
で，
休眠まわりと葉死亡確率を分離する，
という方針でいくなら
……
うーん，
樹種ごとの事後分布の平均値を使う，
という方向でいくのが無難そう．
個体差はどこにいったんだというハナシもあるが，
これはおそらく
「休眠とかの個体差」と「葉っぱ死亡の個体差」
をこの観測データからは識別不可能だろうから，
葉っぱ死亡確率の個体差でもって代用すればよい，
と．
-
さて，
これをどう MCMC 計算に組み込めばよいか
……
と考えてみたんだけど，
いやはや
またしても Gibbs sampler が導出できないので，
とりあえず Metropolis-Hastings 法つまりメトロポリスサンプラーで．
-
いま，
「ある芽のある一年のふるまい」として
(「確率」は事後分布の平均値である)
……
-
二度伸びする
(根元方向のみかけの葉齢に与える影響 d = -1)
: 確率 p
-
ふつーに成長 (d = 0)
: 確率 1 - p - q
-
休眠 (d = 1)
: 確率 q
がわかっていたとしよう．
休眠・二度伸び現象の発生は毎年独立であるとする．
-
今度は逆に
「ある芽鱗痕を選んだすると，それはどういう状態にあるのか」
を考えてみよう
……
こういう確率分布で与えられる:
-
二度伸びでできた (d = -1)
: 確率 p
(二度伸びの場合は d = -1 しかありえない)
-
ふつーの成長でできた (d = 0)
: 確率 1 - p - q
-
k 回休眠している (d = k > 0)
: 確率 qk * (1 - q)
(k = 1 ..∞の和をとると q になる)
これらを休眠・二度伸び現象に関する尤度計算につかえばよい
……
これは事前分布
(prior distribution)
ではない
(この部分に関しては事前分布なし?)．
たぶん．
-
さて，
よくわからないところが各芽鱗痕の状態遷移だ．
愚直にやるなら，
各 step において
全芽鱗痕ひとつひとつについて
(d をランダムに ±1 増減させて)
メトロポリスサンプリングすることになる
……
これは計算時間から考えて不可能．
-
次におもいつくのが，
いくつかの芽鱗痕をもっている「枝ごとに」めとろぽりすする，
という策なんだが
……
このときに，
新しく「提案」するパターンを生成するときに
全芽鱗痕で d を ±1 変動させたらほとんど
「却下」
されるような気がする
(やってみないとわからないけど)．
-
そこでこれを改善する妥協策として，
ある枝を選び
(枝といっても長さはさまざまなんだけど)，
その中のどれかひとつの芽鱗痕での d を ±1 変動させたものを
「提案」してメトロポロス殿に審査していただく，
という計算方法がありえるんではないかな，
と考えている．
-
これのまずそうなところは枝の長短
(正確には枝がもつ芽鱗痕数の大小)
になンか影響されそうなあたり，
かなあ
……
まあ，
こんなんで許されるんではなかろーか．
ひとまず昼飯．
-
お，
Dell 機にやらせてた計算が終了している．
昨日の SWIG による高速化で
計算時間は 10 時間 40 分から 7 時間 30 分に「短縮」されてるわけだが
……
あまり短縮感ないなぁ．
-
図の保管場所としてここが便利なんで，
またしても陳列しておく
(501-1500 番目のサンプルのみを使った事後分布図)
……
毎日毎日同じよーな図ばかりを，
と思われるかもしれないが，
これは昨日とまったく同じ条件
(ただし今回は SWIG 使ってるけど)
なのに，
「だいたい同じだけど，びみょーに違う」
というところが気になっている．
観測データの交絡が原因なのか，
MCMC サンプリングのやりかたに問題があるのか，
それとも尤度まわりの定式化と数値計算にばぐが残存しているのか．
-
上の図で，
Hyperspecies 階層の
Narea
なんてのは，
おそらく MCMC 計算がうまくいってなくて
(0.4 あたりに「ひっかかっている」ように見え，
これは明らかに hyperspecies のパラメーター値として適切ではない)，
その原因は
……
観測データとモデルの両方?
つまり「よくわからん」ということです．
-
よくわからんまま，
Hyperspecies の層におけるパラメーターの遷移
(Narea を使ったモデル)
の図を置いとく．
-
上の図をみると，
step 200 あたりですごい「大ジャンプ」が生じてるよな
……
これって MCMC 計算もうちょい工夫しろってことかしらん．
パラレルテンパリング法 (拡張アンサンブル法)，
とか．
計算がますます遅くなるな．
-
パラレルテンパリングはしんどいので，
simulated annealing 法もどき，
とか?
最初は「温度」を高めにしておいて，
徐々に下げていく，
とか．
いや，
これではあまり改善されんような気がするな．
-
「ごくまれに」
M-H ルールを無視させてみる?
アヤしげではあるけど
……
-
ついでに，
全体の fixed effects の事前分布の分散逆数の遷移
(Narea を使ったモデルの 501-1500 番目のサンプルのみ)
の図も．
ふーむ
……
