スライド 1

Size: px

Start display at page:

Download "スライド 1"

えのみょうだに
5 years ago
Views:

1 WinBUGS 入門水産資源学におけるベイズ統計の応用ワークショップ 2007 年 8 月 2-3 日, 中央水研遠洋水産研究所外洋資源部鯨類管理研究室岡村寛

2 WinBUGS とは BUGS (Bayesian Inference Using Gibbs Sampling) の Windows バージョンフリーのソフトウェア Gibbs samplingを利用した事後確率からのサンプリングを行う

3 WinBUGS プログラムコードは R とよく似ているそれゆえ,Rでプログラムが書けると理解が早いし, 何かと便利 R2WinBUGSというRからWinBUGSを呼び出して実行するプログラムも開発されている様々な状況に自動で対応して適当なサンプリング法を選んでくれる例. log-concave adaptive rejection sampling 範囲を制限したときスライスサンプラー範囲制約がないとき Metropolis 法

4 WinBUGS 画面

5 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. } }

6 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)

7 WinBUGS プログラミングの注意 2 最初に少ないパラメータからスタートして ( いくつかのパラメータをfix), 段々複雑にする R2WinBUGSを使うとデータの扱いなど便利 bugs(data, inits,, debug=true, ) Trapというエラーメッセージが出ることがある初期値がおかしい事前分布を少し狭く ~ 徐々に広げる

8 WinBUGS プログラミングの注意 3 まず最尤法を用いた簡単なプログラムであたりをつけるプロット! プロット! プロット! あまり洗練されたプログラムを書こうとしないことよく使う人と友達になる仲間を増やす

9 例題 1( 検定 ) 池 A, 池 B の中のある場所をランダムに選んで毎回同じ時間調査し, その中にいるメダカの数を数えた結果, 下のようなデータが得られた ( 池 B は行きにくいので 3 回しか調査できなかった ). 池 A と池 B のメダカの数は違うであろうか?

10 検定仮説検定したい仮説 : 池 Aの平均メダカ数 = 池 Bの平均メダカ数池 Aと池 Bの平均メダカ数は, 池 A: 平均 1.5, 分散 1.7 池 B: 平均 0.0, 分散 0.0

11 t 検定池 B はメダカのサンプルなしなので, 池 A の平均が 0 と有意に違うかどうかを見る t 検定を行う > x <- c(1,1,0,4,2,2,0,2) > t.test(x) t = , df = 7, p-value = 有意

12 通常の検定の何が問題か? 池 Bはすべて0の観測値であるが,3 回 0が続いただけである. これは10 回 0が続いたのとは証拠の大きさが異なるはずである. 先の検定は, そこのとこを考慮していない. 3 回とも0だったので分散は0となり,0 確実である. しかし,3 回しかやってないので, たまたま0が続いたということもあるのでは? > rpois(3,1.5) [1] なんていこうともありうる.

13 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)

14 GLM 結果 > summary(res) Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) areab 収束していない. これは池 B が全部 0 だからである

15 こんなふうに考えてみる池 Aから3 回続けてサンプルしたとき, すべて0となる確率はどうなるか? これがすごく小さかったら, 池 Bは池 Aよりメダカが少ないと考えてもおかしくないであろう. 確率が大きければ, 池 Aと池 Bが違うと考える根拠は希薄であることになる.

16 事後予測分布ベイズ推定で考えると, 池 A のデータが与えられたとき, 新たな調査を行って 3 回続けて 0 となる確率を推定すれば良い y = { 新たに調査したら3 回続けて0} とするとき, 事後予測分布 P(y x) = P(y θ)p(θ x)dθを推定するポアソン分布を仮定すると3 回続けて0 = {exp(-θ)} 3 = exp(-3θ)

17 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

18 WinBUGS 実行 1. プログラムを WinBUGS 上で開く 2. ツールバーの Model の Specification を選択 3. プログラムの model 部分を選んで, Specification Tool の check model をクリック 4. データ先頭の list を選んで,Specification Tool の load data をクリック 5. compile をクリック

19 WinBUGS 実行 6. 初期値を自動生成するため, 今回は gen inits をクリック 7. ツールバー Inference をクリックして, Samples を選択 8. Sample Monitor Tool の node に mu と入れて, set をクリック 9. Sample Monitor Tool の node に Prob と入れて,set をクリック

20 WinBUGS 実行 10. 再びツールバーの Model をクリックし, 今度は Update を選択 11. Update Tool の update ボタンをクリック 12. Sample Monitor Tool に行って,node に mu と入れる. 下のボタンが active になる 13. density や stats ボタンを押して結果を見る 14. node に Prob を入れて同じことをする

21 WinBUGS 実行と文章や口で言っても分かりにくいので, 実際にやってみる実演 WinBUGS the movie! というページもある. 動画で紹介.

22 事後分布平均 %CI[0.72, 2.38] 平均 %CI[0.001, 0.11]

23 まとめベイズ推定によって, うまいこと問題 ( 池 Bのサンプルサイズを考慮, 観測したどのデータも0, など ) を回避していることに注意池 Bの結果を, 池 Aの結果が与えられたという条件のもとで, 池 A= 池 Bとしたとき, 池 Bの結果と同じようになる確率で評価した. 他にどんなモデリングが考えられるか?( 例. 池 Bの方も池 Aと同じようにモデル化して, 二つの平均を比較するとか )

24 例題 2( 階層ベイズ ) 大西洋ビンナガマグロの漁獲量と CPUEのデータにプロダクションモデルをフィットし, そのパラメータの推定結果を用いて,MSYを計算せよ. 最近の漁獲量は適切な水準にあると言えるか?

25 CPUE (Catch Per Unit Effort) CPUE = 漁獲量 / 努力量 CPUE = qp q: 漁獲効率,P: 資源量 ( 個体数 or 個体重量 ) CPUE を使って個体群の増減が分かる年トレンドの推定 CPUE の標準化 (GLM 等 ) CPUE = qp x (Hilborn & Walters 1992)

26 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 になるように魚を獲っていたら, 最大の収穫を得ながら, 個体群を健全に維持できる!

27 状態空間モデル 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 )

28 個体群動態モデルのパラメータ推定ロジスティックモデル 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 を使う

29 状態空間モデルの特徴観測誤差とプロセス誤差を同時に推定することが可能計算機集約的な方法と組み合わせることにより, 柔軟なモデリング推定が可能カルマンフィルタは状態空間モデルの特殊な場合と考えられる一般に計算大変

30 データ : 漁獲量と CPUE 南大西洋ビンナガマグロのデータ

31 WinBUGS 実行実演

32 事後分布 110,000 万回計算して, 最初の 10,000 回を burn-in として取り除いた

33 R と K の相関 K 相関係数 = r

34 自己相関 thinning なし thinning=50 thinning=100 thinning=200

35 2 本の chain でやってみる V R W R = V/W V = pool した chain の変動 W = 個々の chain の平均変動

36 2 本の chain でやってみる

37 R2WinBUGS RからWinBUGSを呼び出して解析を実行し,Rの中でoutputを操作データや初期値の作成が楽ちん即 codaを使える最初は少ない繰り返し数で

38 R2WinBUGS でやってみる

39 coda

40 まとめ個体群動態モデルでは, 従来観測誤差かプロセス誤差のどちらか一方だけを仮定することが多かったが, ベイズ型状態空間モデルを使うことによってどちらの誤差も考慮した分析が比較的容易に行えるいろいろプロットして, うまくいっているか確認することが大事 R2WinBUGSを使うといろいろ便利なことが

41 参考文献 Meyer, R., and Millar, R. B BUGS in Bayesian stock assessments. CJFAS 56: Plummer, M CODA: convergence diagnosis and output analysis for MCMC. R News 6(1): Spiegelhalter et al WinBUGS user manual version 1.4.

kubo2015ngt6 p.2 ( ( (MLE 8 y i L(q q log L(q q 0 ˆq log L(q / q = 0 q ˆq = = = * ˆq = 0.46 ( 8 y 0.46 y y y i kubo (ht

kubo2015ngt6 p.2 ( ( (MLE 8 y i L(q q log L(q q 0 ˆq log L(q / q = 0 q ˆq = = = * ˆq = 0.46 ( 8 y 0.46 y y y i kubo (ht kubo2015ngt6 p.1 2015 (6 MCMC kubo@ees.hokudai.ac.jp, @KuboBook http://goo.gl/m8hsbm 1 ( 2 3 4 5 JAGS : 2015 05 18 16:48 kubo (http://goo.gl/m8hsbm 2015 (6 1 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 2 /