ぎょーむ日誌 2006-01-06

2006 年 01 月 06 日 (金)

• 0750 起床． ねぶそく． ねむい． コーヒー． 朝飯．
• 寝る前に始動した 50 MCMC step とばしの試験運転的計算， ひととーり終わってるな． うーむ， なかんぢだけど記録のためここにのせておく． 事前分布の逆分散を決める hyperprior (ガンマ分布てぬきの指数分布) のばらつきは scale = 3.0 (つまり IG(1, 3))． 事後分布描画の標本範囲は 501:1500．

シュート伸長の休眠・二度伸び

葉面積あたり窒素量モデル

葉重量あたり窒素量モデル

• 今朝もいろいろ考えてみたんだけど， 事前分布 (prior) と超事前分布 (hyperprior) まわりについてもうちょいよく考えてみる必要あるような気がする．
• なおこのあたりのモデリングでは 「階層」というコトバがむやみに使われているんだけど， ちょっと整理してみると
  • 階層ベイズモデルの意味での階層: hyperprior - prior (- parameter)
  • 生物 (樹木) の階層 (むしろ nested の意味): hyperspecies - species - tree
  となる． 両方について考える．
• まず prior に関してはどれも 「平均ゼロの正規分布」 としており， この正規分布の分散だけ (正確には分散の逆数だけ) が hyperprior で与えられていて， これがガンマ分布てぬきの指数分布になっている．
• ここで， 平均ゼロの正規分布を事前分布とする理由について， 生物の階層のがわから考えてみる．
• Tree の層においてはこれは明確で， 「個体差」をあらわす random effects である， という説明になる．
• Speceis の層は個体の層より難解かも． つまり「樹種差」なるものが random effects である， としているわけだ． この屋久島常緑樹の観測データ一式と研究目的からすると， この設定でよいように思う． 一樹種あたりのデータが少なく， さらに種間をこえた共通点を探索する研究であるからだ．
• さて問題は hyperspecies の層である． ここで事前分布を「平均ゼロの正規分布」とか置いてしまってよいのか? これはよくわからない． ここでは事象 (休眠・二度伸び・葉っぱの死亡) が発生する確率を logit のたぐいのモデルで表現しているので， 線形予測子 == ゼロは確率 0.5 つまり 「発生するともいえるし発生しないともいえる」 状況と対応している (このあたりの説明やや強引)．
• これが logit モデルでなければ (たとえばポアソン回帰みたいなことをやっているのであれば)， 事前分布が平均ゼロの確率分布， というのは明らかにまずいんじゃないの， という場合が多いだろう． 少なくとも定数部分 (いわゆる「切片」) はゼロではあるまい ……
• Hyperspecies の層における事前分布とパラメーターの関係は， やや奇妙なものに見える． hyperspecies 用のパラメーターは一個だけ生成されるんだけど， これを生成する事前分布は 超パラメーターを一個または二個必要とする． パラメーターが生成され， 事前分布はそれを追いかけ， またそこからパラメーターが生成され， 事前分布はまたそれを追いかける …… となるわけだよな． まあ， species と tree の階層においては， 独立同分布の事前分布から種数ぶん・個体数ぶんの パラメーターが生成されるんだけど， 事前分布が自分の生成した「点々」をおっかけていく， というのは同じか．
• さてどうすりゃいいのかよくわからぬまま， 1020 自宅発． 今日も雪． 1035 研究室着．
• 屋久島シュート伸長休眠の推定計算， hyperprior で scale = 1 とかダメっぽいような気がする． 個体差とかが無意味にでかくなる． scale = 20 とかにしてみる (上で示してる計算結果では scale = 3)．
• また事前分布についてうだうだと考えてみたり …… ふと研究室ネット雑用一個を放置してたのに気づいたんで (昨日のさーばークラッシュで中断されたアカウント作成)， さっさと片づける．
• ここ数日， 屋久島まわりの問題にいわば「かまけてる」状況なんだが …… 冷静に自分のおかれた状況をみると危機の包囲網 ぢりぢりとせばまりつつあるような気がする． 整理してみよう．
  • 苫小牧樹木直径成長モデルの原稿 …… ともかく苫小牧ボスと連絡をとらなくては (そして私の担当すべき部分を書く，と)
  • 生態学会新潟大会要旨 …… 11 日 (水) が締め切りなので， さっさと文案つくって小林さん・加藤さんに連絡
  • 19 日 (木) の東京出張準備． はあ．
  • 23 日 (月) と 25 日 (水) の M1 授業準備．
  これだけでなく他の修論したうけとか塩寺さん計算したうけとか ……
• 1240 ごろ hyperprior で scale = 20 とした休眠・二度伸び計算やりなおしが終了． 幸か不幸か， 超事前分布のパラメーターを変えたぐらいでは結果はほとんど同じ． 「多峰」状態は多少マシか? マシになったのがいいことなのかマズいのかよくわからんけど． よくわからぬまま， この結果を使って葉っぱ生き死に確率の推定計算もやらせる． あいかわらず試験運転的な 50 MCMC step とばしのサンプリング．
• 昨年 plotnet あぶく銭で買った (買わされた) ArcGIS 9.0， これを 9.1 無償アップデイトうんぬんはいいとして， そのライセンスキーとやらを書いた紙きれがどこかにあるはず， と探し回ることに (私が使うわけでもないのに)． 意外にも， というかわかりにくいところながらきちんと保存していたので， なんとかなった． これだから有償ソフトウェアは ……
• 昼飯．
• 雑用少々． A802 のプリンター LBP-2810 のトナーリサイクルだしを 雪野さんにおねがいする． 黒と青． みなさん夜中にトナーをスプーンですくって食べてるんぢゃないかと思うほど， ちかごろは減りが速い． ついでに年末大掃除のときにでた備品番号つき蛍光灯スタンドのあつかい についても教えていただく． これは備品番号シールを紙にはって提出すればよいとのこと． 捨てる前に会計係に連絡， とか．
• 屋久島計算進行中． 超パラメーター変えたぐらいではあまり変わらんのかな． それはそれでたいへんけっこうではあるんだけど．
• 苫小牧樹木成長モデルの原稿， これ以上放置しとくのはすごく危険というかコワくなってきたので， いまのうちになんとかしなくては ……
• しかし， あまり内容のないメイルを苫小牧に送信することになってしまった． どうなることやら …… 時間はかかってしまってすでに 1610． 嗚呼．
• 事前分布の hyperprior で scale = 20 とした試験運転的計算 (つまり 50 MCMC stepとばしサンプリング …… MCMC バグをかたづけてからはこれぐらいの間隔で十分というかんぢ)， 葉重量窒素濃度モデルに関しては終了． まあ， あまりかわりばえしない． あいかわらず窒素が多いと長命化 …… ということで， 元データから次のような図を作ってみる．

[明るさと窒素の交絡]

重量あたり窒素濃度． R の density() で推定した経験確率密度分布． 明るいほうがびみょーに窒素が多かったり．

[あやしげ図]

LMA なるヘンてこな割り算値と (じつはこれまた割り算値にすぎない) 重量あたり窒素濃度の関係． まあ， 上と同じく明るさと交絡があるってことぐらいかな …… ということでタテ軸・ヨコ軸のあいだの関係はどうでもよい．

• 蛇足ながら， (なんと割り算値どうしの割り算値である) 「面積あたり窒素濃度量」 なる数量は もうめちゃくちゃというか， 事実上 LMA の逆数にすぎない， とわかった (上の図をみれば明らかであろう) …… C-N コーダーなる測定器の吐き出す数値って何かヘンなんじゃないか? まあ， よくわからんけど， いま計算させている二種類の葉っぱ生き死にモデルのうち 「葉面積窒素濃度モデル」 ってのは単に「LMA の逆数」 (この割り算値は SLA などと呼ばれていたりするんだが) を要因にしました， ぐらいの意味に考えたほうがいいな．

  後記: ここは記述がまちがっていて (さいわいにも計算プログラム内ではまちがっていなかったんだけど)， Narea = Nmass * LMA が正しい． つまり， Narea なる量は LMA とほぼ同じ， ということ．

• いやはや．
• 院生密度が高くなると， お茶部屋に拘束されてしまう時間が長くなってしまう．
• 夜になってから， 逆分散の超パラメーターを scale = 20 とした計算結果ひととーりでた． これは 50 MCMC step とばしのサンプリング． あまりかわりばえしないけれど， ここに記録しておく．

シュート伸長の休眠・二度伸び

葉面積あたり窒素量モデル

葉重量あたり窒素量モデル

• で， 本日から登校している計算発注者の牛原さんに上の図ひととーり説明．
• 事前分布のいれかたを変えたいところなんだが …… 本日はもう計算プログラム改変とか手をださないほうがよさそう． そこで単純作業としてサンプリング間隔を 100 MCMC step に倍化させた計算を命じておいて， 2110 研究室発． 2130 帰宅． 晩飯．
• [今日の運動]
  • 修論したうけ運動ふそくがつづく ……
• [今日の食卓]
  • 朝 (0830): サツマイモ電磁波蒸し．
  • 昼 (1310): 研究室お茶部屋． 米麦 0.5 合． ハクサイ・ニンジン・ブナシメジ・豆腐のシチュー．
  • 晩 (2200): 米麦 0.7 合． ネギ卵炒飯． 海藻スープ．