ぎょーむ日誌 2006-11-(11-20)

2006 年 11 月 11 日 (土)

• 0700 起床． 曇天． 寒い． 10 時間以上ねたので回復 …… 昨日， こんな会話が何度か独立に反復されてたんだよね．
  「久保さん，徹夜ですか?」
  「そういうよれよれ状態に見えるかもしれませんが， そうではありません」
  朝飯． コーヒー． 洗濯． 怠業． ぎょーむ日誌生成 Perl スクリプトの修正 (とはいえ browser 見た目は何もかわらんが)， とか．
• 1050 自宅発北大構内走． 1150 帰宅． 体重 72.2kg． 昼飯．
• 1310 自宅発． 雨． 1325 研究室着．
• アカマツ野外データのほうは， よくみると 7 箇所の異なる場所から得られた 21 個体だと気づいたので， このサンプリング場所の random effects を考慮している 非線形混合モデル， とゆーあまりみこみのないデータ解析をつづけてしまった． というか， library(nlme) の使いかたの練習， というか． これ， すごくややこしくなってますよ (と思って nlme 本 ながめてみたら， じつは今に始まったことではなく昔からややこしいシロモノだった， とわかった)．

[光合成速度 vs 葉重量]

nlme() に微分不可能な関数を計算させるのは面倒 (不可能?) なので， みょーな漸近線をもつ関数をでっちあげてみた． そいつは selfStart() 関数として定義されねばならぬ． さらにデータも面倒で， ふつーの data.frame() はダメで groupedData() なるものを使わねばならぬ (非線形関数などの場合)．

• で， これだけ努力したのに， サンプル場所によって切片の異なる 2 パラメーター線形混合モデルに AIC 的に敗北 (「惜敗」?) しました， と．
• ならばもうちょい単純化してみたらどうだ， と．

[光合成速度 vs 葉重量]

こんどは formula じたいがまがりなりにも線形なので， selfStart 関数つかわず推定計算できた． この状況だと， なぜかしら (互換性維持?)

  random = a + b + c ~ 1 | plot

といった古い nlme() 書式も使えた． で， この二分割モデルは AIC 的にはさらに劣悪， と． こちらの 2 パラメーター線形混合モデルでは サンプル場所によって切片・傾きが異なるとしているので， 上とはびみょーに異なった推定結果となっている．

• とはいえ 「よくわからんどうし」 である非線型 vs 線形混合モデルの比較も無意味というか …… もっと単純なモデルとしては 「葉重量つねに一定， ただし観測場所ごとの random effects あり」 というもので， さすがにこれは AIC 的にあまりよろしくない， と．
• といったことなど簡単報告メイルしてから， 1920 研究室発． 雨． 1940 帰宅． 晩飯．
• [今日の運動]
  • 北大構内走 60 分間．
• [今日の食卓]
  • 朝 (0740): 米麦 0.4 合． ダイコン・豆腐・ワカメ・煮干の味噌汁．
  • 昼 (1245): 明太子スパゲッティー．
  • 晩 (2100): マカロニ． ハクサイ・ネギ・エリンギのスープ．

2006 年 11 月 12 日 (日)

• 0740 起床． 雨． まだ雪になってない． 朝飯． コーヒー．
• 怠業ぎみに R で遊ぶ． ふーむ ……

[しつこく光合成速度 vs 葉重量]

昨日でっちあげた， 奇妙な漸近線をもつ 3 パラメーター非線形モデルを 実験アカマツ集団データにあてはめてみる． これは random effects を仮定できる構造がないので (この正規分布モデルでは「個体差」は測定誤差と識別不能)， nls() を使ったパラメーター推定となる． (昨日 nlme まわりでなやまされた) くだんの selfStart 関数などは使いまわせる．

• そしてどろなわ的な 生態学会自由集会 準備． 次は MCMC & ベイズモデルの一番初歩的なあたり， をやろうかと考えているんだけど， 話題提供してくださるヒトがなかなか見つからなくて， ですね．
• 雨音にまじって窓をぱしぱしと叩く音が． 札幌市内， まあ初雪といってもいいんぢゃないでしょうか．
• 昼飯． 1310 自宅発． 雨， かと思ってたら途中から霰になった． 風つよし． 1330 研究室着．
• 自分では使わないのに HTML editor について調べねばならぬことに …… SeaMonkey は無難にインストールできるけど， ホントに必要最小限の機能だけ， というかんぢ．
• なんかもうちょっと良さげなものないかね， と調べてみると …… CSS まわりとの関連で言うと Amaya とかが W3C 謹製おカタいかんぢで良さそうな? 日本語つかえなくてかまわんわけだし． Linux 用 RPM 版あるんだけど， Vine 3.2 ではへっぽこなことにインストールできんかった (libstdc++ が古い!)． 自分で使うわけでもないので build する根性もなし．
• Vine で apt-get install できる Bluefish はなかなか興味ぶかいエディターだな (server upload 機能とかないみたいだけど)． なぜか R にも対応してるし． あ， Windows 版だけないのか． こりゃ， 楽しい． しかし却下．
• まあ， いろいろあるわけだし， ご本人に選んでもらいましょうかね． 私自身は vim 上で HTML タグを書いてたりします． こんなのススめたら人権侵害とか言われるかもね．

• 土日の A 棟はスチームヒーターが切られていて寒い …… これを本日の「仕事が進捗せぬ理由」としよう． 窓の外は， まあ雪といってよい状態でしょう． 知床調査隊とかどうなってることやら．
• とかうだうだしてる場合でもないんで， 母子里林冠作文の修正にとりくむ．
• 次の生態学会松山大会の データ解析自由集会， どうやら粕谷さんと私だけの話題提供になりそうだ． まあ， 二時間 なんで， これぐらいの話題提供者数が適切な気もする． 内容は階層ベイズモデル， MCMC 計算， 経験ベイズ法， といったべいづまわりである．
• 自由集会の採否は抽選できまるんで， 「概要」作文をがんばってもしょうがないんだが …… とりあえず 500 字ほど概要説明の作文．
• 1940 研究室発． 2000 帰宅． 晩飯．
• さて， 生態学会大会といえば一般講演のほうだが …… ふと思いついて昨年度のぎょーむ日誌などを見なおしているうちに， アリ研究の発表は無理 と突然にサトる． というのも， 先月・今月にある程度すすめる予定だったわけだが， 下うけ仕事 3-4 個であっさりとその予定が流されてしまったわけで …… そしてこれから三月末にむけてこういう生活がつづきそうだ． 今年度は M2 が 4 名， M1 が 6 名． 昨年度のぎょーむ日誌みて気づいたのは M1 下うけにもかなり 体力・知力・気力・時間を費していたらしい， ということで …… 今年度は昨年度よりも， 少なくとも突発的に難題が発生する rate はそうとうに高まっていると 考えるべきだろう．
• 言わずもがなのコトなんだけど， D なヒトたちも多数いる． ただし D なヒトたちに関しては 「こういうふうにやればいいからあとは自分でやってね」 戦術が使える (先方の交渉能力が高いと通用しないが) …… M なヒトたちあいてにこの術策をもちいると 「めちゃくちゃ」なことになりがちで， とりかえしがつかなくなることがある． もちろん D なヒトたちも「めちゃくちゃ」な状況におちいりがちだが， 自らのぞんでここまで進学してしまったヒトたちには 「自業自得で自滅していく」特権が付与されているように思えるんで．
• こういう状況で「アリコロニー間戦争」モデリング研究みたいな， 観測データが何もない， 先行研究が何もない， 現時点ではほとんど誰も興味をもっていない， 結果が見えないだけでなく「何を示すべきか」もよくわからない， 丸なげ下うけぎょーむが完結する， とは考えにくい．
• えーと， つまり
  • 今年度中にアリ研究に着手し， 計算いんふらの整備は始めなければならない (というか，システムはほぼ完成状態までもっていく?)
  • それっぽい計算結果も出さねばなるまい
  • しかし松山大会にまにあうかどうかは …… まあ， 無理でしょう
  ということだ．
• さてさて， それでは一般講演に何をもってくるか， という問題だけど …… 牛原さんとススめていた 屋久島葉寿命の階層ベイズモデルとかどうだろう? 昨年 11-12 月 から 今年 1 月 にかけてすごく苦闘して計算したのに， そのまま放置されてるあれのことですね．

• まあ， 一般講演もうしこみは 11/17 (金) 17:00 までなんで， 明日あたりに甲山さんと， すでに卒業して就職された牛原さんにお尋ねしてみますか．
• …… てなことで一件落着したつもりになって， 本日はもう寝よう， と思ったんだが …… 「あのややこしい屋久島ベイズモデルを ゐんばぐす に解かせるには， どういう工夫が必要か?」 という問題がアタマの中で全力疾走してしまったので， 真夜中すぎてもぜんぜん眠れなかった．
• [今日の運動]
  • うんどう怠業日
• [今日の食卓]
  • 朝 (0800): バゲット．
  • 昼 (1230): スパゲッティー． ハクサイ・ネギ・エリンギのスープ．
  • 晩 (2130): 米麦 0.7 合． ダイコン・豆腐・ワカメ・煮干の味噌汁． ニラ卵．

2006 年 11 月 13 日 (月)

• 0750 起床． 昨晩は睡眠時間を屋久島べいづモデルの検討にとられてしまって， ねむい． 朝飯． コーヒー． 0910 自宅発． うっすらと積雪． 晴． 0925 研究室着．
• いろいろとまたメイルかきとか．
• 輪読会の予習． 昼飯．
• 1300 より 輪読会， 本日は IPPB 第五章， まああえていえば植物の定着過程の個体群動態みたいなところ， か． 担当は小山さん．
• 昨晩， 屋久島葉寿命モデルがどーのこーの， と書いた事象とは独立に， 牛原さんが北大に来られたので驚く (何やら甲山さんに用事があったとかで) …… ついでに松山大会まわりのハナシも．
• 1845 研究室発． 1900 帰宅． 体重 72.4kg． 晩飯．
• そしてまた真夜中ちかくなると， なぜかアタマの中に計算問題が「わいて」きて， ですね．
• [今日の運動]
  • そしてまたうんどう怠業日
• [今日の食卓]
  • 朝 (0820): 米麦 0.5 合． ダイコン・豆腐・ワカメ・煮干の味噌汁．
  • 昼 (1250): 研究室お茶部屋． 食パン．
  • 晩 (2130): 米麦 0.7 合． ダイコン・豆腐・ワカメ・煮干の味噌汁． コマツナあえもの．

2006 年 11 月 14 日 (火)

• 0815 起床． 朝飯． コーヒー． 0915 自宅発． 曇． 0930 研究室着． 輪読会予習．
• 0900 より IPBB 輪読会， いろいろな予定変更のあおりをうけて二日連続． 本日の担当は小林君． 今回は行列モデルみたいなところなんだが …… この IPBB は「植物の個体群動態の教科書」と銘うっておきながら， 個体群動態に関する Silvertown の説明はぼろぼろで (とくに数式まわり)， ここだけ読んでもほとんど意味が取れない． 一般化もヘンであるし， ちっとも具体的でもない．
• そして植物の個体群動態モデルが行列モデルだけってのはどうよ …… この章で「つぎはぎ」的に環境変動性だの密度依存性だのもちこむと， 行列モデルのご利益 (「固有値とそれに対応する左右固有ヴェクトルで何もかもわかりまーす」 とかいうアレ) がまったくなくなるわけで． いまどき生態学会大会の一般講演でも 「固有値計算したら，この集団は無限大にむかって 指数関数的に増大していくとわかりました」 とかないよ．
• さらに植物生態学では行列モデルなど個体群動態モデルの推定に関して， 固有の問題があって …… (たとえばプランクトンの個体群動態研究とは異なり) 植物の場合は
  • 個体が動かない → 空間相関・「ブロック差」の発生
  • 動かないので個体識別が可能 (な場合が多い) → 「同じ個体」と識別されるものから何度もサンプリング (経年変化など) → 「個体差」ある状況での擬似反復 (pseudo replication)
  動かない・個体識別可能という特性が意外な面倒を惹起している (パラメーター推定にあれこれと注意が必要)， ということを行列まにあなヒトたちは気づいてないんだよね． 以前に Caswell さんがこちらに来たときに質問してみたんだけど， 「まるこふだからいいんだ」 とのことで (これはまちがい!)， このへんにある問題を理解していなかった．
• 北大構内昼飯調達の旅． 本日は比較的あたたかい． 研究室にもどって昼飯．
• 夜までなぜか井田君の光合成曲線推定問題につきあう． これはなかなかたいへんで， しかも本日中には部分的にしか解決できなかった．
• 要するに LI-6400 で測定したデータから， 2 処理に合計 7 個体， 各個体 3 シーズンの測定値あり， 光合成曲線を推定しなさいという問題なんだけど， 以下のような面倒がある．
  1. 光合成曲線の関数型の問題 (パラメーター数おおすぎ)
  2. 一個体から何度も何度も (一度の測定で数十個の測定値を一個体からとり， それを同一個体で 3 シーズンくりかえす) → またまた擬似反復 (pseudo replication)
  3. このめんどうな状況下で 「処理が光合成曲線にカタチを変える」 を統計モデルの中に組みこまねばならない
• まず関数型の問題． これは私なども母子里林冠モデリングでも (惰性で) 使っている 非直角双曲線 (Thornley 1976) とやらを使おうとしたんだが …… この 3 パラメーター (呼吸もいれると 4 パラメーター) モデルを表現している数式， 関数としての性質は劣悪で， 私はくわしく知らないんだけど， 光合成測定ぎょーかいではずいぶんとアヤしげな fitting がよく使われているらしい．
• R の nls() はなかなか強力な非線形最小二乗法関数なんだけど …… このタチの悪い Thornley モデルにあてはめるには start = c(...) をうまく選んでやらねばならない (こんなに劣悪な関数型でなければ， 初期値はもっといいかげんに設定しても問題ない)． しかし 21 セットもいちいち初期値を吟味してられない． とくにこのあと「個体差」考慮した mixed effects モデルがひかえてる場合には．
• もちろん R には optim() とかあるので Newton 法では苦闘させられる問題であっても， Nelder-Mead 法など 「ゆっくり探しまわる」 方法がいろいろあるわけだが …… これも mixed effects がらみつまり nlme() とのかねあいで却下．
• ということで， まったく十分に現象を説明できている 2 パラメーター (「初期勾配」「最大値」) モデルを使うことにした． 統計学的なモデル選択の観点では， こちらのほうが良い． さらに Pmax * (1 - exp(-f * I)) といった関数型なので， 関数としての性質はだんぜん良い． nls() ならものの数秒で 21 本の光合成曲線をぴたっと推定してくれる (下の図はそのうち 12 本)．

• さて， 次にこれを (たまたまここしばらくあれこれと勉強することになった) nlme() で非線形モデルの mixed model 推定やるわけで …… とりあえず処理・シーズンを考慮しない， 「個体差」だけ考えた推定は何の問題もなく始末できた．
• しかし処理だのシーズン変化だのを 光合成曲線モデルに組みこむところでいきづまった …… くだんのややこしい nlme 本 とかながめてみたんだけど， 複数の fixed effect 要因 (光量 + 処理 + しーずん) あるモデルの計算ができないんでは? という疑惑が ……
• 本日はここまで． もうちょっと調べてみてわからなかったら …… この問題も階層ベイズモデル送りだ． これなら必ず計算できる．
• 1900 研究室発． とある事情 (?) でお供をもうしつけれられて， 地環研最近傍「サンクス」北側の洋食屋「ねこや」に初めて入ることに． 当店は (べつに悪く言ってるわけではなく) つねに「ゆっくり」がぽりしーのようで． しかも本日はピザの日で店内は満席． カウンター席目の前の厨房で何枚ものピザが次々と作られ， オーブンで焼かれて …… 別の席に運びさられていくのをのんびりと見おくる． 一時間ほどして巨大なピザがでてきて， これまたゆっくりと食べることに． 2130 帰宅．
• [今日の運動]
  • うーむ ……
• [今日の食卓]
  • 朝 (0840): 米麦 0.5 合． ワカメスープ．
  • 昼 (1330): 研究室お茶部屋． 北大生協の巻きずし．
  • 晩 (2010): 北 12 西 4 (5?) 「ねこや」 ピザの日， 24cm (?) ピザ． トッピングは生ハム， 生バジルなど． たいへんおいしかったです． ポタージュ， コーヒーつき 850 円．

2006 年 11 月 15 日 (水)

• 0800 起床． 0830 自宅発． 晴． 0845 研究室着． 朝飯． コーヒー．
• なンとなく断わりにくかったので， また査読ひきうけてしまった． そして生態学会 松山大会 の参加 & 一般講演もうしこみ．

==================================================
大会登録番号：0377
お名前：久保 拓弥 (くぼ たくや)
生態学会員（会員番号  3134）
...
公募シンポ企画・講演予定：なし
一般講演：ポスター発表
一般講演の演題：屋久島照葉樹の葉寿命推定の階層ベイズモデル
一般講演の著者：牛原阿海, *久保拓弥, 甲山隆司 (北大・地球環境)
希望セッション（第１/第２/第３）：植物生活史 / 植物個体群 / なし
==================================================

------------------------ 会費情報 ----------------------
大会参加費　2006年11月17日（金）までに払い込み　6500円
2006年11月18日（土）以降　　　　　　7500円
※学会員でない場合は納入時期を問わず7500円

懇親会費　　2006年11月17日（金）までに払い込み　6000円
2006年11月18日（土）以降　　　　7000円
※学会員でない場合は納入時期を問わず7000円
--------------------------------------------------------

• えーと， 金曜日までに 12500 円を郵便局で支払わねばならぬ， と．
• 大会登録番号: 0377 ってことはまだ 377 人しかもうしこんでいないのか! 明後日 17 日 (金) 17:00 が〆切なんですが． この申しこみ速度って (私にとってはイヤな雑事にふりまわされた) あの 2004 年釧路大会のとき より少し低めかも． 2004 年のときは締め切りすぎになぜか大量に追加されたわけだが， 今回はそのテは使えないはず， だと思うんだけど．
• 〆切といえば自由集会はまだ申しこんでいなかったな．
• …… とお茶部屋でうろうろしてると， 問答無用で江川さん机に連行されて発芽実験 (また発芽実験!) データの結果作図 R プログラミング． お， また overdispersion らしきものが …… しかも植物種によってそのでかたが異なるような．
• 昼飯．
• 日没ごろまで， 午前中の発芽データを使って江川さんあいてのお茶部屋 R 個人講習． 一日目はだいたいいつも同じ順序で
  • data.frame というデータ構造がある
  • そこから行・列をえらぶ方法がいろいろある
  • 列は vector
  • R のデータ型にはいろいろある: numeric, integer, character, factor, logical, ...
  • data.frame からさらに条件つきデータぬきとり
  • ラベルつき vector
  • 作図あれこれ
  といったところ． 重要なのは 「本人がとったデータを使う」 というあたり． すでに講習開始時点で江川さんのアタマの中には そのデータ構造がはいっているので， 「それを R であつかうには」 というところにハナシを集中できる． ということで， 私は「R 教えるなら講義形式より個人講習したほうが速い」 と考えているのである．
• 次は平尾君 (苫小牧の平尾君ではなくこちらの平尾君) あいてのデータ解析こんさる． 例によって集団遺伝学的なモデルで， その最尤推定を R でやるにはどうしたらよいか， というような． けっこうややこしそうなので， まずは対数尤度関数をどう書くか， というところにハナシを限定する． ざっと設計してあとは宿題 (これは私ではなく平尾君の! やはり PD はラクだ……)， ということで．
• 1945 研究室発． 2000 帰宅． 晩飯．
• [今日の運動]
  • …… 明日からは昼運動にきりかえよう ……
• [今日の食卓]
  • 朝 (0900): 研究室お茶部屋． 食パン．
  • 昼 (1330): 研究室お茶部屋． 食パン．
  • 晩 (2100): うどん．

2006 年 11 月 16 日 (木)

• 0650 起床． 朝飯． コーヒー． 0815 自宅発． 晴． 0830 研究室着．
• どうも R を 2.4.0 を変えて以来， update.packages() に失敗する場合が多い． この version up 時に， なんか CRAN package 依存性に齟齬が生じているような気がするな． install.packages() で手動で補正すればいいだけなんだけど．
• Trendy-ml 経由で今月末のセミナー案内． なんと今回は埼玉の崎尾均さんがみえるのか． ということで， セミナー案内 を更新する．
• 一昨日の井田君の光合成速度「個体差」問題， 魔宮的な nlme 本 を探索するのはあきらめて， 階層ベイズモデル化 → R2WinBUGS で MCMC 計算， という戦略にきりかえる． 「めんどうになったらベイズ」 というのは今後の下うけ生活でありがちな状況になりそうだ．
• 個体差があって， 個体の処理 (状態) が途中でかわって， 同じ個体から異なる時期に三回サンプリングして …… という状況で測定される光合成曲線はどうなるか， というあたりをまずは BUGS 言語 (R みたいなもんです) で書いていく． この問題， けっこうややこしいよ ……
• R2WinBUGS 経由で呼びだす (Wine 実行環境上の) WinBUGS が出力する意味不明ぎみの error message どもを解読しつつばぐとり．
• ようやく走ったと思ったら， 計算にえらく時間かかる． 1000 MCMC step で 80 秒ぐらい? 地衣類繁殖の階層ベイズの数倍， といったところか …… で， よーやく結果でたと思ったら MCMC 計算がちっとも混交してないし． WinBUGS を使ってなおこの状況とは．
• 北大構内走するつもりがシューズを自宅におきわすれてきた． まぬけだ． しかし， いずれにせよ雨． 傘をさして， 郵便局へ． 生態学会松山大会参加費 + 懇親会費 12500 円ふりこむ． 手数料 150 円． 研究室にもどって昼飯．
• 小山さんの R + Tinn-R が呪われぎみのようだったので， どちらも最新版に更新しておはらい．
• 昨日のつづきで平尾君に R の optim() を使った最尤推定法のやりかたの説明． 要点は
  • 対数尤度を計算する関数を作る (これは平尾君が準備・点検ずみだった)
  • help(optim) をよく読む: とくに， control と ... つまり上記の対数尤度関数にデータなどわたすあたり
  • 「最小化ではなく最大化だ」 と optim() に指示する control = list(fnscale = -1) 指定を忘れない (忘れると最尤推定ならぬ最悪尤度推定になるので)
  といったところか．
• そして川合さん研究室雑用 Wiki (PukiWiki) 使いかたの説明．
• 宮田さんに R の plot() の説明．
• …… という旅路のはてに自分の机によーやく戻れた． 井田君の光合成曲線の階層ベイズモデル， これはあきらかに収束してないよね． してるパラメーターもあるけど． 何か工夫が必要だ． いや， それよりも前にまちがい点検を．

• 窓の外は雪．
• で， 夜中まで (計算まち時間に) A 棟 8F あちこちをうろうろしつつ， 「チェインごとにみょーなところに収束してしまう」 問題にとりくむ ……
• いろいろ試行錯誤してみてついにわかった答えはぢつに簡単なもので， (WinBUGS で実装されてる強力かつアヤしげな) Gibbs sampler とて尤度の「劣悪すぎる地形」 にいったんトラップされたら自力脱出は不可能なので， パラメーターの初期値は適切に選んでやりましょう， というだけのことだった．
• もちろん R2WinBUGS を使ってるときでも R の中でパラメーター初期値を生成して， それを WinBUGS にわたすことができる． ただし， 今回はいくつかのパラメーターの初期値が指定できない場合があって， ですね …… s[1], s[2], s[3] といったパラメーターがあったときに BUGS コードの中で

  s[1] <- 0.0;
  s[2] ~ dnorm(0.0, iv.noninformative);
  s[3] ~ dnorm(0.0, iv.noninformative);
  

  てなことをやってると (s[1] だけが定数であとは 無情報事前分布をもつパラメーター)， どの s[] にも初期値が設定できないんだよね (設定するとエラーがでる)．
• で， WinBUGS では「初期値はおまかせ」モードがあるんだけど， おまかせしてしまうとトンでもない初期値をほうりこまれてしまい， さらに Gibbs sampler ではそのトンでもない値周辺から脱出できなくなる， と． いやはや， 考えてみるまでもなくあたりまえのことだ．
• ということで， 解決策としては統計モデルをとりあえず少しばかり変更して， アブないパラメーターにはことごとく無難な初期値を設定してやる．
• ちなみにこういった MCMC 計算， Gibbs sampler ではなく Metropolis-Hastings 法を使っていれば， かくのごとき窮地には陥らないですんだのであろうか? これも M-H 法の「さんぷりんぐ方法」に依存しているはずで， random walk sampler を使っていればやはり 「自力脱出不可」 におちいっていたであろう． 独立サンプラーならだいじょうぶかな?
• あるいは母子里林冠モデルで使ったパラレルテンパリング法 を使ってみるとか． WinBUGS も複数チェイン計算はやってるわけだけど， どれも同じ「温度」なので窮地からの自力脱出には 何の役にも立たないんだよね …… まあ， 複数チェイン計算のおかげで 「WinBUGS おまかせ初期値が面倒の原因」 と特定できたわけだけど．
• とりあえず長めのサンプリングを命じておいて， 撤退． 2240 研究室発． 雪はやんでいる． 2255 帰宅．
• [今日の運動]
  • なんというか ……
• [今日の食卓]
  • 朝 (0720): 食パン．
  • 昼 (1330): 研究室お茶部屋． 食パン．
  • 晩 (2030): 研究室お茶部屋． 食パン …… けっこう最悪な一日だ．

2006 年 11 月 17 日 (金)

• 0710 起床． 朝飯． コーヒー． 0830 自宅発． 雪がやんで晴． 0845 研究室着．
• さーて， 昨晩かえりがけにしかけておいた計算だが …… 15000 MCMC step でおよそ 1500 秒で終了していたか． そして混交のぐあいだけど …… まだまだだね． これはサンプル長の問題ではなく， 超事前分布なんかの「無情報ぶり」に問題あるというかんぢだ． 昨日あれこれ試行錯誤してるうちに， このへんの variance をかなり小さめにしてしまったからねえ．
• 試行錯誤． なンかみょーなエラーがでると思ったら …… WinBUGS はまぬけなので 1.0e-2 と書いても理解できなくて 1.0E-2 と書いてやらんとダメだったんだよな， とよーやく思いだせた．
• 計算にかまけてると何もかも放擲状態になりがちである． かかるときはぎょーむ日誌をみると， 意外と正気にもどったりする． とりあえず， 自由集会ぎょーむ， 粕谷さんに「これでいいでしょうか」確認メイルを再度だす． なンだか 破滅的にご多忙 のようで． 午後になっても返信なければ当方の独断で申しこんでしまうか．
• 光合成の階層ベイズモデル， 「個体差」の超事前分布や無情報事前分布のばらつきの設定がむずかしい． できるだけ「無情報」にしようとするんだけど， やりすぎると WinBUGS の Gibbs sampler が計算をナゲてしまう， と．
• この計算問題をしちめんどうにしている原因のひとつ， 「季節変化」のパラメーターのくみあわせかたを変えてみる． 幸か不幸か， いろいろと創意工夫の余地があるよね ……
• 現在の計算まち時間は一回 700 秒ぐらい． ところで， R2WinBUGS による Wine 上の WinBUGS の挙動がいまいちヘンだ． しかも私のてもとにある ThinkPad では仕様どーりの動きなんだけど， (別の部屋にあるので) ssh 経由で動かしてる Dell デスクトップ機だとヘンなんだよね． 何が「ヘン」かというと R2WinBUGS で WinBUGS は起動するんだけど， R2WinBUGS が生成したファイルを自動的によみこんで自動的に実行開始， とならない． ということで script.txt を open して Model menu → Script を選んで 計算開始を手動で命じてやる必要がある． ナゾだ．
• じつは昨日はもっとコワれてしまったので (WinBUGS が起動しなくなった) …… 「こりゃ Wine が壊れたのか?」 と思ってまた 30 分かけてコンパイルして再インストールしてみた． …… けど治らなかったんだよね． またいろいろ調査してみてわかったんだけど， ここでじつは腐っていたのは wine 本体ではなく， $HOME/.wine/ 以下だったと判明した (れぢすとりがヘンになっていく， というゐんどーづの仕様もシミュレイトしているのか?)． これを全部削除してから winecfg などで再構築してやると， ともかく WinBUGS は起動するようにはなった． 手動実行ばぐはあいかわらずだけど． ナゾだ．
• Wine + WinBUGS といえば， 森林総研の伊東さんが MacOS X 上の挙動を調べてくださった (Taglibro 2006-11-13) …… 嗚呼， 現時点では Darwine はまだまだでしたか． 参考になります．
• とりぷるチェイン 8000 MCMC step でまち時間 700 秒のはず， だったんだけど 1200 秒ちかく費してよーやく終了． これは上記変更にともなって WinBUGS 内で構成される グラフィカルモデル の構造が複雑になり， そのため尤度評価に時間かかるようになった， と． かんぢんの混交のぐあいはぢりぢりと改善されつつあるように見える．
• ちょーご多忙の粕谷さんからお返事いただいたので， 1225 自由集会もうしこみ …… うーむ， なんとすでに 31 番目なのか? これはキビしい抽選になりそうだ (ハズれると今回は開催できない)． はてさてどうなることやら．

==========================================
【登録番号】 W-031
【タイトル】 データ解析で出会う統計的問題 - ベイズ統計学の考えかた・使いかた
【略称】 ベイズ統計データ解析
【共同企画者】 粕谷英一 (九州大・理), 久保拓弥 (北海道大・地球環境)

【概要】 これまでこのデータ解析自由集会では一般化線形モデル (GLM)・モデル選
択・一般化線形混合モデル (GLMM) をあつかってきた．これらとの関連を検討しなが
ら，今回はデータ解析の強力な道具として近年普及しつつある「ベイズ統計学」を生
態学研究に役だてていく方法について参加者と議論したい．ベイズ統計学は事前分布
と呼ばれる確率モデルを導入することで観測対象を現実的かつ柔軟に表現する統計モ
デリング手法である．計算機のハードウェア・ソフトウェアの発達によって，事前分
布を「客観的」に規定できる階層ベイズモデルや経験ベイズ法が誰にでも手軽にあつ
かえるようになってきている．

この自由集会では「ベイズ統計，最初の一歩」となるような，階層ベイズモデルと経
験ベイズ法に関する入門的な話題を提供し議論の材料としたい．久保はごく簡単なデー
タ解析例題 (前回あつかった「架空植物の種子数の統計モデル」) を階層ベイズモデ
ル化していく過程を示し，その考えかたとマルコフ連鎖モンテカルロ (MCMC) 法によ
る推定計算を紹介する．粕谷は経験ベイズ法を使ったデータ解析の例を紹介する．
==========================================

• なぜかその粕谷さんから宿題が …… 「同じ座標系 (XY座標) に 2 つの多角形があったとき， 重なっている部分の面積を求めるアルゴリズムを 探しています． Rの既存の関数でもかまいません」 …… まあいろいろな解法ありそうだけど， 何か良さげな R package はないかしらんとネット上で検索してみると， gpclib package が見つかった． install.packages("gpclib") して library(gpclib) して help("gpc.poly-class") 読みつつ， 以下のような例題を作ってみる …… うん， うまく計算できてるみたいだ． 便利な世の中になりましたね． といったことをメイルに書いて宿題終了． 所要時間 30 分．

library(gpclib)
x1 <- c(0, 2, 2, 0)
y1 <- c(0, 0, 2, 2)
x2 <- c(1, 3, 3, 1)
y2 <- c(1, 1, 3, 3)
xy1 <- structure(c(x1, y1), .Dim = c(length(x1), 2))
xy2 <- structure(c(x2, y2), .Dim = c(length(x2), 2))
p1 <- as(xy1, "gpc.poly")
p2 <- as(xy2, "gpc.poly")
plot(append.poly(p1, p2))
cat("# area of p1 =", area.poly(p1), "\n")
p12 <- intersect(p1, p2)
plot(p12, poly.args = list(col = "#ff8000"), add = TRUE)
cat("# area of p12 =", area.poly(p12), "\n")

[重なった領域の面積]

解法の手順は， gpc.poly クラスのオブジェクトつまり ポリゴン p1, p2を生成し， その交わり (積集合) ポリゴン p12 を
  intersect(p1, p2)
で作って面積を計算させる， というもの．

• 昼飯． しかし光合成速度の階層ベイズモデル改良の 「次の一手」 が整理できん．
• バックアップをとってから， 「個体差」のいれかたを整理してみる． 「季節」がみっつあるわけだが
  1. 「季節」ごとの「平均的な」光合成曲線がある
  2. 個体ごとに 1. からずれる (ここまでは今までと同じ)
  3. 個体ごとの 2. のずれかたは「季節」よって異なる (今回の変更)
  さーて， 改善となるか改悪となるか …… ともあれグラフィカルモデルはまた複雑になったような気がする (ので計算時間が増えるだろうという憶測はハズれて， 今までと変わりなくおよそ 1200 秒後に結果が得られた)．
• そして R2WinBUGS でよびだす WinBUGS の挙動がまともになった． すなわち， 自動的に起動 & 計算開始． うーむ， 今までのよびだしかたに問題あったのか? 依然としてナゾ．
• 計算まち時間に発注者の井田君と相談． 午前中の計算で「チェイン間差」の少ない結果がでてくるようになったんで， だいたいどういう推定結果が得られそうかわかってきたため．
• さて， 「個体差 & 季節差」 random effects 的にした計算結果だが …… どう考えるべきかな． deviance は「今までの計算は何だったの?」 と言いたくなるぐらい改善された． つまりあてはまりは格段に断固としてよくなった (DIC 的にも)． 処理の効果なんかも収束は悪くなさそう． しかし最大光合成速度はあいかわらずチェイン間差がでかい． これは何か定式化をしくじっているのか， それとも別の原因があるのか ……
• よくわからんけど， 最大光合成速度と直接に「かけひき」をするパラメーターである 呼吸速度から「季節」変化を外してみる． 実際には「季節」変化してるんだろうけど， このデータからそれを推定するのは困難だろう， という理由で．
• この計算結果だけど， deviance & DIC 的には悪くなる， ただし呼吸速度の混交はよく， しかし最大光合成速度はチェイン差があいかわらず． 呼吸速度の混交がよくなったのは， 最大光合成速度からの影響が小さくなったからだよなあ ……
• 最大光合成速度， うごきが鈍いのはしょうがない． というのも， この植物はわりとすぐに光飽和するので， このパラメーターを推定する材料がたくさん供給されてるからだ． 問題はチェイン間差だ． 原因は何か簡単なことだと思うんだけどよくわからん．
• もはや策もつきたので， とりあえず無情報事前分布をさらに無情報にしてみる． これで収束が速くなるなら誰も苦労はしないんだが ……
• そしてこんなふうに計算にかまけていると他のコトは何もできないわけで． 事務雑用． 12 月出張 の振替休日を 12/6 (水) とする． 実際には休めないけど． いや， 何か有効活用できないかな?
• 無情報つよめた計算， ちょっとだけ混交はマシになったけど， まあ実際のところはあまり変わらずだな．
• BUGS コードよくみなおすと「個体差」+「季節差」 の入りかたが少しヘンだと気づいた． 修正して再計算を命じる． いやはや， これからの植物生理生態データの統計モデリングは たいへんなことになりそうな予感． 今までのが「てきとー」すぎたんでは， とも思うんだけど．
• 請求書を (分担者ぶん科研費を管理してもらっている) 九大の事務室に転送する手紙かき．
• うーむ， つくづくこのパラメーターは 「最後に帳尻があってればよい」 といったかんじででたらめな値をとりやがるな …… 関数型をかえて制約を試みる． また再計算．
• …… どうもヘンだと思ったらまたまた初期値問題だった． モデルをかえたら初期値を再確認． われながら進歩がないなあ．
• 初期値をてきとうに設定してやると， 最大光合成速度まわりのでたらめ計算はずいぶんとマシになった …… 初期値なんぞに依存してる点でかなりおそろしいハナシなのではあるが． とりあえず DIC = 2026 ぐらいと記録しておこう．
• 次は …… また variance のいれかたのみなおしか? しかしこれ以上の無情報化は無理． よくわからんので， まあまったく節操ないわけだが， 呼吸まわりを複雑にしてみる．
• 無節操計算終了． 冬に近づくとびみょーに呼吸量が低下する? といった推定も含むようになると， (ホントにいいことなのかどうかはともかく) やはり deviance や DIC の値は改善され DIC = 1850 ぐらいか． 混交のぐあいは下の図のごとし． ぼこぼこになってる山はチェイン間差．

• さてさて次の策だが …… 最大光合成速度を構成する components の一部を exp(...) からとり出してみるか? これで場合によってはマイナスの光合成速度とかありえてしまうわけだが， adaptive rejection sampling のアルゴリズムで駆動される Gibbs sampler はそっちに値をもっていかぬはず， と信頼することにして， ですね …… これで変動にびんかんな指数関数のくびきを逃れて， しゃかしゃか MCMC 計算が進行するようになればよいのだが， さて． ついでに呼吸速度も．
• この策はそこそこにアタって， 「ぼこぼこ山」の個数はへった． DIC = 1848． ちょーしに乗ってあつかいの難しそうな 初期勾配パラメーターも同様に 部分的脱指数関数化してみる? …… いきなりエラー吐くかと思いきや， 粛々と MCMC 計算が継続している． まち時間およそ 20 分間．
• 初期勾配の混交はよくなった (たぶん前よりは振幅そのものは小さいような気がする)． Gibbs sampler ってうまく作れるもんなんだなぁ． DIC = 1848． 最大光合成速度のまわりの混交は悪くなった． このあたり， 「個体差」も exp(...) から解除してみるか? この変更は解くに問題あるわけではないよな …… 混交がより劣悪になるかもしれんけど．
• 最大光合成速度・呼吸の混交はよくなった． DIC = 1843． 初期勾配にチェイン間の差がある． また初期値問題か? いや， これは意味がないよな …… 初期勾配の「個体差」事前分布の variance， 当然ながら極小になっている． これって何か関係あるのか? ためしに明るさの単位をかえてみる． 超事前分布のカタチが影響している， とは考えたくないが ……
• 結果は少し変わった． 初期勾配の混交はマシになる． DIC = 1843． そしてあとは burn-in (直訳は焼き付き，だけど計算上は「焼き捨て」みたいなかんぢ?) とサンプリングの期間を増やせばなんとかなるような 気がしてきた．

[計算だんだん改善]

あちこちに「個体差」が入っている状況でなお， 「処理の効果」 とやらが見えてきましたよ， ということをあらわす図．

• 20000 MCMC step 計算の命じておいて撤退． 昨日・本日の教訓は，
  • (大域的には) Gibbs sampler は劣悪な尤度地形から自力脱出できないので (パラレルテンパリングとかやらないかぎりは)， 初期値の与えかたに気をつける
  • (局所的には) しかしながら， 関数型の工夫とかで パラメーターの範囲を制約するよりも， むしろ Gibbs sampler を猟犬のごとく自由に動きまわらせたほうが， よく混交していて適切な範囲をもつ事後分布が得られたりする
  • つまり Gibbs sampler にとって働きやすい環境を整備してやる
  といったところか． 2140 研究室発． 2200 帰宅． 晩飯．
• いま使ってる 2 パラメーター (呼吸もいれると 3 パラメーター) 光合成曲線にちょっとヘンなところがある， と気づいた． 院生まんが PhD のセシリアほど 深刻な状況 ではないけど． 現在， 「初期勾配」と称しているパラメーターの正体は， 実際のところ 「初期勾配」/「最大光合成速度」 なる値になっているな …… どうりでヘンだと思った．
• [今日の運動]
  • すみません， 今日はシューズも持ってきたんですが， 計算にかまけていて走るのを忘れてました ……
• [今日の食卓]
  • 朝 (0730): 食パン．
  • 昼 (1350): 研究室お茶部屋． 食パン．
  • 晩 (2300): 米麦 0.7 合． ニラ・ショウガ麻婆豆腐．

2006 年 11 月 18 日 (土)

• 0740 起床． 朝飯． コーヒー． 0915 自宅発． 曇． 0930 研究室着．
• 昨晩の帰りがけに命じておいた 20000 MCMC step 計算はおよそ 2800 秒で終了していたか． そしてトリプルチェインがよくまざっている． たいへんよろしい． DIC = 1834．
• 昨晩帰宅後に気づいた f = 「初期勾配」/「最大光合成速度」 問題を修正するべく， 関数型を Pmax[j] * (1 - exp(-f[j] / Pmax[j] * PARi[i])) と変更して再計算を命じておく． これでたぶん PARi[i] の単位を変換といったことは不要になるんではないかな?
• 計算まち時間は作図プログラムかきとか．
• 関数型変更した計算， 1400 秒ではなく 1980 秒かかって終了． 予測より長くかかったのは上記のわずかな変更によって モデルが格段に複雑になったため． そしてこの割算導入による複雑怪奇さ増大の 副作用・悪影響のあおりをうけて収束がまったくダメダメになった …… ということで， 「何でも割り算するな!!」 教義は尤度評価における統計モデルの関数型にも応用可能とわかってしまった．
• 関数型をもとにもどす． ただしパラメーターの「名前」を「初期勾配とやらではありません」 を明示するために Pmax[j] * (1 - exp(-fcff[j] * PARi[i])) とする． 「初期勾配」なるものを示す必要があるなら， MCMC サンプル値を組み合わせて計算すればよい． また再計算を命じる．
• 計算まち時間に井田君と関数型の議論． ついでに現在えられつつある計算結果の説明．
• その間に再計算のほうは 25 分ほどで終了． DIC = 1828． しかし n.burnin = 6000, n.iter = 10000 ではまだまだ収束してないな． やはり 50 分ちかくの計算時間は必要とされているわけか．
• で， 再計算やらせてみたんだけど …… またまた初期値問題でしくじる． われながらこりないというか．
• 再計算命じておいて昼飯のため自宅にもどる． 昼飯． 洗濯． 1420 研究室にもどる．
• 20000 MCMC step 計算， 2800 秒後に終了してた． DIC = 1793． しかし絶対値の小さいパラメーターは収束悪いな． 計算機の精度の問題， adaptive rejection sampler のアルゴリズムとかと関係ありそう．
• ついでにあの (計算が不必要に複雑になる) 割算くみこみモデルも良い初期値を与えなおして再計算させてみる． こいつは時間かかるので 10000 step で終了させることに．
• 計算まち時間はまた作図プログラムかきとか．
• 1900 秒でさきほどの計算終了． DIC = 1828． やっぱり収束してないな． ということで， またまた明るさの単位変換して計算せてみる．
• 同じく 1900 秒で終了． こんどは最大光合成まわりの収束がよろしくない． DIC = 1764， こういう値だけは改善されているんだが． ぼこぼこ三峰は三本のチェインがからんでいないことをあらわす．

• 最大光合成速度が二箇所に影響あたえたりしてるんで， やはり割算モデルはダメか? ということでもとにもどして再計算． えーと 10000 MCMC step なので， 1400 秒ぐらいで終了するはず …… 終わった． DIC = 1778． やはり 10000 ではだめだな． そして初期値の選びかたもびみょー．
• で， 2800 秒かけてやりなおしてみると今朝えられた結果と同じく DIC = 1835． 意外と安定してるな． いまいち収束がよくないんだけど．

• ためしにもっと長くしてみるか …… ということで 4200 秒計算を命じてみる．
• 本日も計算まち時間 (20-50 分) は A 棟 8F 内をうろうろして， R わんぽいんと講習とか． 大学院生はしばしば難しいお客さんとなるもので ……
  • 自分がどこまでわかっているかよくわからない， あるいはわからないけどそれ自体を隠匿しようとする
  • ココロザシ高く勉強してるかと思えば， わりと安易に姑息な方向に逃げこんだりすることも
  • 言われたことは何も考えずに実行する， あるいは逆にぜんぜんちがうモノにすりかえてしまう
  …… まあ， こういう行動じたいは他ならぬこのぎょーむ日誌に記録されている 私自身の行動でもあったりするわけだが． 帰結がこうなりがちであったとしても， そこに至る経緯の「個体差」はたいへん大きく， また個体ごとにその「個体差」が定着しているわけでもない …… ということで， つねに注意ぶかい観察が必要とされる． なかなかうまくいかないけど．
• 4220 秒． DIC = 1824. さすがに 30000 MCMC step (最後の 4000 MCMC step) だとぼこぼこしないよな， と納得しそうになるところにワナがあるわけで …… ♪ 誰も知らない宇宙には 誰も知らない道がある ……

• 準備しておいた 「個体・季節ごとのあてはまり」 をチェックしてみる R の作図プログラムにながしこむと …… 7 個体中 2-3 個体の 1 シーズンが「びみょー」におかしい．
• その図の印刷出力を片手に A 棟 8F をうろうろ． お茶部部屋の片すみで曜日に関係なく展開されている院生たちの 「食」情報の高速やりとりを拝聴しつつ図をながめているうちに …… 発注者・井田君のしかけた「ワナ」にはまっていることに気づいた． いやいや， ワナというのはひどくて， いくつかの個体の「状態」が季節ごとに変化するという説明は きちんとうけていて， 自分のアタマにも入っていたんだけど， この推定結果はおそらく …… と部屋にもどってみると， そのワナ回避処理がなぜかうっかり脱落していることに気づいた． いやはやー
• このあたり一部の個体にしか関係ないわけだが， これが全個体の推定に影響およぼしてるのはまちがいあるまい …… ということで， (もう帰宅しようとおもっていたんだけど) 修正してとりあえず 8000 MCMC step 再計算． まち時間は 20 分ちょいぐらいか? また A 棟 8F R 苦闘院生の見舞．
• …… やっぱり， というか MCMC 計算の収束は改善されてしまった． いままでの 8000 step 計算とは別ものだね， これは． とりあえずあてはまりぐあいチェック作図． これでヨシ．

[個体差考慮した光合成曲線推定]

階層ベイズモデルの本領発揮 …… と言いたいところだが， 今回のモデリングでわかったのは， 私はこのあたりの現象についてまだまだ理解が足りん， ということ． 今回の問題をすごく端的にいえば「葉の老化」をモデリングせよ， ということなんだけど， 推定やったり図を描いたりしてるうちに 今つかっている統計モデルとはまた別のモデルがアタマにうかんできた． それは「葉の老化曲線」といったものがあり， 個々の葉によってかなりカタチが異なっている． 早めに老化するものもあれば， 他よりも老化が遅れるものもある …… といった描像だ． さらに， 単なる憶測なんだけど， その老化曲線なるものはどこまでもはてしなく下降するものではなく， たとえば最大光合成速度に関して 「最低の最大光合成速度」 みたいなものがあり， 「自然状態で葉として機能しているならば， これよりも劣悪な光合成性能にはならない」 といった水準が存在するような気もする (たとえば葉内の窒素を全部ひきぬくことは不可能だから， といった理由で) ……

• 20000 MCMC step 計算を命じておいて撤退 …… のつもりがまたデータ解析こんさるやってしまって， 2140 研究室発． 2200 帰宅． 晩飯．
• [今日の運動]
  • 明日はうんどうを!
• [今日の食卓]
  • 朝 (0810): 食パン．
  • 昼 (1330): 米麦 0.7 合． ニラ・ショウガ麻婆豆腐．
  • 晩 (2240): スパゲッティー． ニラ・ショウガ麻婆豆腐．

2006 年 11 月 19 日 (日)

• 0840 起床． 朝飯． コーヒー． 昨晩かえりがけに命じておいた計算結果が気になるんだが …… まあうまくいってるだろうと楽観することに． 洗濯． 怠業．
• 1200 自宅発北大構内走． 晴． 今日は気温たかめ． 走りつつ植物個体内資源分配モデル， とか「次の下うけ計算」になりそうな問題について検討してみる． 重量で観測される個体内資源分配はわりとすなおなモデリングになるはず． 問題は重量ではなくカタチで表現される部分だよな …… と考えてるうちに， なぜか査読報告かきの作文を検討していたり． 1250 帰宅． 体重 72.6kg． 昼飯．
• 北大ネットのへっぽこな spam filter が「破られた」気がする． いやいや， いつものごとく単なる故障かも．
• 1500 自宅発． 曇． 1515 研究室着．
• 昨晩かえるまえに命じておいた計算 …… うーむ， 処理とかのパラメーターがえらくきれいに収束してるぢゃん， と思ったら， 帰宅直前にどたばたやったプログラム修正にばぐがあって， このパラメーターはデータのくびきから開放されて MCMC ならぬ単なるモンテカルサンプリング (つまり事前分布 == 事後分布) になってしまっているじゃないか ……
• 脱力． 修正． 再計算．
• 1400 秒で終了． DIC = 1837． 収束は 10000 MCMC step ならこんなもんか． 2800 秒再計算を命じる．
• これも順当に終了して DIC = 1833． 収束は下の図のようなかんぢで， まあこのあたりでよいかな， と．

• 特にアイデアもないので， 4200 秒計算を命じてみる …… 計算結果， あまりかわったかんぢがしない． ということで， この計算はいったん終了． A801 の Dell 機からは撤収．
• 今日の夕方はずっと推定結果とかをながめて改善策を考えていたんだけど， なかなか思いつかず …… 「処理の効果のかかりかた」 に個体差があるように見えるんだけど， こいつをどう表現すればよいのかよくわからない． その理由は 「季節ごとの最大光合成の個体差」 と区別がつかなくなるようにおもえるので．
• まあ， バックアップとかはすでに取ってるんで …… 実験的に 「処理の効果のかかりかた」 の個体差とやらをいれてみるか． これがうまくいかなかったら， 今までの計算結果を使うことにしよう． どうなることやら． とりあえず 2000 + 4000 MCMC step で計算させてみる．
• うーむ， マシな結果が得られてるような気がする (無駄パラメーターが 4×2 個ふえたので DIC 的にはあまりよくないけど)． じゃあまた 20000 MCMC step 計算やらせてみるか． 計算時間は一時間ぐらいかしらん?
• 晩飯調達にでる． 研究室お茶部屋にもどって晩飯．
• 2800 秒で終わった． やっぱり収束は悪い． しかし推定された曲線のカタチはよい …… ということで定式化にもうひと工夫ひつようか． 無駄パラメーターをなくせということだろ， これは．
• 本日の「計算まち時間副業」は小山さんとサロベツデータ解析のハナシ． これまた例によって例のごとくややこしいデータ構造なので， 「ベイズ統計やりましょうよー」 と例によって例のごとくベイズ勧誘．
• うーむ， 光合成モデルのほうに「処理への応答の個体差」 をくみこむのは難しい， とわかった． ということで， 実験用プログラムは廃棄． その前に保存してあったやつを明日あたり納品してしまおう．
• 2245 研究室発． 2300 帰宅．
• [今日の運動]
  • 北大構内走 50 分間．
• [今日の食卓]
  • 朝 (0930): 米麦 0.6 合． ネギ・ブナシメジ・ワカメ・卵のスープ．
  • 昼 (1300): 食パン． ネギ・ブナシメジ・ワカメ・卵のスープ．
  • 晩 (2010): 研究室お茶部屋． ほっかほっか亭， 牛とじ重． 450 円． 牛肉くったのはずいぶんとひさしぶりな気がする．

2006 年 11 月 20 日 (月)

• 0825 起床． 朝飯． コーヒー． 0925 自宅発． 晴． 0940 研究室着． せみなー予習．
• 研究室セミナー， 本日は Annika さんで極地のコケの OTC & 施肥実験． point-grid 法でやってるんだから， カウントデータになってるわけで， だったらいろいろモデリングは簡単なんだが …… と思いつつ (いわゆる) ANOVA に「流しこまれて」 ややこしくなってしまった結果をながめる．
• 井田君にここまでの階層ベイズ光合成モデルのあれこれをひきわたし． モデルの説明をまたくりかえしてみたんだが …… やっぱり自分でも いまいちすっきりしないモデリングだと思う．
• 昼飯調達のため北大生協にいき， ついでに書籍部とかもうろうろと． 研究室にもどって昼飯．
• そして夜まで小山さんとサロベツデータ解析．
• また「やちパワー地図」 (例: 以前つくったもの) をつくりなおす． 今回は今年のデータつかったカウントデータ版． どの植物種がどこに何個体あらわれてきたか， をあらわす地図．
• そしてぢりぢりと階層ベイズモデル説明を． とりあえずは一般化線形混合モデル (GLMM) でも解けるような簡単なところから出発． というか， 小山さんの場合は「GLMM は省略」みたいなかんぢですな．
• で， 小山さん PC に WinBUGS と R の R2WinBUGS package をインストールして， bugs ファイルなど書いていく．
• …… しかし， 本日は WinBUGS の解読不能な「Trap 窓」出現にやられたまま終了． つづきは明日， ということで． どうしてこんなに簡単なモデルがうまく計算できないのかな?
• 2030 研究室発． 2050 帰宅． 晩飯． ばてた．
• [今日の運動]
  • やっぱ運動するなら「昼」だよなあ …… このしーづんは
• [今日の食卓]
  • 朝 (0840): 米麦 0.6 合． ワカメスープ．
  • 昼 (1340): 研究室お茶部屋． 食パン．
  • 晩 (2130): スパゲッティー． レトルトパウチドカレー． ワカメの酢のもの．