ぎょーむ日誌 2006-11-17

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 合． ニラ・ショウガ麻婆豆腐．