kubostat1g p.1 1 (g Hierarchical Bayesian Model kubo@ees.hokudai.ac.jp http://goo.gl/7ci The development of linear models Hierarchical Bayesian Model Be more flexible Generalized Linear Mixed Model (GLMM Incoporating random effects such as individuality parameter MCMC MLE Generalized Linear Model (GLM Always normal distribution? That's non-sense! MSE Linear model : 1 7 8 : kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7 MCMC 1 MCMC Gibbs sampling GLMM GLMM JAGS kubostat1g (http://goo.gl/7ci 1 (g / 7 1. MCMC and logit link function kubostat1g (http://goo.gl/7ci 1 (g / 7 : MCMC seed survivorship, again : 8 : N i = 8 i y i = : 1 7 MCMC 1 7 8 1 1 1 y i kubostat1g (http://goo.gl/7ci 1 (g / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7
kubostat1g p. MCMC binomial distribution q MCMC : i N i y i p(y i q = ( Ni y i q y i (1 q N i y i, q {y i } q likelihood q L(q {y i } = i=1 p(y i q 1 7 8 1 1 1 kubostat1g (http://goo.gl/7ci 1 (g 7 / 7 kubostat1g (http://goo.gl/7ci 1 (g 8 / 7 MCMC MCMC maximum likelihood (MLE L(q ˆq log L(q = log i=1 ( Ni + {y i log(q + (N i y i log(1 q} i=1 q kubostat1g (http://goo.gl/7ci 1 (g 9 / 7 y i L(q q log L(q q ˆq log L(q / q = q - - - ˆq = = 7 1 =..... kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 * MCMC 8 y i ˆq =. ( 8 y. y. 8 y. y i kubostat1g (http://goo.gl/7ci 1 (g 11 / 7 kubostat1g (http://goo.gl/7ci 1 (g 1 / 7
kubostat1g p. : MCMC ( Markov chain Monte Carlo (MCMC (Metropolis method :?? MCMC - - - log L(q *.... : q - - -.... q ( kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 q -9-8 -7 - - - -7. -.8 -..8.9..1. q 1 q q (1, ( kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 : 1 q ( q. q ( q qnew q new L(q new L(q L(q new L(q (: q q new L(q new < L(q (: r = L(q new/l(q q qnew 1 r q. (q =.1 q =.99 kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 q -9-8 -7 - - - -7. -.8 -. - - -.8.9..1. (MCMC..1.... kubostat1g (http://goo.gl/7ci 1 (g 17 / 7 q ( MCMC - - -....? q kubostat1g (http://goo.gl/7ci 1 (g 18 / 7
kubostat1g p. q q.... 8 1? q q.... 8 1? q MCMC MCMC?? kubostat1g (http://goo.gl/7ci 1 (g 19 / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7 MCMC? q.... 8 MCMC q q? log L(q - - - *.... L(q p(q = L(q L(q q kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7 MCMC q ( MCMC p(q q 9%! kubostat1g (http://goo.gl/7ci 1 (g / 7 q MCMC q kubostat1g (http://goo.gl/7ci 1 (g / 7
kubostat1g p. kubostat1g (http://goo.gl/7ci 1 (g / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7 : p(y q p(q p(q Y = p(y p(q Y (Y (q ( p(q q ( p(y q ( p(y Y ( ( (? q? p(q Y L(q p(q q.....8 1. kubostat1g (http://goo.gl/7ci 1 (g 7 / 7 kubostat1g (http://goo.gl/7ci 1 (g 8 / 7 Gibbs sampling MCMC. Gibbs sampling kubostat1g (http://goo.gl/7ci 1 (g 9 / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7
kubostat1g p. 統計ソフトウェア R 簡単な GLMM なら R だけで推定可能 今回の例題の事後分布 (Y = {yi } はデータ http://www.r-project.org/ p(a, {ri }, s Y 1 i=1 p(yi q(a + ri p(a p(ri s p(s 積分で 個体差 ri を消して 周辺尤度を定義する L(a, s Y = 1 i=1 p(yi q(a + ri p(ri sdri これを最大化する a と s を推定すればよい 経験ベイズ法 (empirical Bayesian method kubostat1g (http://goo.gl/7ci 1 / 7 kubostat1g (http://goo.gl/7ci / 7 しかし R だけ では限界があるかも そこで MCMC による事後分布からのサンプリング! R にはいろいろな GLMM の最尤推定関数が準備さ library(lme の lmer(.. しかし もうちょっと複雑な GLMM たとえば個体 差 + 地域差をいれた統計モデルの最尤推定はかなり 難しい (ヘンな結果が得られたりする... 事前分布 p(q. 尤度 L(q.97 library(nlme の nlme( (正規分布のみ 事後分布 p(q Y 1.98 library(glmmml の glmmml( れている.1.........8 1. 生存確率 q アルゴリズムにしたがって 乱数を発生させていくだけで OK 積分がたくさん入っている尤度関数の評価がしんどい kubostat1g (http://goo.gl/7ci / 7 kubostat1g (http://goo.gl/7ci / 7 どのようなソフトウェアで MCMC 計算するか? 再確認: 事後分布からのサンプル って何の役にたつの? > post.mcmc[,"a"] # 事後分布からのサンプルを表示 [1] -.79 -.789 -.98-1.1 -.89-1.8 -.81 -.987 [9] -.8 -.89 -.9 -.981 -.79 -.819 -.9 -.91 [17] -.7-1.11-1. -1.17 -.9 -.87 -.88-1.... (以下略... 1 自作プログラム 利点: 問題にあわせて自由に設計できる 欠点: 階層ベイズモデル用の MCMC プログラミング けっこうめん どう これらのサンプルの平均値 中央値 9% 区間を 調べることで事後分布の概要がわかる R のベイズな package 利点: 空間ベイズ統計など便利な専用 package がある 1. 欠点: 汎用性 とぼしい.8 BUGS で Gibbs sampler なソフトウェア. 利点: 幅ひろい問題に適用できて 便利. 欠点というほどでもないけど 多少の勉強が必要 -1. -... 1. N = 1 Bandwidth =.88 kubostat1g (http://goo.gl/7ci 1. えーっと Gibbs sampler って何? / 7 kubostat1g (http://goo.gl/7ci / 7
kubostat1g p.7 さまざまな MCMC アルゴリズム Gibbs sampling とは何か? MCMC アルゴリズムのひとつ いろいろな MCMC メトロポリス法: 試行錯誤で値を変化させていく 複数のパラメーターの MCMC サンプリングに使う 例: パラメーター β1 と β の Gibbs sampling MCMC 1 メトロポリス ヘイスティングス法: その改良版 ギブス サンプリング: 条件つき確率分布を使った MCMC β に何か適当な値を与える β の値はそのままにして その条件のもとでの β1 の MCMC sampling をする (条件つき事後分布 β1 の値はそのままにして その条件のもとでの β の MCMC sampling をする (条件つき事後分布 複数の変数 (パラメーター 状態 を効率よくサンプリング.. をくりかえす 教科書の第 9 章の例題で説明 kubostat1g (http://goo.gl/7ci 図解: Gibbs sampling 7 / 7 記述できるソフトウェア -. 7. β. 7 8 1 1. 1.8... β1 -. 8 1 step 7. β. 7 WinBUGS ありがとう さようなら? OpenBUGS 予算が足りなくて停滞? JAGS お手軽で良い どんな OS でも動く Stan たぶん 次 はこれ 今日は紹介しませんが 8 1 リンク集: 1. 1.8... β1 kubostat1g (http://goo.gl/7ci 7 -. 8 1 step BUGS 言語 (+ っぽいもの でベイズモデルを 8 1 8 / 7 β のサンプリング 8 1 便利な BUGS 汎用 Gibbs sampler たち (統計モデリング入門の第 9 章 1. 1.8... β1 step 1 kubostat1g (http://goo.gl/7ci β1 のサンプリング MCMC. β. 7 http://hosho.ees.hokudai.ac.jp/~kubo/ce/bayesianmcmc.html えーと BUGS 言語って何? 9 / 7 kubostat1g (http://goo.gl/7ci / 7 なんとなく使われ続けている WinBUGS 1.. このベイズモデルを BUGS 言語で記述したい データ Y[i] 種子数８個のうちの生存数 二項分布 dbin(q,8 おそらく世界でもっともよく使われている Gibbs sampler BUGS 言語コード BUGS 言語の実装 -9-1 に最新版 (ここで開発停止 OpenBUGS ソースなど非公開 無料 ユーザー登録不要 生存確率 q Windows バイナリーとして配布されている 無情報事前分布 歴史を変えたソフトウェアだけど 開発も停止していることだし まあ もう ごくろうさま ということで 矢印は手順ではなく 依存関係をあらわしている BUGS 言語: ベイズモデルを記述する言語 Spiegelhalter et al. 199. BUGS: Bayesian Using Gibbs Sampling version.. kubostat1g (http://goo.gl/7ci 1 / 7 kubostat1g (http://goo.gl/7ci / 7
kubostat1g p.8 Gibbs sampling Gibbs sampling OS JAGS.. R JAGS (1 / R core team Martyn Plummer Just Another Gibbs Sampler C++ R Linux, Windows, Mac OS X R : library(rjags library(rjags library(rwinbugs # to use write.model( model.bugs <- function( { for (i in 1:N.data { Y[i] ~ dbin(q, 8 # } q ~ dunif(., 1. # q } file.model <- "model.bug.txt" write.model(model.bugs, file.model # kubostat1g (http://goo.gl/7ci 1 (g / 7 # kubostat1g (http://goo.gl/7ci 1 (g / 7 Gibbs sampling Gibbs sampling R JAGS ( / R JAGS ( / load("data.rdata" list.data <- list(y = data, N.data = length(data inits <- list(q =. n.burnin <- 1 n.chain <- n.thin <- 1 n.iter <- n.thin * 1 model <- jags.model( file = file.model, data = list.data, inits = inits, n.chain = n.chain # kubostat1g (http://goo.gl/7ci 1 (g / 7 # 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 # kubostat1g (http://goo.gl/7ci 1 (g / 7 Gibbs sampling Gibbs sampling burn in?.... 1 MCMC step kubostat1g (http://goo.gl/7ci 1 (g 7 / 7 MCMC ˆR = 1.19 MCMC ˆR =. MCMC MCMC step! kubostat1g (http://goo.gl/7ci 1 (g 8 / 7
kubostat1g p.9 Gibbs sampling Gibbs sampling ˆ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.. Bayesian Data Analysis. Chapman & Hall/CRC... Trace of q 8 1 Density of q 1 1 18 Iterations.... N = 1 Bandwidth =.88 kubostat1g (http://goo.gl/7ci 1 (g 9 / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7 GLMM GLMM GLMM GLMM!. GLMM 1 8. GLMM hierarchical Bayesian 1 1! 8 y i? ( 1 kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7 GLMM GLMM GLMM GLMM (overdispersion : 1 1 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,. : overdispersion q i kubostat1g (http://goo.gl/7ci 1 (g / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7
kubostat1g p.1 GLMM と階層ベイズモデル GLMM のベイズモデル化 GLMM と階層ベイズモデル 個々の個体差 ri を最尤推定するのはまずい GLM わざ: ロジスティック関数で表現する生存確率 生存確率 qi = q(zi をロジスティック関数 q(z = 1/{1 + exp( z} で表現..8 1.... q(z - - - 1 個体の生存確率を推定するためにパラメーター 11 個 (a と {r1, r,, r1 } を推定すると 個体ごとに生存数 / 種子数を計算していることと同 じ! ( データのよみあげ と同じ z 線形予測子 zi = a + ri とする そこで 次のように考えてみる パラメーター a: 全体の平均 パラメーター ri : 個体 i の個体差 (ずれ kubostat1g (http://goo.gl/7ci GLMM と階層ベイズモデル GLMM のベイズモデル化 / 7 kubostat1g (http://goo.gl/7ci GLMM のベイズモデル化 GLMM と階層ベイズモデル suppose {ri } follow the Gausssian distribution {ri } のばらつきは正規分布だと考えてみる / 7 GLMM のベイズモデル化 ひとつの例示: 個体差 ri の分布と過分散の関係 (A 個体差のばらつきが小さい場合 (B 個体差のばらつきが大きい場合 s = 1. p(ri s が生成した 個体ぶんの {ri } s =. s = 1. - - IIII IIIIII IIIIIIII - ri s =. 確率 qi = - - - s =. I-I I -I I III-IIIIIIIIIIIII II I IIII ri 1 1+exp( ri の二項乱数を発生させる すればよいでしょう ri がゼロにちかい個体はわりと ありがち で ri GLMM と階層ベイズモデル 1 1 8 生存種子数 yi の絶対値が大きな個体は相対的に あまりいない kubostat1g (http://goo.gl/7ci 標本分散 9.9 1 1 p(yi qi が 生成した生 存種子数の 一例 この確率密度 p(ri s は ri の 出現しやすさ をあらわしていると解釈 標本分散.9 ( 1 r p(ri s = exp i s πs 観察された個体数 個体差 ri 7 / 7 kubostat1g (http://goo.gl/7ci GLMM のベイズモデル化 GLMM と階層ベイズモデル これは ri の事前分布の指定 ということ 8 生存種子数 yi 8 / 7 GLMM のベイズモデル化 ベイズ統計モデルでよく使われる三種類の事前分布 たいていのベイズ統計モデルでは ひとつのモデルの中で複数の種類の 前回の授業で {ri } は正規分布にしたがうと仮定したが ベイズ統計モデリングでは 1 個の ri たちに 共通する事前分布として正規分布 を指定した ということになる 事前分布を混ぜて使用する (A 主観的な事前分布 s = 1. (B 無情報事前分布 (C 階層事前分布 (できれば使いたくない! s = 1. s =. - - - 信じる s によって 変わる わからない 個体差 ri ( 1 r p(ri s = exp i s πs kubostat1g (http://goo.gl/7ci. 9 / 7....8 1. kubostat1g (http://goo.gl/7ci.....8 1......8 1. / 7
kubostat1g p.11 GLMM GLMM GLMM GLMM r i! s = 1. s = 1. r i s =. p(r i s = 1 πs exp ( r i s kubostat1g (http://goo.gl/7ci 1 (g 1 / 7 全データ 個体個体 のデータのデータ個体 1 のデータ個体個体 のデータのデータ個体 のデータ {r 1, r, r,..., r 1 } s local parameter a global parameter? kubostat1g (http://goo.gl/7ci 1 (g / 7 GLMM GLMM GLMM GLMM {r i } s (B (C a, s {r i } s s = 1. s = 1. s =......8 1......8 1. global local kubostat1g (http://goo.gl/7ci 1 (g / 7 s s (non-informative prior < s < 1 kubostat1g (http://goo.gl/7ci 1 (g / 7 GLMM GLMM GLMM GLMM a :..1... ( ; 1 ( ; 1-1 - 1 (logit a a 種子 8 個のうち Y[i] が生存 生存 q[i] 全個体共通の 平均 a r[i] s hyper parameter kubostat1g (http://goo.gl/7ci 1 (g / 7 kubostat1g (http://goo.gl/7ci 1 (g / 7
kubostat1g p.1 JAGS JAGS. 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(, 1.E- for (i in 1:N.data { r[i] ~ dnorm(, tau } tau <- 1 / (s * s s ~ dunif(, 1.E+ } 種子 8 個のうち Y[i] が生存 生存 q[i] 全個体共通の 平均 a r[i] s hyper parameter kubostat1g (http://goo.gl/7ci 1 (g 7 / 7 kubostat1g (http://goo.gl/7ci 1 (g 8 / 7 JAGS JAGS JAGS > source("mcmc.listbugs.r" # > post.bugs <- mcmc.listbugs(post.mcmc.list # bugs 8% 1 interval for each 1 chain 1 R hat 1. + a * r[1] [] [] [] [] [] [7] [8] [9] [1] [11] [1] [1] [1] [1] [1] [17] [18] [19] [] [1] [] [] [] [] [] [7] [8] [9] [] [1] [] [] [] [] [] [7] [8] [9] [] [1] [] [] [] [] [] [7] [8] [9] [] [1] [] [] [] [] [] [7] [8] [9] [] [1] [] [] [] [] [] [7] [8] [9] [7] [71] [7] [7] [7] [7] s [7] 1 1 1 1. + * array truncated for lack of space chains, each with iterations (first discarded..1 a.1. 1 * r 1. s. medians and 8% intervals 1 7 8 91 1 1 1 18 8 8 kubostat1g (http://goo.gl/7ci 1 (g 9 / 7 bugs post.bugs print(post.bugs, digits.summary = 9% chains, each with iterations (first discarded, n.thin = n.sims = iterations saved mean sd.% % % 7% 97.% Rhat n.eff a..1 -.18 -.19.8..1 1.7 8 s.1.9..77.99..79 1. 1 r[1] -.778 1.71-7.19 -.7 -. -.8-1. 1.1 r[] -1.17.88 -.997-1.7-1.118 -.1. 1.1 r[].1 1.7. 1.8 1.9.8.1 1.1 r[].7 1.7.998..8.8 7.9 1.1 r[] -.18 1.111 -.8 -.77 -.7-1. -.1 1.1... ( r[99]. 1.1.18 1.7 1.99.71.1 1.1 r[1] -.88 1.7-7.99 -.89 -. -.88-1.8 1. 11 kubostat1g (http://goo.gl/7ci 1 (g 7 / 7 JAGS JAGS R Trace of a Density of a post.mcmc <- to.mcmc(post.bugs 1.... 1 Iterations Trace of s 1..8..8 1... 1. N = 1 Bandwidth =.79 Density of s.... matrix 1 1 8 Iterations N = 1 Bandwidth =.77 kubostat1g (http://goo.gl/7ci 1 (g 71 / 7 kubostat1g (http://goo.gl/7ci 1 (g 7 / 7
kubostat1g p.1. Yay! Be more flexible The development of linear models Hierarchical Bayesian Model Generalized Linear Mixed Model (GLMM Incoporating random effects such as individuality parameter MCMC MLE Generalized Linear Model (GLM Always normal distribution? That's non-sense! MSE Linear model kubostat1g (http://goo.gl/7ci 1 (g 7 / 7 kubostat1g (http://goo.gl/7ci 1 (g 7 / 7