WinBUGS 入門 水産資源学におけるベイズ統計の応用ワークショップ 2007 年 8 月 2-3 日, 中央水研 遠洋水産研究所外洋資源部 鯨類管理研究室 岡村寛
WinBUGS とは BUGS (Bayesian Inference Using Gibbs Sampling) の Windows バージョン フリーのソフトウェア Gibbs samplingを利用した事後確率からのサンプリングを行う
WinBUGS プログラムコードは R とよく似ている それゆえ,Rでプログラムが書けると理解が早いし, 何かと便利 R2WinBUGSというRからWinBUGSを呼び出して実行するプログラムも開発されている 様々な状況に自動で対応して適当なサンプリング法を選んでくれる 例. log-concave adaptive rejection sampling 範囲を制限したとき スライスサンプラー範囲制約がないとき Metropolis 法
WinBUGS 画面
WinBUGS プログラム 基本的に,model.odc,data.dat,inits.inの3つのファイルからなる. model.odc の中身の例 model{ for (i in 1:N){ Y[i] ~ dnorm(mu[i],tau) mu[i] <- a + b*x[i] # a + b*(x[i]-mean(x)) is better. } }
WinBUGS プログラミングの注意 1 変数を2 度定義できない x <- 1 x <- 3 ベクトルや行列で対応 if 文は使えない X = (3, NA, 2) X=cbind(c(1,3),c(3,2)) X[t,2]~dpois(N[X[t,1]]) ( ただし, この場合はNAのままでもOK) Y ~ dnorm(log(n), sigma) Y ~ dnorm(ln, sigma) LN <- log(n)
WinBUGS プログラミングの注意 2 最初に少ないパラメータからスタートして ( いくつかのパラメータをfix), 段々複雑にする R2WinBUGSを使うとデータの扱いなど便利 bugs(data, inits,, debug=true, ) Trapというエラーメッセージが出ることがある 初期値がおかしい事前分布を少し狭く ~ 徐々に広げる
WinBUGS プログラミングの注意 3 まず最尤法を用いた簡単なプログラムであたりをつける プロット! プロット! プロット! あまり洗練されたプログラムを書こうとしないこと よく使う人と友達になる 仲間を増やす
例題 1( 検定 ) 池 A, 池 B の中のある場所をランダムに選んで毎回同じ時間調査し, その中にいるメダカの数を数えた結果, 下のようなデータが得られた ( 池 B は行きにくいので 3 回しか調査できなかった ). 池 A と池 B のメダカの数は違うであろうか?
検定仮説 検定したい仮説 : 池 Aの平均メダカ数 = 池 Bの平均メダカ数 池 Aと池 Bの平均メダカ数は, 池 A: 平均 1.5, 分散 1.7 池 B: 平均 0.0, 分散 0.0
t 検定 池 B はメダカのサンプルなしなので, 池 A の平均 が 0 と有意に違うかどうかを見る t 検定を行う > x <- c(1,1,0,4,2,2,0,2) > t.test(x) t = 3.2404, df = 7, p-value = 0.01425 有意
通常の検定の何が問題か? 池 Bはすべて0の観測値であるが,3 回 0が続いただけである. これは10 回 0が続いたのとは証拠の大きさが異なるはずである. 先の検定は, そこのとこを考慮していない. 3 回とも0だったので分散は0となり,0 確実である. しかし,3 回しかやってないので, たまたま0が続いたということもあるのでは? > rpois(3,1.5) [1] 0 0 0 なんていこうともありうる.
GLM してみる カウントデータであるし, 池 A の平均と分散はおよそ等しいのでポアソン分布を仮定して GLM してみる a <- factor(c(rep("a",8),rep("b",3))) x1 <- data.frame(area=a, med=c(x,rep(0,3))) res <- glm(med~area,family=poisson,data=x1)
GLM 結果 > summary(res) Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) 0.4055 0.2887 1.405 0.160 areab -19.7081 5442.46-0.004 0.997 収束していない. これは池 B が全部 0 だからである
こんなふうに考えてみる 池 Aから3 回続けてサンプルしたとき, すべて0となる確率はどうなるか? これがすごく小さかったら, 池 Bは池 Aよりメダカが少ないと考えてもおかしくないであろう. 確率が大きければ, 池 Aと池 Bが違うと考える根拠は希薄であることになる.
事後予測分布 ベイズ推定で考えると, 池 A のデータが与えられたとき, 新たな調査を行って 3 回続けて 0 となる確率を推定す れば良い y = { 新たに調査したら3 回続けて0} とするとき, 事後予測分布 P(y x) = P(y θ)p(θ x)dθを推定する ポアソン分布を仮定すると3 回続けて0 = {exp(-θ)} 3 = exp(-3θ)
WinBUGS コード model{ Jm ~ dunif(0,100) mu <- pow(jm,-2) for (i in 1:N){ x[i] ~ dpois(mu) # likelihood } Prob <- exp(-3*mu) } list(x = c(1,1,0,4,2,2,0,2), N=8) # data
WinBUGS 実行 1. プログラムを WinBUGS 上で開く 2. ツールバーの Model の Specification を選択 3. プログラムの model 部分を選んで, Specification Tool の check model をクリック 4. データ先頭の list を選んで,Specification Tool の load data をクリック 5. compile をクリック
WinBUGS 実行 6. 初期値を自動生成するため, 今回は gen inits をクリック 7. ツールバー Inference をクリックして, Samples を選択 8. Sample Monitor Tool の node に mu と入れて, set をクリック 9. Sample Monitor Tool の node に Prob と入れて,set をクリック
WinBUGS 実行 10. 再びツールバーの Model をクリックし, 今度は Update を選択 11. Update Tool の update ボタンをクリック 12. Sample Monitor Tool に行って,node に mu と入れる. 下のボタンが active になる 13. density や stats ボタンを押して結果を見る 14. node に Prob を入れて同じことをする
WinBUGS 実行 と文章や口で言っても分かりにくいので, 実際にやってみる 実演 WinBUGS the movie! というページもある. 動画で紹介. http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/winbugsthemovie.html
事後分布 平均 1.43 95%CI[0.72, 2.38] 平均 0.03 95%CI[0.001, 0.11]
まとめ ベイズ推定によって, うまいこと問題 ( 池 Bのサンプルサイズを考慮, 観測したどのデータも0, など ) を回避していることに注意 池 Bの結果を, 池 Aの結果が与えられたという条件のもとで, 池 A= 池 Bとしたとき, 池 Bの結果と同じようになる確率で評価した. 他にどんなモデリングが考えられるか?( 例. 池 Bの方も池 Aと同じようにモデル化して, 二つの平均を比較するとか )
例題 2( 階層ベイズ ) 大西洋ビンナガマグロの漁獲量と CPUEのデータにプロダクションモデルをフィットし, そのパラメータの推定結果を用いて,MSYを計算せよ. 最近の漁獲量は適切な水準にあると言えるか?
CPUE (Catch Per Unit Effort) CPUE = 漁獲量 / 努力量 CPUE = qp q: 漁獲効率,P: 資源量 ( 個体数 or 個体重量 ) CPUE を使って個体群の増減が分かる 年トレンドの推定 CPUE の標準化 (GLM 等 ) CPUE = qp x (Hilborn & Walters 1992)
MSY の考え方 個体群の増加率 ( ロジスティック曲線 ) dp/dt = rp(1-p/k) P: 資源量,r: 内的自然増加率, K: 環境収容量 MSY(Maximum Sustainable Yield): dp/dt が最大になるのは,P=K/2 のとき 0 K P P=K/2 になるように魚を獲っていたら, 最大の収穫を得ながら, 個体群を健全に維持できる!
状態空間モデル State-Space Model X t-1 X t X t+1 Y t-1 Y t Y t+1 プロセス方程式 : X t = G(X t-1, w t ) 観測方程式 : Y t = F(X t, v t )
個体群動態モデルのパラメータ推定 ロジスティックモデル B t+1 = {B t + rb t (1-B t /K)-C t }exp(u t ) I t = qb t exp(w t ) MSY = rk/4 B MSY = K/2 E MSY = r/(2q) u t ~ N(0,s A2 ), w t ~ N(0,s O2 ) プログラム中では,r と K の相関を弱めるため,P t = B t /K を使う
状態空間モデルの特徴 観測誤差とプロセス誤差を同時に推定することが可能 計算機集約的な方法と組み合わせることにより, 柔軟なモデリング 推定が可能 カルマン フィルタは状態空間モデルの特殊な場合と考えられる 一般に計算大変
データ : 漁獲量と CPUE 南大西洋ビン ナガマグロの データ
WinBUGS 実行 実演
事後分布 110,000 万回計算して, 最初の 10,000 回を burn-in として取り除いた
R と K の相関 K 600.0 500.0 400.0 相関係数 = -0.84 300.0 200.0 0.0 0.2 0.4 r
自己相関 thinning なし thinning=50 thinning=100 thinning=200
2 本の chain でやってみる V R W R = V/W V = pool した chain の変動 W = 個々の chain の平均変動
2 本の chain でやってみる
R2WinBUGS RからWinBUGSを呼び出して解析を実行し,Rの中でoutputを操作 データや初期値の作成が楽ちん 即 codaを使える 最初は少ない繰り返し数で
R2WinBUGS でやってみる
coda
まとめ 個体群動態モデルでは, 従来観測誤差かプロセス誤差のどちらか一方だけを仮定することが多かったが, ベイズ型状態空間モデルを使うことによってどちらの誤差も考慮した分析が比較的容易に行える いろいろプロットして, うまくいっているか確認することが大事 R2WinBUGSを使うといろいろ便利なことが
参考文献 Meyer, R., and Millar, R. B. 1999. BUGS in Bayesian stock assessments. CJFAS 56: 1078-1086. Plummer, M. 2006. CODA: convergence diagnosis and output analysis for MCMC. R News 6(1): 7-11. Spiegelhalter et al. 2003. WinBUGS user manual version 1.4.