ぎょーむ日誌 2006-06-(01-10)

2006 年 06 月 01 日 (木)

• 0750 起床． 朝飯． コーヒー． 朝からアカマツ資源分配モデルの検討に没頭してしまう． 0915 自宅発． 晴． 0930 研究室着．
• メイル書きとかネット雑用とか．
• GLMM 解説記事のしあげについてはもうちょいコメント待ちすることに． アカマツデータながめる．
• 昼飯．
• アカマツモデリングのやり直しの可能性を検討する データみなおしつづく．
• この樹木個体成長モデルは (ぱいぷ樹木 などとは異なり) たいへん単純化されたものになっていて，
  1. 各軸が直径成長 (age 依存)
  2. 枯れ上り (断面積差依存)
  3. 新しい軸の生成 (random resampling)
  でハナシがおわる．
• で， いま検討しなおしてるのは，
  1. 各軸が直径成長 (age 依存) …… age 依存ってのもなぁ
  2. 枯れ上り (断面積差依存) …… これは面白いんだけどちょっと強引かも
  3. 新しい軸の生成 (random resampling) …… これは他に方法がないので， この方式でいくしかない， ただし「どれだけ作るか」には検討の余地がある
  といったあたりだ． 1. と 2. をまとめて変更できないかな， と．
• そしてこれらの過程で 3 つの未知パラメーターの推定値が 必要になる:
  • 軸のふとり速度
  • 枯れ上りの速度
  • 「どれだけ作るか」 …… 当年における膨張
  こいつらのバランスでバイオマス蓄積がきまるんだよねえ． モデルが変わると当然これらの推定方法と推定値が変わるわけで (あるいは必要とされるパラメーターが変わる)．
• につまってしまったので， 論文読み． 来訪中の Caswell さんとその PD だった Fujiwara さんの クジラ生残推定論文 をススめられたのでダウンロードしてみる …… 「見つかった・見つからなかった」データから生残確率の最尤推定とモデル選択， か． ここで「このクジラは観測場所 X 付近を好む・好まない」 という「個性・個体差」を個体ごとの割算推定値でくみこんでしまっているな． ここは何かしら mixed model あるいは Bayes model 的に考えるべきところだ．
• ともかく行列モデルは推定すべきパラメーターが多すぎる．
• 1800 すぎから地環研前芝生でジンパ (ジンギスカン用肉とよばれるラム肉をあまり調理むけではない使い捨て アルミ鍋にいれて七輪の上で加熱してみる儀式もしくはパーティー)． なんだかひたすらウチワでぱたぱたと炭火に送風していたような気がする．
• 2000 すぎおひらき． 研究室にもどってまたアカマツデータをひねくる． 2210 研究室発． 2225 帰宅．
• [今日の運動]
  • また運動休養日．
• [今日の食卓]
  • 朝 (0830): 米麦 0.6 合． コマツナあえもの．
  • 昼 (1250): 研究室お茶部屋． 米麦 0.6 合． コマツナあえもの．
  • 晩 (1830): ぢんぱ．

2006 年 06 月 02 日 (金)

• 0740 起床． 朝飯． コーヒー． またアカマツ問題の検討． 起床直後はアタマの中が整理されているような気がする． そして整理されたアタマだからこそいささか正気から逸脱してる策がだせたりする …… よーな気がする． 0930 自宅発． 曇． 0945 研究室着．
• 午前中はアカマツ検討で終わった． まずは昨晩つかった七輪などかたづけつつ， アタマの中であれこれと． しかし七輪の底板に修理が必要な気もするな．
• 有用かもしれないけどまだ使ってないわざを整理してみる．
  • 加重つき resampling
  • EM アルゴリズム使った 「可能な cohort 構成」 集合のもとでの最尤推定
  • MCMC 計算による cohort 間の接続の集合の生成 (Gibbs assignment?)
  などなどいろいろあるよな …… ということで 「もしこいつらが有効に使えるとするならば」 とホワイトボードにらくがきしつつ検討．

[ホワイトボードらくがき]

今朝から「問題をどこまで簡単化できるか」を考えている． この段階ですでに 4 つある stem cohort のうち， 新しい二つ (cohort1999 と cohort2000) の関係だけを検討するモデルに単純化している．

• 飯つかいまわしローテイションがずれたので， 北大生協まで昼飯調達の旅． 歩きながら考えるのはなかなか快適である …… 少なくとも煮つまり感はウスい． 昼飯かってお茶部屋にもどると， 「なンで生協なんかでかってくるんですか」 と院生にまた叱られてしまった． 昨日から日曜日まで北大大学祭で， そこらへんの屋台 (特に留学生がお国料理を売ってるところ) で何か良いもの (おそらく院生のぶんもふくめて，という意味だろう) 買ってきなさいよとの趣旨であった． 叱られつつ昼飯．
• 今度は紙の上にモデルの概念図をかき， いきづまると A 棟 8F 内をうろうろする．

[紙の上のらくがき]

単純化はさらにつづく． モデルから将来予測の能力を剥奪すると， さらに推定計算が簡単になる …… と思いつつ図を描いてながめてるうちに ……

• 紙の上に stem 成長の概念図を描き， その推定計算用モデリングの定式化をかきこんでいると， また単純化に気づいた． cohort1999 と cohort2000 を接続する Gibbs assignment は興味ぶかくやってみたい計算なんだけど， ぢつはこれすらもいらん， ということに気づいた． cohort2000 の生産力は cohort1999 の「支持力」に依存しているからだ．
• ここまで単純化すると Bayes モデルの推定計算はここまで簡単になる:
  • stem set からの sampling ~ (個体ごとの stem set)
  • 「処理」個体の場合はサイズの補正
  • cohort1999 の 2000 年成長期における成長量 ~ (ばらつき指定の事前分布)
  • ばらつきパラメーター ~ (無情報な超事前分布)
  • 成長速度パラメーター全個体共通部分 ~ (いきなり無情報事前分布)
  • 成長速度パラメーター「処理」の効果 ~ (ここも無情報事前分布)
  • 成長速度パラメーター「個体差」 ~ (ばらつき指定の事前分布)
  • 「個体差」のばらつきパラメーター ~ (無情報な超事前分布)
  これで MCMC 計算を駆動すればよい． すごく簡単だ …… あの 屋久島多樹種系 とかに比べたら．
• さて， 計算の実装はどうしたもんかな …… と．
• 帰宅 30 分ばかり前ごろから， お茶部屋からしきりにがやがやわーわーと …… と思っていたら， やはり予想どーりというか大学院生たちをここまで興奮させるハナシは といえば食べものに関する議論であった． 1935 研究室発． ふらふらと歩いてるとみょーな考えが． このアカマツ Bayesian 推定計算， じつは MCMC (Markov Chain Monte Carlo) 法ではなく MCEM (Monte Carlo Expectation-Maximization) 法で計算できるんぢゃないかしらん， という気がしてきた …… えーとこれって正規分布ばかりなんだけどふつーの 一般化線形ではないから nlme()? いや待て， mean = exp(……) * sなる関数形だから exp(…… + offset(log(s))) として family = gaussian(link = "log") の lmer() か? …… といった方策を検討してみたんだけど， 対数尤度の期待値の最大化というのも迂遠な気がしてきたので (いや， 生成した MCMC sample を流し込んで計算すればいいのか?) まあとりあえずは MCMC かな， ということで． 2000 帰宅． 運動． 体重 72.8kg． 晩飯．
• [今日の運動]
  • エアロバイク 60 分間．
• [今日の食卓]
  • 朝 (0810): 米麦 0.6 合． 海藻スープ． コマツナあえもの．
  • 昼 (1340): 研究室お茶部屋． 北大生協のりまき．
  • 晩 (2200): 米麦 1.0 合． ゆでサケ． ハクサイ・ニラ・シイタケの炒めもの．

2006 年 06 月 03 日 (土)

• 0740 起床． 朝飯． コーヒー． 洗濯たくさん． データ解析こんさるメイル & R プログラム書き． 怠業．
• あ． いきなり親知らず (左上) の一部が「欠落」してしまった． 痛みとかは今のところない． 欠落した破片をつくづくとながめてみるに …… どうも欠けた原因は虫歯で (自覚症状なし)， しかも歯石とか成層してるし． いやはや．
• しょうがないからひさしぶりに歯医者でも行くか． この親知らず， ちゃんとはえてないからたぶん抜歯だろうな …… ということで， 自宅近傍の歯科を Google maps と デンターネット で検索する．
• で， 何件か電話してみたんだけど， 土曜日は混んでるね． そもそも休診日とか午前中だけのところもおおいし． ということで， 本日みてもらうことはあきらめて， 自宅から歩いて 1 分のところに電話予約． 月曜日 1030． 半休とりますか．
• 昼飯くってもくだんの歯まわりで特に問題ないのは， 不幸中のさいわいというべきか． 1450 自宅発． 晴． 北 9 東 5 ほ Homac などふらふらして 1620 研究室着．
• さて， アカマツの件はいったんおいて， まずは生態学会和文誌の GLMM 解説記事原稿を 月曜日に提出できるようにせんといかん． この解説記事の共著者の粕谷さん， 何やらいそぎの調査か何かで五島列島の南西の無人の 男女群島 (地図 で見るとまさに絶海の孤島だ; 写真) に行ってしまわれたので数日前から通信途絶状態． 原稿に関する最終コメントはもらえそうにない．
• ということで， まああとは当方で直せるだけ直して提出しますか， と作業に着手．
• というか， 本日はおもに体裁の修正だけでオワった． 原稿依頼時点では提出原稿の体裁については無指定だったので， てきとーでいいのかと思っていたんだが， 先日とどいた 和文誌最新号みると投稿原稿の体裁の規定みたいなのがあったので， それに準じた変更をしてみる．
• つまり何がやりたいのかというと …… 原稿を書いてた段階ではこういう体裁， つまり和文誌の刷り上がりにちかいかんぢで出力していた． 読者が目にするであろうレイアウトに近い状態で推敲できるのは便利だ．

• これをまあ下のごとく (投稿規定にしたがって) わざわざ読みにくい体裁に変更する作業である．

• こういうときに TeX & LaTeX は便利なのである． 文章の内容と見てくれを分離 できているので．
• まず体裁ファイルを二種類作る．
  • 見てくれのよいレイアウト (上の画像)
  • 見てくれの悪いレイアウト (下の画像)
  これらのファイルは次の 「内容に関するファイル」
  1. 共通設定に関するファイル (データ定義その他のマクロなど)
  2. 要旨ファイル
  3. 本文ファイル
  4. 図 (scalable な EPS ファイル) のファイル
  5. 図のよみこみかた (紙上でのサイズなど) 指定のファイル
  6. 図の説明に関するファイル
  を読みこみ， それぞれ指定したレイアウトにそって出力できるのである．
• まあ LaTeX だけで完結すれば一番よかったんだけど， 図に関する 5. と 6. についてはちょっとややこしい部分があるので， Perl による助けが必要だった． すなわち， fig.tex なる 4. と 5. の情報を統合したファイルを作っておき， ここで簡単な自作 Perl スクリプトでフィルタリングして 二種類の体裁にそれぞれあわせた「図に関するファイル群」 を生成したのである．
• まあ， このあたりのシステム化作業に多少の時間がかかったのは認めます …… が， こういう文書まわりに関しては LaTeX と何かテキストのあつかいにすぐれたプログラミング言語があれば 「何とでもできる」 感はすごく大事だと思う．
• 「阿呆らしい不本意なことをやらされてる」 感を軽減するためにとりくんだ この数時間の作業， つまり

  内容はいっさい変更せず見てくれだけまったく変える

  を「わーぷろ」のたぐいでやろうとすると …… おそらくそれはかなり絶望的な試みだよね． 内容と見てくれがぜんぜん分離されてないから． それゆえに私は「わーぷろ」は使えないわけで．
• わざわざこういうこと書いてるのは， 「何でもわーど」 ぽりしーに寸毫の疑いすらもたないヒトたちにちょっとでも 「うーむ……」 と思わせることができれば楽しいからである．
• 本日は， 甲山さん・鍋嶋さん・宮田さん・田辺さんが 苫小牧 まで「試掘」に言ってたそーで …… 7F 実験室で掘り取られた実生についていろいろと教えていただいた．

[広葉樹とかの実生]

オオモミジは細いのが広くひろがるけれど， ミズナラは太いのが地中にささっているといった 「地下部のカタチ」 はどう特徴づけたもんかな? 細根は重要だけど完ペキに採取してくるのは難しい． となると， (ある程度以上の太さの) パーツに分離してその重量・長さをはかって 地下部における「資源分割」をみるとか?

• 皆さん土曜日の夜おそくまで働いているけれど， 私はこそこそと撤退． 2050 研究室発． 2105 帰宅． 運動． 体重 73.0kg． 晩飯．
• [今日の運動]
  • エアロバイク 40 分間．
• [今日の食卓]
  • 朝 (0830): 米麦 0.6 合． ゆでサケ． ハクサイ・ニラ・シイタケの炒めもの．
  • 昼 (1230): 米麦 0.5 合． 海藻スープ．
  • 晩 (2300): スパゲッティー． タマネギ・シイタケ・豆腐の炒めもの．

2006 年 06 月 04 日 (日)

• 0700 起床． 朝飯． コーヒー． 洗濯． 怠業． 1105 自宅発． 晴． 1120 研究室着．
• そうか， resume 時の APM 制御 (USB まわりとか) は /etc/sysconfig/apm-scripts/apmcontinue だったな (ぎょーむ日誌 2004-10-23)．
• 地環研 A 棟 8F 非常階段から見える最高峰である手稲山の山頂 (標高 1023m)， すでにほとんど残雪ナシか …… 昨年 6 月初頭 のほうが残ってるようにみえる．

• GLMM 解説記事の原稿， よーやく内容のみなおし． けっこうひどい文章だ， といまさらながら気づいた． 修正．
• 昼飯調達の旅． 研究室にもどって昼飯．
• また原稿みなおしのつづき．
• 説明の順番， またかなりいれかえてしまう．
• まあ， ひととーりはなおせたかな， ということで撤退． 2110 研究室発． 2125 帰宅． 晩飯．
• 月曜日の午前中は歯医者．
• [今日の運動]
  • うんどう休養日．
• [今日の食卓]
  • 朝 (0830): バゲット．
  • 昼 (1420): 研究室お茶部屋． ほっかほか亭幕の内弁当 470 円．
  • 晩 (2200): 米麦 0.7 合． タマネギ・シイタケ・豆腐の炒めもの．

2006 年 06 月 05 日 (月)

• 0740 起床． 朝飯． コーヒー． 午前中は歯医者に行くため半休． 洗濯．
• GLMM 解説原稿みなおし． どうも歯というより歯医者が気になって集中できない．
• 1020 自宅発． 晴． 歩いて一分で おおつか歯科 着． 問診票記入． すぐに診察室へ．
• 土曜日に欠落した左上の親知らず， けっきょく 抜歯 することに． いろいろと詳しく説明してもらったんだが …… どうもはえかたがヘンで， 下の歯茎とかとも干渉してるため． あと私の上下あごは左右にずれてんだけど， いわれてみるとたしかにずれを補正しようとすると， くだんの親知らずがひっかかる．
• 口蓋全周というか半周の X 線写真撮影． 左上親知らずは頭蓋骨についてないので抜くのは問題なかろう， とのこと． こいつはヘンなふうにはえてるので下の歯茎を痛めつけているとか． ついでながら他の親知らずの現状:
  • 右上: 完全に歯茎の中に完全に埋没 (X 線写真にその様子がくっきりとうつってる)
  • 左下: とくに問題ナシ
  • 右下: 大学院生のころ九大ちかくの歯医者でぬいた
• 「また抜歯の日を予約しますか? 今日いきなり抜くってのも『覚悟』が……」
  「いや，今日は抜歯かもしれない，と思っていましたんで」
  「ならハナシが早い．5 分か 10 分で終わりますよ」
  ということで， 抜歯そのものは麻酔注射の時間もふくめて， たしかに 5 分もかからずに終了． 鉗子でごりごりと． 右下ぬいたときに比べると， しごくあっさりとぬけたかんじ．
• 抜けた歯をみせてもらったんだけど， 歯根がやや未発達なまま? というような． そして磨きにくい位置にあったんで， 咬合面は虫歯にやられてしまっている． やれやれ．
• 1105 診療終了． 化膿止めと痛み止めの内服薬をもらう． 診療費は
  • 基本診療料: 180 点
  • 投薬処方等: 180 点
  • 検査処置等: 787 点
  負担率 30% の保険点数 1147 点の定率負担金額 3440 円． 明朝も消毒のため通院．
• 1107 帰宅． まだ局所麻酔がきいてるせいか， ぜんぜん痛くない． たしかに例の左上親知らずがないと， アゴの動きは少し快適になったような気がする． 出血は 30 分もしないうちにだいたい止まった．
• 昼飯． 抜歯あとちょっとだけ痛む． だしてもらった化膿どめ内服薬のむ． そういえば， 半年ほど前にとある院生が抜歯したんだけど， その直後は顔がはれあがっていたな． 1240 自宅発． 晴． 1255 研究室着．
• あいかわらず GLMM 解説記事みなおし． いつまで続くのやら …… って本日中には提出せんといかんわけだが．
• 不思議なことに抜歯部分はほとんど痛くない． 麻酔の影響なのか， アタマがややふらふらする．
• 原稿みなおしもいいかげんイヤになってしまったので， 投げだすことに． 1627 日本生態学会誌編集長の大串さんあてにメイルで PDF ファイル 「[解説記事] 『個体差』の統計モデリング」 (久保・粕谷) を送信． 一件落着 …… なのだろうか?
• 私の横で web server 機の冷却ファンが不調． A801 であまってる電源ユニットをいただいてくる． あとで交換しよう．
• GLMM 解説記事さぽーとぺいじでも作るか …… といっても現時点ではナゾの登場人物 「Q 氏」の観測データファイルぐらいしか置くべきものがないわけだが．
• さてさてアカマツべいぢあん不思議世界の検討に復帰するか ……
• ここ数日， そのへんふらふらと歩いたり飯くってるときとかに， このへんのモデリングについては検討を続けていた． まあ， お望みの結果がでるかどうかはともかく， モデリングそのものは問題なさそう， という気がしている．
• ということで， 計算の実装を検討すべき段階なんだが …… うーむ， JAGS にやらせるには， 私の JAGS 理解がぜんぜんおよんでいないな． こういうので計算できればラクなんだろうけど．
• いや， そう簡単にあきらめずにもう少し調べてみるか …… 難所は stem の Gibbs sampling なんだよな． 皮肉なことに， というべきかこの部分を Gibbs sampler として書くにはどうしたらよいか， といったことはわりと簡単に思いつく． ただし， それはむちゃくちゃに計算コストがかかるサンプラー (おそらく他の方法はない)， というあたりが問題なんだが． JAGS の adaptive rejection sampling 法とはあいいれないシロモノだ．
• 1811 web server 機の電源ユニット交換作業開始． つけかえそのものは 10 分ちょいでできたんだけど， 再起動すると例によって HDD のファイルシステムチェックに時間がかかる． /dev/hda は bad block だらけなんで． そろそろ買い替えんといかんのかな． 1902 web server 復旧．

[web server 機現状]

電源ユニットの空冷ファンがときどき止まるので， 筐体よこの「ふた」をあけて温度が上がらぬよう運転 (左上)． 電力ケイブルはこのようにわしゃわしゃと …… まあ， 見た目ほどにはややこしくないんだけど (右上)． そして時間かかる bad block 修復 (左下)．

• 1930 研究室発． 1950 帰宅． 晩飯．
• [今日の運動]
  • 抜歯うんどう休養日
• [今日の食卓]
  • 朝 (0830): 米麦 0.6 合． ネギ・卵の炒めもの．
  • 昼 (1220): 米麦 0.5 合． 海藻スープ．
  • 晩 (2100): 米麦 1.0 合． 麻婆豆腐．

2006 年 06 月 06 日 (火)

• 0730 起床． 朝飯． コーヒー． 昨日の抜歯あとの消毒のため歯医者にいくので， 今朝の午前中も有給休暇の半日を消費する． 0920 自宅発． 曇． 歩いて一分でおおつか歯科着．
• 0930 から診療． 抜歯あととくに問題なし． すでに麻酔ふらふら感もないし， ぜんぜん痛まないので助かっている． これはくだんの親知らずがちゃんと成長してなかったおかげか …… よくわからんけど， 歯の中の神経がすでに死んでた， とか?
• 抜歯に関しては， かくのごとく damage ほとんどなかったし， アゴの動きが自由になったし (歯一本でこんなに変わるのか， というのはなかなか興味ぶかい)， みがきにくい歯は無くなったし， ということで今のところは良いことばかり． 出費はしてるけど．
• 本日は， 「じゃ，ついでに」 ということで歯石とりとか． 歯科衛生士さんが電動歯ブラシの親玉みたいなやつとスケイラーで …… まあ， ようするに歯の表面を研磨し歯間を掘削してるわけで． 以前に歯石とったときには， 口中けっこう血まみれになったけど， 本日はなぜか出血なし．
• で， 本日は Crest whitening (お， Wikipedia にも) なる過酸化水素水いり歯ミガキぺいすとの試用をすすめられて終わった． 経緯が少しひねくれてて
  • 歯石除去してみると右上の犬歯にまたしても 小さな虫歯が発見された (いやはや)
  • ところで私の歯はコーヒーのせいか， かなり茶変している (これまたいやはや)
  • いまこの虫歯を治療するとつめものも まわりの歯の色にあわせねばならない
  • ところでこの医院では， その Crest なるペイストで歯をみがけば 簡単に白くなりますよとススめている
  • ということで， 少し白くしてから治療しては?
  という流れなわけで． なんといいますか， ちょっとまにあっくなかんぢの歯科医先生に 当方の 「じゃ， 実験してみますか」 なノリを見すかされてるわけで．
• 歯の色まめ知識: ヨーグルト食べたあとにコーヒー飲むと歯が茶いろくなる． 理由は乳酸が歯の表面をとかすので， コーヒーに染色されやすくなるとのこと． 果物に入ってるクエン酸も同様らしい．
• 本日の会計． 基本診療料 38 点， 投薬処方等 59 点， 検査処置等 384 点， の定率負担額は 1440 円． 投薬はなぜか「うがい薬」イソジンで， これは口中の消毒のためとのこと． あ， それから例の Crest 2000 円 …… えーと， 値段が高いのは， 日本では認可されてなくて医薬部外品として輸入してるから， とか． さて， 私の歯には効力を発揮するのか?
• 金曜日にまた抜歯あと消毒のため来院せよとのご指示． 1005 おおつか歯科発． 1020 研究室着．
• 1030 から 研究室セミナー． 本日は林君のアラスカ森林火災あとのミズゴケ上の black spruce 定着とか調べますというハナシと， 矢澤さんのマレイシア pasoh で熱帯林の「傾き」を測定しますというハナシ． アラスカのほうではミズゴケの「もりあがり」を三次元的に測定して， 生産力とか推定してみるとか． いっぽう熱帯ではレーザー測距で樹高測定か …… どっちもたいへんそうだ．
• 1140 ごろ終了． メイルかきとか．
• 昼飯． 昼飯後は， なんとなくお茶部屋ながしの掃除． 田辺さんと樹木実生の地下部測定相談．
• アカマツ検討のつづき． 今朝も歯なぞをごりごりと研磨されつつ考えていたことは， 改良アカマツ成長モデルを Bayesian ふうに考えるとして，
  • MCEM 法: 加重 resampling + (nlme() か lmer() による最尤推定)
  • MCMC 法: ぜんぶ Metropolis-Hastings 法による sampling
  といった計算方法があるんだけど， 私にとっての利害得失がよくわからん， といったことども． ここでてきとーに整理してみると (「最尤推定」と書いてるのは厳密には対数尤度の期待値の最大化， か)，

  	MCEM 法	MCMC 法
  計算速度	最尤推定が遅そう	収束むちゃくちゃ遅そう
  予測される計算難所	最尤推定	ややこしーバグ
  事前の調査・実験	かなり必要かも	特になし
  実装	まだマシか?	めんどくさそう
  計算の説明	めんどくさそう	まだマシか?
  知名度	どマイナー	ややマイナー

  うーむ ……
• どちらの方法をとるにせよ， 尤度関数まわりはよく検討せんといかんので， そのヘンの数式を紙に書いてひねくってみたり …… やっぱあまり簡単にならんよな． もし MCEM 法でいくとしても， できあいの線形混合モデル計算関数におしこむのもたいへんか? まあ， R なら optim() + integrate() ワザとか使えるけどなぁ (計算速度はいよいよ遅くなるわけだが)．
• といった検討だけで本日はオワってしまった． 1845 研究室発． 1900 帰宅． 体重 73.4kg． 運動． 抜歯あと化膿どめの薬のせいか， 体調いまいち． 晩飯．
• [今日の運動]
  • エアロバイク 30 分間．
• [今日の食卓]
  • 朝 (0750): 食パン．
  • 昼 (1310): 研究室お茶部屋． 米麦 0.5 合． 麻婆豆腐．
  • 晩 (2150): 米麦 0.7 合． ネギ・ブナシメジ・コンニャク・鶏ササミの煮物．

2006 年 06 月 07 日 (水)

• 0700 起床． 朝飯． コーヒー． 0835 自宅発． 曇． 0850 研究室着．
• アタマが整理されてきたような気がするので， アカマツ成長の MCMC 計算プログラムかきに着手してみる． 例の屋久島葉死亡モデルの計算中核部の改造． なんと Perl+SWIG でかかれているこのコード， 階層ベイズモデルを計算するために けっこうややこしい階層処理が入っていて (当時はそうするしかないと考えていた)， このへんはたぶんいらないだろうからりすとら． 使える部品 (class) を整理する．
• で， どんどん実装していくわけだが …… コード書きなおしてるうちに， いまさらながら 今回のモデリングに関して， ひじょーにまぬけなしくじりに気づいた． cohort2000 の成長が悪い， ということと cohort1999 の成長がよろしくない， というのは数式の上では識別がつかない! つまり推定不可能だな．
• しょうがないんでモデルを再検討．
• …… にとりかかろうとしたところ， 塩寺さんの熱帯葉死亡モデリングのデータ解析こんさる． うーむ， Cox の比例ハザードモデルを R で推定計算しておられるんだけど， 「与えられた期間」 より長い base line hazzard 関数を計算できないのか …… ということで， coxph() による推定はあきらめていただき， survreg() でワイブルな死亡曲線とかを使ってもらうことに．
• アカマツ検討にもどる． どうも茎直径分布を考えるハナシはダメそうだな …… ということで， パイプモデル系の前提にもとづくモデリングはやめてみる．
• 枝さきにおける重量分配モデル， みたいなのを考え， いわば「欠測」データに該当するのが成長前の cohort1999 の大きさ， と考える． cohort1999 の 2000 における成長， そして cohort1999 と cohort2000 における資源分割を考えて， と．
• で， ここで (欠測をうめるための) 私のおきにいりアイデアである resampling はあきらめて， なにかパラメトリックな事前分布におきかえてしまおう． ここで興味ぶかいのは， この欠測をうめるための事前分布は全個体共通としてよいところだ． なぜかといえば， 「成長の遅さ」の個体差と 「初期状態の貧弱さ」の個体差は識別できないからである． つまり今度はさきほどの識別不可能性を逆に悪用してやろう， と．
• さてさて， 問題がここまでさらに単純化されてしまったら (というか resampling をあきらめたので) JAGS とかでもさらさらと MCMC 計算できてしまうモデルが構築されつつあるような気がする． とりあえず昼飯．
• Perl 版の自作 MCMC 計算プログラムの改善はとりあえず中断． JAGS 初めての実戦投入へのシークエンスを開始する． まず JAGS のコードを読んで BUGS 文法の研究． R とみかけが似ているわけだが …… 文末に ; がいるのかいらんのかよくわからん (無くてもいちおう動作する， とわかった)．
• よしやってみるか， とコード書いてみる． ついでに架空データも． 30 分もかからずにひととーりできた． けど， いきなり syntax error ではじきかえされる …… ~ と <- の使いわけに注意． 下のごとくでとりあえずは動く．

model {
	for (i in 1:n.tree) {
		bi[i] ~ dnorm(0, iv.bi);
		y[i] ~ dnorm(
			exp(b0 + bt * treatment[i] + bi[i] + x[i]),
			iv.g
		);
		z[i] ~ dnorm(
			exp(a0 + at * treatment[i]) * (y[i] - exp(x[i])),
			iv.a
		);
	}
	b0 ~ dnorm(0, 1.0e-6);
	bt ~ dnorm(0, 1.0e-6);
	a0 ~ dnorm(0, 1.0e-6);
	at ~ dnorm(0, 1.0e-6);
	iv.bi ~ dgamma(1.0e-3, 1.0e-3);
	iv.x ~ dgamma(1.0e-3, 1.0e-3);
	iv.g ~ dgamma(1.0e-3, 1.0e-3);
	iv.a ~ dgamma(1.0e-3, 1.0e-3);
}

• 架空データつかった JAGS の MCMC 計算試運転 …… あっさりできた． burn-in 10000 step の sampling 10000 step で計算時間 15 秒 (PentiumM 1.6GHz)． 作図は RjpWiki の説明 とかが役にたつ． 架空データに設定しておいたとーりのパラメーター事後分布がでたな …… なかなかすばらしい．

• さて， あとはこいつに本ものデータをどう食わせるか， どういう推定結果になってしまうのか， というあたりだが …… ちょいひと休み．
• しかし， この JAGS の Gibbs sampler むちゃくちゃ速いなあ …… まあ， この架空データは実データにあわせて個体数 8 という状況なんだけど． 私が rejection sampling の関数とか書いたら， すごく遅くなりそうな気がするんですが．
• さて， 実データを JAGS 用データファイルを作るには， R をもちいる． JAGS はそのように設計されているんで． R でデータを読みこみ， 必要な部分をぬきとり， dump() をつかって JAGS 用データファイルを作る． たとえば， こういうかんぢになる．

  "n.tree" <-
  as.integer(8)
  "treatment" <-
  c(1, 1, 1, 1, 0, 0, 0, 0)
  "mx" <-
  c(2.16615384615385, 2.32212389380531, 2.83478260869565, 5.38701298701299, 
  2.96280991735537, 5.14838709677419, 3.2787610619469, 5.78448275862069
  )
  "y" <-
  c(28.2, 15.2, 31.8, 37.2, 26.1, 22.9, 21.5, 25.1)
  "z" <-
  c(17.6, 16.4, 16.3, 24.4, 23.9, 22.8, 24.7, 30.5)
  

• で， 2000 ごろまでかかってべいづなモデリングの 工夫・試行錯誤・悪戦苦闘・七転八倒 …… JAGS が使えてよかった． 自作の Perl 版 Metropolis-Hastings sampler でやってたら， おそらくさらにきちがいじみた状況を楽しまねばならぬところであった．
• けっきょく「欠測」部分の事前分布をどうしたらよろしかろう， というのが面倒の焦点であった． resampling まではいかなかったけど …… 全個体共通の prior というアイデアはぜんぜんダメであった． 今晩得られた最良のモデルは， (やや不本意なものがあるんだが) 事前分布を個体ごとに準備し， その平均値を cohort2000 で決めてやる …… まあ， いったん放棄した resampling アイデアの一部を流用してるわけで．
• そして cohort1999 に関する処理の影響は「原理的に」 (とでもいうべきか) 依然として推定計算できない． にもかかわらず (なぜ cancel out できるのか説明するのが面倒なので) それを含めたまま計算させている …… まあ， いづれにせよあってもなくてもよいパラメーターなんだけど．
• それに対して， いちばん先端部の cohort2000 への資源分配については (個体差その他を考慮してなお) 明らかに減っているといえる …… なんとか 「使いものになる」結果が得られたところで本日は終了． ちょっと無理矢理な部分もあるんだけど．

model {
	b0 ~ dnorm(0, 1.0e-1);
	bt ~ dnorm(0, 1.0e-1);
	a0 ~ dnorm(0, 1.0e-1);
	at ~ dnorm(0, 1.0e-1);
	iv.x ~ dgamma(1.0e-1, 1.0e-3);
	iv.g ~ dgamma(1.0e-3, 1.0e-3);
	iv.a ~ dgamma(1.0e-3, 1.0e-3);
	for (i in 1:n.tree) {
		x[i] ~ dnorm(
			mx[i] * exp(-(at + bt) * treatment[i]),
			iv.x
		)
		y[i] ~ dnorm(
			exp(b0 + bt * treatment[i]) * x[i],
			iv.g
		);
		z[i] ~ dnorm(
			exp(a0 + at * treatment[i]) * (y[i] - x[i]),
			iv.a
		);
	}
}
------------------------------------------------------------------------
1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:
        Mean      SD Naive SE Time-series SE
b0    1.5707   0.186  0.00587        0.01693
bt    0.3129   0.857  0.02709        0.11263
a0    0.2603   0.124  0.00391        0.00788
at   -0.5195   0.173  0.00546        0.01491
x[1]  4.7551   2.801  0.08856        0.35565
x[2]  3.0440   1.914  0.06052        0.23131
x[3]  5.6167   3.294  0.10416        0.41052
x[4]  7.5072   4.469  0.14133        0.58324
x[5]  4.6935   1.564  0.04946        0.16381
x[6]  4.9561   0.930  0.02940        0.09669
x[7]  4.0829   0.977  0.03090        0.09581
x[8]  5.4080   0.923  0.02918        0.09687
iv.x 68.4803 279.713  8.84529       17.38780
iv.g 65.0620 249.301  7.88358       13.45898
iv.a  0.0832   0.572  0.01810        0.02830
2. Quantiles for each variable:
         2.5%     25%    50%     75%   97.5%
b0    1.20958  1.4580  1.603  1.6914   1.861
bt   -0.78595 -0.2827  0.167  0.6841   2.914
a0    0.00788  0.1910  0.260  0.3308   0.512
at   -0.87161 -0.6331 -0.517 -0.4075  -0.187
x[1]  0.24130  2.7066  4.458  6.1832  11.481
x[2]  0.22400  1.6915  2.667  4.3029   6.661
x[3]  0.31323  3.2326  5.214  7.5022  13.006
x[4]  0.48048  4.3256  6.666 10.5067  15.893
x[5]  2.75829  3.1660  4.726  5.5114   7.423
x[6]  3.54275  4.3730  4.994  5.2866   6.750
x[7]  2.77661  3.3414  3.926  4.6285   6.007
x[8]  3.81784  4.7601  5.506  5.8456   6.987
iv.x  0.06327  0.2717  0.597  3.8881 791.402
iv.g  0.00583  0.0205  0.293 15.5221 670.739
iv.a  0.00985  0.0309  0.047  0.0717   0.147
------------------------------------------------------------------------

• 私のぎょーむ時間を日々監視してるらしい院生たちから 久保さんもたまには日没後も仕事やってるんですね， といった勤務評価うけつつ 2015 研究室発． 2035 帰宅． 晩飯の準備． 運動． 体重 73.2kg． 晩飯．
• [今日の運動]
  • エアロバイク 40 分間．
• [今日の食卓]
  • 朝 (0730): 米麦 0.6 合． ネギ・ブナシメジ・コンニャク・鶏ササミの煮物．
  • 昼 (1340): 研究室お茶部屋． 米麦 0.6 合． ネギ・ブナシメジ・コンニャク・鶏ササミの煮物．
  • 晩 (2240): 米麦 1.0 合． タマネギ・ジャガイモ・豆腐のカレー．

2006 年 06 月 08 日 (木)

• 0755 起床． 朝飯． コーヒー． 0920 自宅発． 曇． 0935 研究室着．
• さて， アカマツ成長 Bayesian 計算の山場 (少なくとも最初の山場) はおそらく昨日に終了したので， このあたりの作図あれこれだの説明準備が必要だな． 説明とか書いてるうちにモデルのよろしくない点に気づいたりするだろうし．
• ところで昨日， 各個体の初期状態を生成する事前分布を全個体共通にするアイデア ってのはダメだったと書いたけどそれについて少しばかり補足． 全個体共通の事前分布をあたえても， Gibbs sampler はたいへん賢いことに (というか当然ながらというべきか) 各個体の「実情」にみあった事後分布を生成してくれる． この点では全個体共通の事前分布というアイデアは悪くなさそうだ． しかしながら 問題は私がつくったモデル全体がまぬけであることで， 「初期サイズをあるていど自由に選ばせる」 とした場合に
  • 初期サイズ 1 で成長によって 4 倍になる
  • 初期サイズ 2 で成長によって 2 倍になる
  の区別がつかない (どっちでも同じ)， となってしまうことだ．
• それでも Gibbs sampler は疾走してくれるわけだが …… 値がいつまでたっても「うねり」続けるんですよね， こういういいかげんな統計モデリングでは．
• そこで「初期サイズはこれぐらいの範囲」 というのを事前分布にくみこんでやる必要あるんだけど， これに関しては全個体共通の「これぐらいの範囲」を憶測するのは かえって面倒で， 個体ごとに「ここを中心とする事前分布」 を設定してやるほうが説明がやりやすいのでそちらを使った， ということである．
• まずは Tgif でモデルの概念図を．

• Google image 検索とかやると， ぎょーむ日誌とかから昔の図とかがみつかる．

• cohort の重量分布図の見てくれとかを改善．

• 昼飯．
• 1330-1730 横浜の FRSGC (ぎょーむ日誌用語「独房群」) の 生態系変動予測研究プログラム の 陸域生態系モデル研究グループ の伊藤さん・佐藤さん・加藤さん・稲富さんがこちらに来てくださり， グループ会議 (兼任「小ボス」甲山さんが横浜に行かないので皆さんがこちらにやってきた)． 私もなんとなく参加． たまにこういうハナシをきくのはけっこう楽しい (とくに自分がそういう仕事をやらなくてもよい場合は)． というか， 皆さん楽しそうに研究しておられるように見える．
• 1810 研究室発． 地球フロンティアからこられた皆さん， 甲山さん・高田さんと ひだまり庭 で晩飯． 2145 帰宅．
• 粕谷さんたちは無事に絶海の孤島 男女群島 (地図, 写真) から無事に戻られたよーで．
• 明日の午前中も歯医者． 毎年， 有給休暇はあまるんだけど， このように歯医者で使うと多少は有効活用してるような気分になる …… というのはまちがっている?
• [今日の運動]
  • うんどう休養日
• [今日の食卓]
  • 朝 (0820): 米麦 0.6 合． タマネギ・ジャガイモ・豆腐のカレー．
  • 昼 (1255): 研究室お茶部屋． 米麦 0.5 合． タマネギ・ジャガイモ・豆腐のカレー．
  • 晩 (1830): ひだまり庭． ここはけっこう焼酎が充実していて， 三岳とか伊佐大山とかどんどん飲まれてる． なぜか横浜でぎょーむ日誌とか読んでおられるかたから 「久保さんってホントに飲まないんですねえ!」 と．

2006 年 06 月 09 日 (金)

• 0640 起床． 朝飯． コーヒー． アカマツ計算の改良．
• 0930 自宅から徒歩一分のおおつか歯科で抜歯あとの消毒． どうも回復がおそいらしく， また来週の木曜日も通院しなければならなくなった． そして抜いた親しらず (左上) に対応する 親しらず (左下) もいつの日かぬかねばならぬ可能性あるかも， とのこと …… 「上」がなくなったので「下」がさらに突出して伸びる可能性があるため， という理由だそうで． 本日は 10 分間で終了． 基本診療料 38 点， 定率負担額 110 円． 帰宅．
• 0955 自宅発． 曇． ちょっと雨． 1010 研究室着．
• アカマツ Bayes 成長モデル， サイズに関する事前分布をガンマ分布に変更． 推定結果は (当然ながらというか) あまりかわらないけれど， まあ首尾一貫してるよーな気がしたので．
• R で作図プログラム作り， いろいろと．

[「欠側値」事後分布]

これは「観測してない重量」 の個体ごとの事後分布． 実験処理された個体たちのほうが， その実験処理のパラメーターの影響をうけて， ばらつきが大きくなる．

• アカマツモデル (簡単版) 説明かき．
• 1300 研究室発． 自宅にもどって昼飯． アカマツモデルの成長方法あれこれ検討しながら再登校してたら， 赤信号の交叉点をわたってしまった …… 1350 無事に研究室もどる．
• 説明かきのつづき …… ひととーり できた ので共著者の小林さんにおうかがいメイルをだしてみる． 時刻は 1542． ひとまずは一件落着， か．
• 脱力． ぼけーっとネット雑用など．
• いかん， アカマツ説明であれこれと書き忘れが．
• 1850 研究室発． 1900 帰宅． 運動． 体重 73.4kg． 晩飯．
• [今日の運動]
  • エアロバイク 40 分間．
• [今日の食卓]
  • 朝 (0730): 食パン．
  • 昼 (1320): 米麦 0.7 合． タマネギ・ジャガイモ・豆腐のカレー．
  • 晩 (2210): 米麦 0.7 合． ネギ・エリンギの卵とじ．

2006 年 06 月 10 日 (土)

• 0800 起床． 朝飯． コーヒー． 洗濯怠業．
• アカマツモデルをちょっとひねってみる …… やはりこれ以上はちょっと改良のしようがなさそう， とわかった． JAGS ペイジを分離， しかし現時点ではまだ特になにも追加してない．
• 昼飯くったらなぜかアタマがふらふらしてきたので， 夕方まで横になってくたばる．
• うーむ， 何かヘンなものでも食っただろうか …… と考えても思いあたらず． 昼飯後にコーヒーを飲んだのがよくなかったのか? なぜか植物の防御物質を処理する能力が落ちてしまったので， コーヒーは一日一杯まで， とか?
• 横になってうだうだしてると， アタマの中では幾何学的な問題がうかんでは消える． ちかごろ， 「投影された樹冠の重心」 なる概念をよく聞くわけだが …… これってホントに平面図形の重心になってるのだろうか． 三角形の重心は三角形の面積を等分割する． しかし頂点数が 4 以上の多角形の重心は等分割しない， のかな …… つまり多角形等分割点は重心の十分条件にすぎない， と (三角形なら必要十分条件)．
• アタマふらふらがかなり回復してきたので， 1710 自宅発． 晴． 1725 研究室着．
• りはびりと称してだらだらと 計算生態学 VikiWiki ひねくる．
• そういや， アカマツ写真いろいろ置いてたな．

[実験アカマツ個体]

写真は ここ とか ここ に．

• お茶部屋でなんとなくうだうだと． 北大生協の本屋の値引きぐあいがよろしくないのは， 輸送コストがかかってるからだと正当化されてしまっている と院生が述べているので， いやーそれをいったら 福岡あたりだって同じぐらい輸送コストはかかってるハズですよ， とのべてみたり．
• お. (モデル選択規準 AIC の) 赤池弘次さんが今年の京都賞受賞者に決定 ， か (google news とか)． すごく妥当だと思う． 赤池さんって 1927 年うまれだったんだ． そういや最後あたりの江田島の海軍兵学校生徒だったということだから， たしかにこのへんか．
• 仕事は何もススまぬまま撤退． 2120 研究室発． 2140 帰宅． 晩飯．
• [今日の運動]
  • うんどう休養日．
• [今日の食卓]
  • 朝 (0830): 食パン．
  • 昼 (1230): 蕎麦．
  • 晩 (2210): タラコスパゲッティー． サンマかばやき． リンゴ．