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 / 70 ( 階層ベイズモデル もっと自由な統計モデリングを! 線形モデルの発展 (HBM 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM 推定計算方法 MCMC 最尤推定法 (GLM 最小二乗法 線形モデル 1. ( kubo (http://goo.gl/m8hsbm 2015 (6 3 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 4 / 70 ( ( : ( 0 1 2 3 4 5 6 7 8 1 2 1 3 6 6 1 0 0 y i q i N i y i p(y i q = ( Ni y i q y i (1 q N i y i, q kubo (http://goo.gl/m8hsbm 2015 (6 5 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 6 / 70
kubo2015ngt6 p.2 ( ( (MLE 8 y i L(q q log L(q q 0 ˆq log L(q / q = 0 q -50-45 -40 ˆq = = 73 160 = 0.456 * ˆq = 0.46 ( 8 y 0.46 y 0.54 8 y y i kubo (http://goo.gl/m8hsbm 2015 (6 7 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 8 / 70 ( ( 2. kubo (http://goo.gl/m8hsbm 2015 (6 9 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 10 / 70 : MCMC ( Markov chain Monte Carlo (MCMC (Metropolis method :?? -50-45 -40 log L(q * : q -50-45 -40 q ( kubo (http://goo.gl/m8hsbm 2015 (6 11 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 12 / 70
kubo2015ngt6 p.3 q -49-48 -47-46 -45-44 -47.62-46.38-45.24 0.28 0.29 0.30 0.31 0.32 q 1 q 2 q 3 (1, (2 kubo (http://goo.gl/m8hsbm 2015 (6 13 / 70 : 1 q ( q 0.3 2 q ( q q new 3 q new L(q new L(q L(q new L(q (: q q new L(q new < L(q (: r = L(qnew/L(q q qnew 1 r q 4 2. (q = 0.01 q = 0.99 kubo (http://goo.gl/m8hsbm 2015 (6 14 / 70 q -49-48 -47-46 -45-44 -47.62-46.38-45.24-46 -44-42 0.28 0.29 0.30 0.31 0.32 (MCMC 0.30 0.31 0.32 0.33 0.34 0.35 kubo (http://goo.gl/m8hsbm 2015 (6 15 / 70 q ( MCMC -50-45 -40? q kubo (http://goo.gl/m8hsbm 2015 (6 16 / 70 q q 0 20 40 60 80 100? q q 0 200 400 600 800 1000? q MCMC MCMC?? kubo (http://goo.gl/m8hsbm 2015 (6 17 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 18 / 70
kubo2015ngt6 p.4? q 0 20000 40000 60000 80000 MCMC q q? log L(q -50-45 -40 * L(q p(q = L(q L(q q kubo (http://goo.gl/m8hsbm 2015 (6 19 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 20 / 70 q ( MCMC p(q q 95%! kubo (http://goo.gl/m8hsbm 2015 (6 21 / 70 q MCMC q kubo (http://goo.gl/m8hsbm 2015 (6 22 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 23 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 24 / 70
kubo2015ngt6 p.5 最尤推定と MCMC はちがう! 同じような推定を MCMC でやってみる q の事前分布は一様分布 と考えるとつじつまがあう? 事後分布 p(q Y 0.00 (事後分布 3.97 0.3 0.4 0.5 事前分布 p(q (尤度に比例 p(y はデータ Y が得られる確率 (単なる規格化定数 尤度 L(q 1.98 p(q Y = 0 p(y q p(q p(y p(q Y は何かデータ (Y のもとで何かパラメーター (q が 得られる確率 (事後分布 p(q はあるパラメーター q が得られる確率 (事前分布 p(y q パラメーターを決めたときにデータが得られる確率 ベイズの公式 最尤推定と MCMC はちがう! ベイズ統計にむりやりこじつけてみると? 0.10 ベイズモデル: 尤度 事後分布 事前分布 0.05 同じような推定を MCMC でやってみる 0.6 0.3 0.4 0.5 0.6 0.0 0.2 0.4 0.6 0.8 1.0 生存確率 q 尤度 事前分布 (データが得られる確率 事前分布ってのがよくわからない 尤度 事前分布 kubo (http://goo.gl/m8hsbm 統計モデリング入門 新潟大 2015 (6 同じような推定を MCMC でやってみる 25 / 70 最尤推定と MCMC はちがう! kubo (http://goo.gl/m8hsbm 統計モデリング入門 新潟大 2015 (6 MCMC のためのソフトウェア 26 / 70 Gibbs sampling などが簡単にできるような 以上の説明は MCMC によって得られる結果 3. MCMC のためのソフトウェア は ベイズ統計でいうパラメーターの事後分布 と考えると解釈しやすいかも といったことを Gibbs sampling などが簡単にできるような 事後分布から効率よくサンプリングしたい ばくぜんかつなんとなく対応づける ひとつのこころみでありました 厳密な正当化とかそういったものではありません kubo (http://goo.gl/m8hsbm 統計モデリング入門 新潟大 2015 (6 MCMC のためのソフトウェア 27 / 70 Gibbs sampling などが簡単にできるような 統計ソフトウェア R kubo (http://goo.gl/m8hsbm 統計モデリング入門 新潟大 2015 (6 MCMC のためのソフトウェア 28 / 70 Gibbs sampling などが簡単にできるような 簡単な GLMM なら R だけで推定可能 今回の例題の事後分布 (Y = {yi } はデータ http://www.r-project.org/ p(a, {ri }, s Y 100 i=1 p(yi q(a + ri p(a p(ri s p(s 積分で 個体差 ri を消して 周辺尤度を定義する L(a, s Y = 100 i=1 p(yi q(a + ri p(ri sdri これを最大化する a と s を推定すればよい 経験ベイズ法 (empirical Bayesian method kubo (http://goo.gl/m8hsbm 統計モデリング入門 新潟大 2015 (6 29 / 70 kubo (http://goo.gl/m8hsbm 統計モデリング入門 新潟大 2015 (6 30 / 70
kubo2015ngt6 p.6 R R library(glmmml glmmml( library(lme4 lmer( library(nlme nlme( ( GLMM + (! p(q Y L(q p(q q OK 0.0 0.2 0.4 0.6 0.8 1.0 kubo (http://goo.gl/m8hsbm 2015 (6 31 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 32 / 70 :? MCMC? > post.mcmc[,"a"] # [1] -0.7592-0.7689-0.9008-1.0160-0.8439-1.0380-0.8561-0.9837 [9] -0.8043-0.8956-0.9243-0.9861-0.7943-0.8194-0.9006-0.9513 [17] -0.7565-1.1120-1.0430-1.1730-0.6926-0.8742-0.8228-1.0440... (... 95% 0.0 0.4 0.8 1.2-1.0-0.5 0.0 0.5 1.0 1.5 N = 1200 Bandwidth = 0.0838 kubo (http://goo.gl/m8hsbm 2015 (6 33 / 70 1 : : 2 R package : package : 3 BUGS Gibbs sampler : Gibbs sampler? kubo (http://goo.gl/m8hsbm 2015 (6 34 / 70 MCMC Gibbs sampling? MCMC : MCMC : : MCMC ( : β 1 β 2 Gibbs sampling 1 β 2 2 β 2 β 1 MCMC sampling ( 3 β 1 β 2 MCMC sampling ( 4 2. 3. 9 kubo (http://goo.gl/m8hsbm 2015 (6 35 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 36 / 70
kubo2015ngt6 p.7 : Gibbs sampling ( 9 BUGS Gibbs sampler MCMC step 1 4 6 8 10 β 1 3 4 5 6 7 1.6 1.8 2.0 2.2 2.4 β 1 4 6 8 10 β 2-0.02 0.02 0.06 β 2 3 4 5 6 7 BUGS (+ WinBUGS? step 2 4 6 8 10 1.6 1.8 2.0 2.2 2.4 β 1 4 6 8 10-0.02 0.02 0.06 β 2 OpenBUGS? JAGS OS step 3 4 6 8 10 3 4 5 6 7 3 4 5 6 7 1.6 1.8 2.0 2.2 2.4 β 1 4 6 8 10 3 4 5 6 7 3 4 5 6 7-0.02 0.02 0.06 β 2 kubo (http://goo.gl/m8hsbm 2015 (6 37 / 70 Stan : http://hosho.ees.hokudai.ac.jp/~kubo/ce/bayesianmcmc.html BUGS? kubo (http://goo.gl/m8hsbm 2015 (6 38 / 70 BUGS WinBUGS 1.4.3 データ Y[i] 種子数 8 個のうちの生存数 二項分布 dbin(q,8 生存確率 q 無情報事前分布 BUGS BUGS : Spiegelhalter et al. 1995. BUGS: Bayesian Using Gibbs Sampling version 0.50. kubo (http://goo.gl/m8hsbm 2015 (6 39 / 70 Gibbs sampler BUGS 2004-09-13 ( OpenBUGS Windows kubo (http://goo.gl/m8hsbm 2015 (6 40 / 70 OS JAGS3.4.0 R JAGS (1 / 3 R core team Martyn Plummer Just Another Gibbs Sampler C++ R Linux, Windows, Mac OS X R : library(rjags library(rjags library(r2winbugs # to use write.model( model.bugs <- function( { for (i in 1:N.data { Y[i] ~ dbin(q, 8 # } q ~ dunif(0.0, 1.0 # q } file.model <- "model.bug.txt" write.model(model.bugs, file.model # kubo (http://goo.gl/m8hsbm 2015 (6 41 / 70 # kubo (http://goo.gl/m8hsbm 2015 (6 42 / 70
kubo2015ngt6 p.8 R JAGS (2 / 3 load("data.rdata" list.data <- list(y = data, N.data = length(data inits <- list(q = 0.5 n.burnin <- 1000 n.chain <- 3 n.thin <- 1 n.iter <- n.thin * 1000 model <- jags.model( file = file.model, data = list.data, inits = inits, n.chain = n.chain # kubo (http://goo.gl/m8hsbm 2015 (6 43 / 70 R JAGS (3 / 3 # burn-in update(model, n.burnin # burn in # post.mcmc.list post.mcmc.list <- coda.samples( model = model, variable.names = names(inits, n.iter = n.iter, thin = n.thin # kubo (http://goo.gl/m8hsbm 2015 (6 44 / 70 burn in? 0 100 200 300 400 500 MCMC step kubo (http://goo.gl/m8hsbm 2015 (6 45 / 70 MCMC ˆR = 1.019 MCMC ˆR = 2.520 MCMC MCMC step! kubo (http://goo.gl/m8hsbm 2015 (6 46 / 70 ˆR Gibbs sampling plot(post.mcmc.list gelman.diag(post.mcmc.list R-hat Gelman-Rubin ˆR var ˆ + (ψ y = W var ˆ + (ψ y = n 1 n W + 1 n B W : variance B : variance Gelman et al. 2004. Bayesian Data Analysis. Chapman & Hall/CRC 0.35 0.45 0.55 Trace of q 0 2 4 6 8 10 Density of q 1000 1400 1800 terations 0.30 0.40 0.50 0.60 N = 1000 Bandwidth = 0.00838 kubo (http://goo.gl/m8hsbm 2015 (6 47 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 48 / 70
kubo2015ngt6 p.9! 4. 100 800 403 0.50 0 5 10 15 20 25! 0 2 4 6 8 y i? ( 10 kubo (http://goo.gl/m8hsbm 2015 (6 49 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 50 / 70 (overdispersion : 0 5 10 15 20 25 0 2 4 6 8 y i i N i y i p(y i q i = ( Ni y i q y i i (1 q i N i y i, 0.5 : overdispersion q i kubo (http://goo.gl/m8hsbm 2015 (6 51 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 52 / 70 GLM : r i q i = q(z i q(z = 1/{1 + exp( z} 0.0 0.2 0.4 0.6 0.8 1.0 q(z -6-4 -2 0 2 4 6 z i = a + r i a: r i : i ( z 100 101 (a {r 1, r 2,, r 100 } /! ( kubo (http://goo.gl/m8hsbm 2015 (6 53 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 54 / 70
kubo2015ngt6 p.10 {r i } s = 1.0 s = 1.5 s = 3.0 r i 1 p(r i s = ( exp r2 i 2πs 2 2s 2 p(r i s r i r i r i kubo (http://goo.gl/m8hsbm 2015 (6 55 / 70 : r i (A 0 5 10 15-6 -4-2 0 2 4 6 0 2 4 6 8 y i (B p(r i s s = 0.5 50 {r i} s = 3.0 r i 1 q i = 1+exp( r i 0 5 10 15-6 -4-2 0 2 4 6 2.9 9.9 p(y i q i 0 2 4 6 8 y i kubo (http://goo.gl/m8hsbm 2015 (6 56 / 70 r i r i {r i } 100 r i s = 1.0 s = 1.5 r i p(r i s = s = 3.0 1 2πs 2 exp ( r2 i 2s 2 (A (! (B (C s 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 kubo (http://goo.gl/m8hsbm 2015 (6 57 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 58 / 70 r i! s = 1.0 s = 1.5 r i s = 3.0 p(r i s = 1 2πs 2 exp ( r2 i 2s 2 kubo (http://goo.gl/m8hsbm 2015 (6 59 / 70 全データ 個体個体 3 3 のデータのデータ個体 1 のデータ個体個体 3 3 のデータのデータ個体 2 のデータ {r 1, r 2, r 3,..., r 100 } s local parameter a global parameter? kubo (http://goo.gl/m8hsbm 2015 (6 60 / 70
kubo2015ngt6 p.11 {r i } s (B (C s = 1.0 a, s {r i } s s = 1.5 s = 3.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 s s (non-informative prior 0 < s < 10 4 kubo (http://goo.gl/m8hsbm 2015 (6 61 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 62 / 70 a : 0.0 0.1 0.2 0.3 0.4 ( 0; 1 ( 0; 100-10 -5 0 5 10 (logit a a 種子 8 個のうち Y[i] が生存 生存 q[i] 全個体共通の 平均 a r[i] s hyper parameter kubo (http://goo.gl/m8hsbm 2015 (6 63 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 64 / 70 JAGS JAGS 5. JAGS R JAGS BUGS model { } for (i in 1:N.data { Y[i] ~ dbin(q[i], 8 logit(q[i] <- a + r[i] } a ~ dnorm(0, 1.0E-4 for (i in 1:N.data { r[i] ~ dnorm(0, tau } tau <- 1 / (s * s s ~ dunif(0, 1.0E+4 種子 8 個のうち Y[i] が生存 生存 q[i] 全個体共通の 平均 a r[i] s hyper parameter kubo (http://goo.gl/m8hsbm 2015 (6 65 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 66 / 70
kubo2015ngt6 p.12 JAGS JAGS JAGS > source("mcmc.list2bugs.r" # > post.bugs <- mcmc.list2bugs(post.mcmc.list # bugs 80% 10 interval 5 0 for each 5 10 chain 1 R hat 1.5 2+ a * r[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] s [76] 10 5 0 5 10 1 1.5 2+ * array truncated for lack of space 3 chains, each with 4000 iterations (first 2000 discarded 0.02 0.01 a 0 0.01 0.02 10 5 * r 0 5 10 3.5 s 3 2.5 medians and 80% intervals 1 2 3 4 5 6 7 8 910 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 kubo (http://goo.gl/m8hsbm 2015 (6 67 / 70 bugs post.bugs print(post.bugs, digits.summary = 3 95% 3 chains, each with 4000 iterations (first 2000 discarded, n.thin = 2 n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff a 0.020 0.321-0.618-0.190 0.028 0.236 0.651 1.007 380 s 3.015 0.359 2.406 2.757 2.990 3.235 3.749 1.002 1200 r[1] -3.778 1.713-7.619-4.763-3.524-2.568-1.062 1.001 3000 r[2] -1.147 0.885-2.997-1.700-1.118-0.531 0.464 1.001 3000 r[3] 2.014 1.074 0.203 1.282 1.923 2.648 4.410 1.001 3000 r[4] 3.765 1.722 0.998 2.533 3.558 4.840 7.592 1.001 3000 r[5] -2.108 1.111-4.480-2.775-2.047-1.342-0.164 1.001 2300... ( r[99] 2.054 1.103 0.184 1.270 1.996 2.716 4.414 1.001 3000 r[100] -3.828 1.766-7.993-4.829-3.544-2.588-1.082 1.002 1100 kubo (http://goo.gl/m8hsbm 2015 (6 68 / 70 JAGS JAGS R Trace of a Density of a post.mcmc <- to.mcmc(post.bugs 1.0 0.5 2.0 4.0 0 200 600 1000 terations Trace of s 0 200 600 1000 0.0 0.8 0.0 0.8 1.0 0.0 0.5 1.0 N = 1000 Bandwidth = 0.06795 Density of s 2.0 3.0 4.0 5.0 matrix 0 5 10 15 20 25 0 2 4 6 8 terations N = 1000 Bandwidth = 0.07627 kubo (http://goo.gl/m8hsbm 2015 (6 69 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 70 / 70