2006-12-09 1/22 R MCMC R 1. 2. R MCMC? 3. Gibbs sampler : kubo@ees.hokudai.ac.jp http://hosho.ees.hokudai.ac.jp/ kubo/
2006-12-09 2/22 : ( ) : : ( ) : (?) community ( )
2006-12-09 3/22 :? 1. ( ) 2. ( ) 3. or?? MCMC? 1.? 2. MCMC? 3.? ( ) R?!
2006-12-09 4/22 ( ) [ ] Bayes Bayes [ ] fixed effects Bayes (GLMM) [ (GLM)] +, fixed effects [ ] +
2006-12-09 5/22? random effects ( ) nest random effects Plot A Plot B : MCMC : Gaussian Random Field
2006-12-09 6/22 MCMC?! Markov Chain Monte Carlo? Gibbs (sampling) : (sampling!) sample random sample set log likelihood (non constant par 11200 10600 leaf area index 4.2 4.4 4.6 step (burn-in) 0 1000 2000 3000 4000 (K MCMC step) sampling point 0 1000 2000 3000 4000 (K MCMC step)
2006-12-09 7/22 ( ) ( ) ( ) ( ) p(β, α y) p(y β)p(β α)p(α) : Markov Chain Monte Carlo (MCMC) MCMC 1: Metropolis-Hastings MCMC 2: Gibbs sampler ( ) Likelihood(α y) p(y β)p(β α)dβ : α ( ) : (GLMM)? ( :. 2004. )
2006-12-09 8/22 R MCMC / Gibbs sampler MCMC? ( ) R package: library(mcmcpack) ( ) Gibbs sampler (R ) WinBUGS OpenBUGS JAGS WinBUGS OpenBUGS WinBUGS, 2004 OpenBUGS WinBUGS project, GPL
2006-12-09 9/22 : WinBUGS 1.4.1 Gibbs sampler BUGS adaptive rejection sampler 2004-09-13 ( OpenBUGS ) Windows Linux WINE MacOS X Darwine GUI (Linux ) R R2WinBUGS ( )
2006-12-09 10/22 BUGS Spiegelhalter et al. 1995. BUGS: Bayesian Using Gibbs Sampling version 0.50. model { mu dnorm(0, 1.0E-2) tau dgamma(1.0e-3, 1.0E-3) for (i in 1:n.samples) { re[i] dnorm(0.0, tau) p[i] <- 1.0 / (1.0 + exp(-(mu + re[i]))) n.seeds[i] dbin(p[i], n.ovules[i]) } } JAGS ;
2006-12-09 11/22 BUGS (?) : zero-inflated Poisson (ZIP) model BUGS :?
2006-12-09 12/22 GPL WinBUGS : OpenBUGS 2.2.0 Thomas Andrew WinBUGS OpenBUGS is still in development and suffers frequent crashes. Component Pascal BlackBox Component Builder Windows Linux R library(brugs) Windows R!!
2006-12-09 13/22 R (?) Gibbs sampler: JAGS 0.97 R core team Martyn Plummer Just Another Gibbs Sampler C++ R Vine Linux RPM package Windows JAGS MacOS X R : library(rjags)
2006-12-09 14/22 JAGS JAGS directed cycle ( ) model { x dnorm(y, tau) y dnorm(x, tau) }! WinBUGS? R
2006-12-09 15/22 : JAGS WinBUGS
2006-12-09 16/22 : [ ] i {1,, 100} n.ovules[i] = 8 n.seeds[i] p[i]: x = 5 : : logit(p[i]) = 0 + re[i] random effects: re[i] N(0, 1/0.5) = 0.5, = 4
2006-12-09 17/22 : 0 1 2 3 4 5 6 7 8 ( ) 9 9 7 9 14 18 15 15 4 100 0.2 2.2 8.6 19.2 27.0 24.2 13.6 4.4 0.6 100 [ ] = 0.529 p = 0.529 overdispersion ( )! random effects x = 5 ( )
2006-12-09 18/22 JAGS 1. BUGS model 2. R 3. JAGS (foo.cmd) 4. jags foo.cmd 5. JAGS R library(coda) read.coda() (mcmc ) 6. mcmc plot(), summary(),
2006-12-09 19/22 library(coda): MCMC Convergence Diagnosis and Output Analysis for MCMC R package ( S-plus ) Martyn Plummer MCMC mcmc, mcmc.list MCMC
2006-12-09 20/22 mcmc summary > summary(r.mcmc) 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE mu 0.254 0.187 0.00835 0.0254 tau 0.355 0.090 0.00402 0.0101 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% mu -0.108 0.125 0.251 0.372 0.622 tau 0.208 0.290 0.349 0.409 0.560
2006-12-09 21/22 WinBUGS : R2WinBUGS 1. BUGS model 2. R2WinBUGS package R 3. R 2. 4. bugs 5. plot() summary() 6. mcmc / mcmc.list
2006-12-09 22/22 : R MCMC? ( MCMC ) WinBUGS R2WinBUGS? WinBUGS R! ( ) JAGS R ( Plummer?) MCMC geor : Gaussian Random Field GRF MCMC : Langevin-Hastings? : OpenBUGS?... or... Umacs + rv B?