ぎょーむ日誌 2004-04-18

2004 年 04 月 18 日 (日)

• 0900 起床． 朝飯． コーヒー．
• 昨日， 北大構内走しながら， いまやってる苫小牧直径モデル (ややむりやり離散化) で使ってる確率分布すなわち Poisson 分布のいくつかの特徴について考えていた:
  • 異なるパラメーターをもつ Poisson 乱数たち (それぞれ独立) の和もやはり Poisson 乱数である (再生性; reproductive property)
  • Poisson 分布のパラメーターが Gamma 乱数であるとき (Poisson 分布の混合)， 分布は負の二項分布になる
  …… これらよく知られたふたつの基本的な性質を組み合わせて 悪用すれば， 現在の問題解決につながるのではなかろーか， という手がかりを得た．
• すなわち， 成長のよしあしが樹木個体ごとによって異なる (個体は自分のおかれてる「個別に特殊な境遇を記憶」しつづける; かつ， これが人間には観測されてない要因である)， という状況を混合 & 再生 Poisson 分布， つまりごちゃごちゃとややこしいパラメーターの積にて表現される 負の二項分布でモデリング可能では， という奇策だ． しかも， ごくまれに生じるマイナス成長もうやむやのうちに「帳消し」にされる， というオマケつき． こういうモデリング & 推定計算技法が通用するのかしないのか …… 職業的うそつきたる私の欺瞞能力が問われるところだ．

「バレなきゃイカサマじゃあねえんだぜ……」 (空条承太郎, 1989)

• 一番重要なのは， 自分で創作したうそに自分自身がだまされないことである． 次に重要なのは， 他人をだます手間を最小化することだ．
• 分布の平均に関しては首尾一貫してるので問題なさそうだ． 分散に関しては (1 + 1/c) 倍になる (overdispersion を表現)． この c ってのは R ふう にいうと Gamma 分布の rate == 1/scale のことで， 私がよくもちいる α & β notation だと 1/β のことだ …… これも極限をとると， つまり Gamma 分布の分散をゼロにもっていくと ふつーの Poisson 分布に「もどる」ので (あたりまえか)， 問題ないように見える． ということで， このように混合してしまって基本的に不都合なかろう．
• とすると， 次は推定方法の手順か …… のろのろと考えてるうちに難点に気づいてしまいましたよ． えーい． じつは安直に glm.nb(MASS) あたりで始末できるんでは， と楽観していたんだが …… Poisson 乱数の和が Poisson 乱数になる， という再生性なる性質を悪用するためには， 年ごとに特定される パラメーターの和が必要になっちまうよな． 積なんかではないぞ． マズいねえ．
• ちょっと打つ手ナシ． 1250 北大構内走にでかける． 晴． 1340 帰宅． 体重 74.0kg． 昼飯．
• で， 走りながら検討してみたんだけど， 上述のモデルは
  • 基本アイデアは問題なさそう; パラメーター推定可能
  • ただし一般化線形モデルのたぐいには帰着できない; 安直計算に逃げこめない
  ということが理解できた． 結論として， この考えかたは数式とその演算を簡単にするけれど， やはり optim() など力わざでごりごりと推定計算 (尤度最大化) するところは避けられない …… まあ， いいか． この直径成長モデリングの難所たる 「観測してない個体差」&「測定誤差」 どもを一撃で同時にツブせた (あるいはじゅーたんの下に隠しこめた) わけだし． 定式化に関してはとうとう決着ついた， とゆーことで．
• 1520 自宅発． 晴． 1535 研究室着．
• 数式など紙にうだうだと書いてみて検討する …… しかし Poisson 分布のパラメーター (すなわちここでは成長量の平均値) が Gamma 分布にしたがうってのも， なかなかにとりっきーな捉えかたかも． つじつまはあってるけどね． Poisson * Gamma → 負の二項分布， などという関係に気づいた奴って， そうとうにヘンな注意力の持ち主だったんではなかろーか ……
• 進捗しない． すでに夜．
• さて， 実際の推定計算の手順なんだけど …… どうも今までみたいに data.frame() で成長データぽんと渡して， みたいにはいかないような気がするのですよ． 私は． 格納するデータ構造を定義すべきか． それとも data.frame() でなんとかすべきなのか．
• よーわからんのでこのへんに書き散らしてみますかね (in R ふう疑似コード) ……
  • パラメーター値がわたされる
  • 全個体に共通して適用: 定数項， month(?), 気象値項

    lp.common <- function(...)
    	(1, month, (weather terms ...)) * p.common
    

  • 各個体ごとに適用: サイズ依存項， サイズ-気象値「交互作用」項

    lp.size <- function(tree, ...)
    	(tree$dbh.log, tree$dbh.log.sq, (d-w cross terms ...)) * p.size
    

  • 個体の平均成長量 (数年ぶん) = Σyears exp(線形予測子)

    mu <- function(...)
    	sum(
    		sapply(
    			years,
    			function(y)
    				exp(sum(lp.common$y) + sum(lp.size$y))
    		)
    	)
    

  • 全個体の対数尤度の合計 (負の二項分布)

    log.likelihood <- sum(
    	sapply(
    		trees,
    		function(t)
    			log(dnbinom(t$growth, size = theta, mu = mu(t, ...)))
    	)
    )
    

• うう， ややこしー …… どうすれば計算コストが減らせるのかわからん． いっそ R やめるか?
• さらにこまごまとした問題とかがあるし． たとえば optim() とのやりとりだな． c(1:4)[c(0, 1, 0, 1) %in% 1] とすると 2 4 と答える． このへんの性質を使って mapping， とかせんといかんのでしょうなぁ．
• かなり整理整頓してやったというのに， それでもなおとうてい簡単には片づいてくれそうにもない， なる見とおし得られただけで 本日はすごすごと撤退． 2045 研究室発． 2055 帰宅． 晩飯．
• [今日の素読]
  • Salsburg, D. 2001. `` The Lady Tasting Tea -- How statistics revolutionized science in the twentieth century''. Owl Book.
  • Chapter 25. Advice from the Lady in Black
    • Statistics at Guinness
    • Unexpected variability
  • The list of types of collabiration is typical of the work of a statistician in industry. In my own expericence, I have had interacations with chemists, pharmacologists, toxicologists, economists, clinicians, and management (for whom we developed operations research models for decision making). This is one of the things that make the work day of a statistician fascinating. The methods of mathematical statistics are ubiquitous, and the statistician, as the expert in mathematical modeling, is able to collaborate in almost every area in activity.
• [今日の運動]
  • 北大構内走 1250-1340． ストレッチング．
• [今日の食卓]
  • 朝 (0920): パンケイキ．
  • 昼 (1410): スパゲッティー． S&B バジリコソースなる 面妖なる半固形物．
  • 晩 (2140): 米麦 0.9 合． チンゲンサイ・モヤシ・ピーマン・ショウガ・エノキダケの炒めもの．