MCMC MCMC R[26] MCMCpack[22] JAGS[25] MCMCpack Stan[29] * MCMC? MCMC x θ π(θ x) π(θ x) = f(x θ)π(θ) f(x θ)π(θ)dθ (1) f(x θ) π(θ) θ

Size: px

Start display at page:

Download "MCMC MCMC R[26] MCMCpack[22] JAGS[25] MCMCpack Stan[29] * MCMC? MCMC x θ π(θ x) π(θ x) = f(x θ)π(θ) f(x θ)π(θ)dθ (1) f(x θ) π(θ) θ"

付从
5 years ago
Views:

1 MCMC MCMC R[26] MCMCpack[22] JAGS[25] MCMCpack Stan[29] * MCMC? MCMC x θ π(θ x) π(θ x) = f(x θ)π(θ) f(x θ)π(θ)dθ (1) f(x θ) π(θ) θ π(θ x) f(x θ)π(θ) (2) MCMC hiroki@affrc.go.jp *1 R 3.2.2, MCMCpack 1.3-3, JAGS 3.4.0, Stan 2.8.0

2 MCMC=MC+MC MCMC? MCMC MC (Markov Chain) MC (Monte Carlo) Markov chain x t+1 x t : x t 1 Monte Carlo MCMC

3 1.2 MCMC MCMC Metropolis-Hastings Gibbs sampler [1, 10, 34, 37] Hamiltonian Monte Carlo Hybrid Monte Carlo [1, 7, 35, 37] 2 MCMC R : MCMCpack : WinBUGS, OpenBUGS, JAGS, Stan Stan BUGS (Bayesian inference Using Gibbs Sampling) [21, 32] Stan BUGS R Python PyMC C MCMCpack : : R CRAN Metropolis sampler MCMClogit() MCMCpoisson() MCMC MCMChregress() MCMChlogit() MCMChpoisson() 30 MCMCmetrop1R() [22] 2.2 WinBUGS :

4 : Windows Wine *2 OS X Linux BSD GUI R2WinBUGS R WinBUGS key key Windows WinBUGS14.exe Windows Vista C:\Program Files C:\Program Files *3 64 Windows OS X Mac Wine WinBUGS OS X WinBUGS *4 OS X Wine.app *5 MacPors *6 Wine Windows (WinBUGS14.exe) Wine WinBUGS ~/.wine/drive_c/program Files Linux Wine Ubuntu Wine (WinBUGS14.exe) WinBUGS *2 Wine Is Not Emulator Wine 1.2 WinBUGS Wine Windows Windows ( ) Windows *3 * *5 *6

5 BSD UNIX Linux 2.3 OpenBUGS : GPL Black Box Component Builder *7 Windows Linux OS X Wine BRUGS R2OpenBUGS R [33] CRAN Windows Windows 64 Windows Windows XP 8 OS X Linux WinBUGS Wine Windows Linux 2.4 JAGS : : : GPL C C++ *8 BLAS LAPACK rjags R2jags runjags R CRAN WinBUGS *9 WinBUGS * 10 *7 * *9 *10

6 OpenBUGS * 11 Windows OS X Linux 2.5 Stan : : BSD GPL3 RStan R CRAN Pyhton PyStan CmdStan Stan C++ OS X, Windows, Linux Hamiltonian Monte Carlo [1, 7, 35] MCMC BUGS RStan CRAN CmdStan C++ *11

7 3 MCMC 3.1 : MCMC x = (3, 1, 4, 3, 3, 6, 4, 1, 6, 4, 1, 7, 4, 4, 1, 4, 0, 3, 9, 4) λ MCMCpoisson() R MCMCpack MCMCpoisson() R example1.r MCMCpack > + 1 > library ( MCMCpack ) post1 3 1 > post1 <- vector (" list ", 3) MCMCpoisson() MCMCpoisson() log(λ) MCMC 1 > post1 [[1]] <- MCMCpoisson (x ~ 1, beta. start = 1, 2 + burnin = 0, mcmc = 400, 3 + thin = 1, 4 + tune = 0.8, seed = 1117, 5 + verbose = 20) x 1 glm() 1 X Poisson(λ) log λ = C

8 C * 12 improper beta.start log λ burnin burn-in 0 mcmc thin 1 * 13 tune seed verbose 20 Markov chain 2 1 > plot ( post1 [[1]], density = FALSE, col = 1, las = 1) Trace of (Intercept) Iterations 2 MCMCpoissonR() Markov chain 50 burn-in tune [27] *12 b0 B0 : *13 thinning [19]

9 1 > post1 [[2]] <- MCMCpoisson (x ~ 1, beta. start = 5, 2 + burnin = 0, mcmc = 400, thin = 1, 3 + tune = 0.8, seed = 1123, 4 + verbose = 20) 5 > post1 [[3]] <- MCMCpoisson (x ~ 1, beta. start = 10, 6 + burnin = 0, mcmc = 400, thin = 1, 7 + tune = 0.8, seed = 1129, 8 + verbose = 20) mcmc.list 3 1 > post1. mcmc <- mcmc. list ( post1 ) 2 > plot ( post1. mcmc, density = FALSE, ylim = c(0, 10), las = 1) Trace of (Intercept) Iterations 3 Markov chain 3 burn-in 500 burn-in Markov chain = > post2 <- vector (" list ", 3) 2 > burnin < > mcmc < > thin <- 1 5 > tune <- 2 6 > verbose = 0 7 > post2 [[1]] <- MCMCpoisson (x ~ 1, beta. start = 1,

10 8 + burnin = burnin, mcmc = mcmc, 9 + thin = thin, 10 + tune = tune, seed = 1117, 11 + verbose = verbose ) 12 > post2 [[2]] <- MCMCpoisson (x ~ 1, beta. start = 5, 13 + burnin = burnin, mcmc = mcmc, 14 + thin = thin, 15 + tune = tune, seed = 1123, 16 + verbose = verbose ) 17 > post2 [[3]] <- MCMCpoisson (x ~ 1, beta. start = 10, 18 + burnin = burnin, mcmc = mcmc, 19 + thin = thin, 20 + tune = tune, seed = 1129, 21 + verbose = verbose ) 22 > post2. mcmc <- mcmc. list ( post2 ) 23 > 24 > plot ( post2. mcmc ) 4 1 > plot ( post2. mcmc ) Trace of (Intercept) Density of (Intercept) Iterations N = 1000 Bandwidth = Markov chain λ 1 > summary ( post2. mcmc ) 2 3 Iterations = 501: Thinning interval = 1

11 5 Number of chains = 3 6 Sample size per chain = Empirical mean and standard deviation for each variable, 9 plus standard error of the mean : Mean SD Naive SE Time - series SE Quantiles for each variable : % 25% 50% 75% 97. 5% % * 14 log λ % λ exp(1.33) = MCMCmetrop1R() MCMCmetrop1R() MCMCpoisson() MCMCmetrop1R() MCMCmetrop1R() 1 > LogPoisFun <- function ( lambda, x) { 2 + # Poisson distribution : p(x) = lambda ^x exp (- lambda )/x! 3 + if ( lambda >= 0) { # lambda must be non - negative 4 + # prior : dunif (0, 10^4) 5 + log ( ifelse ( lambda >= 0 & lambda < 10^4, 10^ -4, 0)) # log likelihood 7 + sum ( log ( lambda ^x * exp (- lambda ) / factorial (x ))) 8 + } else { 9 + -Inf 10 + } 11 + } λ λ Uniform(0, 10 4 ) *14

12 Pr(λ) = { λ < (, ) improper prior f(x) = λ x e λ /x! λ X L(X λ) 3 N λ X i e λ L(X λ) = X i! i=1 (3) Pr(λ X) L(X λ)pr(λ) (4) { N i=1 λx i e λ L(X λ)pr(λ) = X i! λ < (5) R LogPoisFun() 5 MCMCmetrop1R() fun help(mcmcmetrop1r) LogPoisFun MCMCpoisoon() 3 Markov chain lapply() 1 > chains <- 1:3 2 > inits <- c(1, 10, 20) 3 > seeds <- c (1117, 1123, 1129) 4 > post3 <- lapply ( chains, 5 + function ( chain ) { 6 + MCMCmetrop1R ( fun = LogPoisFun, 7 + theta. init = inits [ chain ], 8 + burnin = 1000, mcmc = 1000, 9 + thin = 1, 10 + tune = 2, seed = seeds [ chain ], 11 + verbose = 0, logfun = TRUE, 12 + x = x) 13 + }) 14 15

13 16 17 The Metropolis acceptance rate was The Metropolis acceptance rate was The Metropolis acceptance rate was > summary ( post3. mcmc ) 2 3 Iterations = 1001: Thinning interval = 1 5 Number of chains = 3 6 Sample size per chain = Empirical mean and standard deviation for each variable, 9 plus standard error of the mean : Mean SD Naive SE Time - series SE Quantiles for each variable : % 25% 50% 75% 97. 5% Gelman-Rubin ( ˆR; Rhat) R coda gelman.diag() 1 chain help(gelman.diag) ˆR 1.1

14 Iterations N = 1000 Bandwidth = > gelman. diag ( post3. mcmc ) 2 Potential scale reduction factors : 3 4 Point est. Upper C.I. 5 [1,] Geweke coda geweke.diag() Z Z > % help(geweke.diag) 1 > geweke. diag ( post3. mcmc ) 2 [[1]] 3 4 Fraction in 1 st window = Fraction in 2 nd window = var [[2]] Fraction in 1 st window = Fraction in 2 nd window = var1

15 [[3]] Fraction in 1 st window = Fraction in 2 nd window = var coda boa (Bayesian Output Analysis) 3.2 MCMC N(=11) (x = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) ) K(=10) x 10 y = (1, 2, 2, 6, 4, 5, 8, 9, 9, 9, 10) x y x y y x 6 3.2

16 3.2.1 MCMCpack MCMCpack R example2.r MCMCpackage 1 > library ( MCMCpack ) data y =0, =1 1 > x <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) 2 > y <- c(1, 2, 2, 6, 4, 5, 8, 9, 9, 9, 10) 3 > k < > data <- data. frame (x = rep (x, each = k), 5 + y = c( sapply (y, function (i) 6 + c( rep (0, k - i), rep (1, i ))))) y x MCMClogit() MCMClogit() MCMCpack chain 3 burn-in 2000 burn-in > chains <- 1:3 2 > inits <- c(1, 10, 20) 3 > seeds <- c(12, 123, 1234) 4 > post <- lapply ( chains, 5 + function ( chain ) { 6 + MCMClogit ( y ~ x, 7 + data = data, 8 + burnin = 2000, mcmc = 2000, 9 + thin = 2, 10 + tune = 1.1, verbose = 500, 11 + seed = seeds [ chain ]) 12 + }) MCMClogit iteration 1 of beta = Metropolis acceptance rate for beta = : 22 :

17 23 24 MCMClogit iteration 3501 of beta = Metropolis acceptance rate for beta = The Metropolis acceptance rate for beta was mcmc.list ( 7) 1 > post. mcmc <- mcmc. list ( post ) 2 > plot ( post. mcmc ) Trace of (Intercept) Density of (Intercept) Iterations N = 1000 Bandwidth = Trace of x Density of x Iterations N = 1000 Bandwidth = coda plot (MCMCpack)

18 Gelman-Rubin (ˆR) 1 > gelman. diag ( post. mcmc ) 2 Potential scale reduction factors : 3 4 Point est. Upper C.I. 5 ( Intercept ) x Multivariate psrf > summary ( post. mcmc ) 2 3 Iterations = 2001: Thinning interval = 2 5 Number of chains = 3 6 Sample size per chain = Empirical mean and standard deviation for each variable, 9 plus standard error of the mean : Mean SD Naive SE Time - series SE 12 ( Intercept ) x Quantiles for each variable : % 25% 50% 75% 97. 5% 18 ( Intercept ) x % x % JAGS JAGS R example2_jags.r R JAGS rjags 1 > library ( rjags )

19 R 1 > k < > x <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) 3 > y <- c(1, 2, 2, 6, 4, 5, 8, 9, 9, 9, 10) 4 > n <- length (x) BUGS example2_model.txt MCMCpack MCMClogit() 1 model { 2 for (i in 1:N) { 3 Y[i] ~ dbin (p[i], K) 4 logit (p[i]) <- beta + beta.x * X[i] 5 } 6 beta ~ dnorm (0, 1.0 E -6) 7 beta.x ~ dnorm (0, 1.0E -6) 8 } 2 5 for for X[i], Y[i], p[i] i i 1, 2,..., N dbin(p, n) Binomial(n, p) p n ~ ( ) logit(p) logit(p) = log(p/(1 p)) Y K p p X 6 7 beta beta.x beta beta.x Normal(0, 10 6 ) BUGS dnorm() 1 2 ( ~ ) ( <- ) Pr(β, β x X, Y ) L(X, Y β, β x ) Pr(β) Pr(β x ) (6) N L(X, Y β, β x ) = Binomial(Y i K, p(x i )) = i=1 N i=1 K! Y i!(k Y i )! p(x i) Y i (1 p(x i )) K Y i (7)

20 p(x i ) = logit 1 (β + β x X i )) = exp(β + β xx i ) 1 + exp(β + β x X i ) Pr(β) = Normal(0, 10 6 ) (9) Pr(β x ) = Normal(0, 10 6 ) (10) (8) chain 1 > ## Number of chains 2 > n. chains <- 3 3 > 4 > ## Initial values 5 > inits <- vector (" list ", n. chains ) 6 > inits [[1]] <- list ( beta = -10, beta.x = 0, 7 +. RNG. seed = 314, 8 +. RNG. name = " base :: Mersenne - Twister ") 9 > inits [[2]] <- list ( beta = -5, beta.x = 2, RNG. seed = 3141, RNG. name = " base :: Mersenne - Twister ") 12 > inits [[3]] <- list ( beta = 0, beta.x = 4, RNG. seed = 31415, RNG. name = " base :: Mersenne - Twister ") 15 > 16 > ## Model file 17 > model. file <- " example2_model. txt " 18 > 19 > ## Parameters 20 > pars <- c(" beta ", " beta.x") n.chains inits rjags chain 3 beta beta.x chain (.RNG.seed) (.RNG.name) (model.file) (pars) beta beta.x jags.mode() JAGS chain MCMC (adaptation) adapt

21 R 100% 1 > model <- jags. model ( file = model. file, 2 + data = list (N = n, K = k, 3 + X = x, Y = y), 4 + inits = inits, n. chains = n. chains, 5 + n. adapt = 1000) % update() burn-in n.adapt n.iter MCMCpack burn-in 1 > update ( model, n. iter = 1000) 2 ************************************************** 100% coda.samples() MCMC coda.samples() coda jags jags.samples() chain > post. samp <- coda. samples ( model, n. iter = 3000, thin = 3, 2 + variable. names = pars ) 3 ************************************************** 100% 8 1 > plot ( post. samp ) chain chain gelman.diag() Gelman-Rubin > gelman. diag ( post. samp ) 2 Potential scale reduction factors : 3 4 Point est. Upper C.I. 5 beta beta. x Multivariate psrf

22 Trace of beta Density of beta Iterations N = 1000 Bandwidth = Trace of beta.x Density of beta.x Iterations N = 1000 Bandwidth = coda plot (JAGS) 1 > summary ( post. samp ) 2 3 Iterations = 2003: Thinning interval = 3 5 Number of chains = 3 6 Sample size per chain = Empirical mean and standard deviation for each variable, 9 plus standard error of the mean : Mean SD Naive SE Time - series SE 12 beta beta. x Quantiles for each variable : 16

23 17 2.5% 25% 50% 75% 97. 5% 18 beta beta. x beta % beta.x % X Y Y 1 > beta <- unlist ( post. samp [, " beta "]) 2 > beta.x <- unlist ( post. samp [, " beta.x "]) 3 > 4 > new. x <- seq (0, 10, len = 100) 5 > logit.p <- beta + beta.x %o% new.x 6 > exp.y <- k * exp ( logit.p) / (1 + exp ( logit.p)) 7 > y. mean <- apply ( exp.y, 2, mean ) 8 > y.975 <- apply ( exp.y, 2, quantile, probs = 0.975) 9 > y.025 <- apply ( exp.y, 2, quantile, probs = 0.025) 10 > y.995 <- apply ( exp.y, 2, quantile, probs = 0.995) 11 > y.005 <- apply ( exp.y, 2, quantile, probs = 0.005) 12 > 13 > plot (x, y, type = "p", ylim = c(0, 10), las = 1) 14 > lines ( new.x, y.mean, lty = 1) 15 > lines ( new.x, y.975, lty = 2) 16 > lines ( new.x, y.025, lty = 2) 17 > lines ( new.x, y.995, lty = 3) 18 > lines ( new.x, y.005, lty = 3) 19 > legend (" bottomright ", 20 + legend = c(" mean ", "95% interval ", "99% interval "), 21 + lty = c(1, 2, 3)) post.samp chain 1 2 unlist() X MCMC Y 9

24 y mean 95% 99% x Y Stan Stan R example2_stan.r RStan 1 > library ( rstan ) example2_code 1 > example2_ code <- " 2 + data { 3 + int < lower =0 > N; 4 + int < lower =0 > K; 5 + int < lower =0 > X[N]; 6 + int < lower =0 > Y[N]; 7 + } 8 + parameters { 9 + real beta ; 10 + real beta_ x ; 11 + } 12 + transformed parameters { 13 + real < lower =0, upper =1 > p[n]; 14 +

25 15 + for (i in 1:N) { 16 + p[i] <- inv_logit ( beta + beta_x * X[i ]); 17 + } 18 + } 19 + model { 20 + Y ~ binomial (K, p); # priors 23 + beta ~ normal (0.0, 1.0 e +3); 24 + beta_ x ~ normal (0.0, 1.0 e +3); 25 + } 26 + " Stan data parameters transformed parameters model. _ (;) chain 1 > n. chains <- 3 2 > inits <- vector (" list ", 3) 3 > inits [[1]] <- list ( beta = -10, beta_x = 0) 4 > inits [[2]] <- list ( beta = -5, beta_x = 2) 5 > inits [[3]] <- list ( beta = 0, beta_x = -2) 6 > pars <- c(" beta ", " beta_x ") stan() Stan fit 1 > fit <- stan ( model_ code = example2_ code, 2 + data = list (X = x, Y = y, N = n, K = k), 3 + pars = pars, init = inits, seed = 123, 4 + chains = n. chains, 5 + iter = 2500, warmup = 500, thin = 2) warmup burn-in [7] 1 > print ( fit ) 2 Inference for Stan model : b08b951f11fe4269ab8d08c9ec. 3 3 chains, each with iter =2500; warmup =500; thin =2; 4 post - warmup draws per chain =1000, total post - warmup draws = mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat 7 beta beta_x

26 9 lp Samples were drawn using NUTS ( diag_e ) at Wed Oct 7 14:10: For each parameter, n_eff is a crude measure of effective sample size, 13 and Rhat is the potential scale reduction factor on split chains ( at 14 convergence, Rhat =1). lp MCMC R MCMC WinBUGS OpenBUGS? ~ <-? ( () {} )? typo??? :? 2 : MCMC MCMC??? Markov chain

27 X 1, X 2, Y Y X 1 X 2 ( ) example3.csv 1 > data <- read. csv (" example3. csv ") 10 1 > head (data, 10) 2 block x1 x2 y > pairs ( data ) JAGS R example3.r BUGS 1 > n. block <- max ( data$block ) # number of blocks 2 > n. data <- nrow ( data ) / n. block # number of data in a block 3 >

28 block x1 x y > x1 <- t( matrix ( data$x1, nrow = n.data, ncol = n. block )) 5 > x2 <- t( matrix ( data$x2, nrow = n.data, ncol = n. block )) 6 > y <- t( matrix ( data$y, nrow = n.data, ncol = n. block )) x1 1 > print (x1) 2 [,1] [,2] [,3] [,4] [,5] 3 [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,]

29 ϵ Bi i σ B Y ij Normal(µ ij, σ 2 ) µ ij = β + β 1 X 1ij + β 2 X 2ij + ϵ Bi ϵ Bi Normal(0, σb) 2 µ ij i j σ β β 1 β 2 X 1 X 2 BUGS example3_model.txt 1 var 2 M, # Number of blocks 3 N, # Number of observations 4 X1[N], X2[N], # Data 5 Y[N], 6 B[N], # Block 7 e.b[m], # Random effect 8 beta, beta.1, beta.2, # Parameters 9 tau, sigma, 10 tau. B, sigma. B; # Hyperparameters model { 13 for (i in 1:N) { 14 Y[i] ~ dnorm (mu[i], tau ) 15 mu[i] <- beta + beta.1 * X1[i] + 16 beta.2 * X2[i] + e.b[b[i]] 17 } 18 for (i in 1:M) { 19 e.b[i] ~ dnorm (0, tau.b) 20 } ## Priors 23 beta ~ dnorm (0, 1.0 E -6) 24 beta.1 ~ dnorm (0, 1.0E -6) 25 beta.2 ~ dnorm (0, 1.0E -6) 26 tau <- 1 / ( sigma * sigma ) 27 tau.b <- 1 / ( sigma.b * sigma.b) 28 sigma ~ dunif (0, 1.0 E +4) 29 sigma.b ~ dunif (0, 1.0 E +4) 30 } BUGS

30 e.b[] 0 1/tau.B (tau.b = 1/sigma.B 2 ) 18 sigma.b [0, 10000] 28 tau.b(sigma.b) e.b[] (hyperparameter) 11 data X 1 X 2 Y process Pr(X 1, X 2,Y β,β 1,β 2,τ,e) parameter β β 1 β 2 τ e hyperparameter τ e 11 rjags JAGS bugs 2 1 > library ( rjags ) 2 > load. module (" glm ") 3 > 4 > ## Model file 5 > model. file <- " example3_model. txt " 6 > 7 > ## Number of chains

31 8 > n. chains <- 3 9 > 10 > ## Initial values 11 > inits <- vector (" list ", n. chains ) 12 > inits [[1]] <- list ( beta = 5, beta.1 = 0, beta.2 = 0, 13 + sigma = 1, sigma. B = 1, RNG. seed = 123, RNG. name = " base :: Mersenne - Twister ") 16 > inits [[2]] <- list ( beta = -5, beta.1 = 10, beta.2 = 10, 17 + sigma = 10, sigma. B = 10, RNG. seed = 1234, RNG. name = " base :: Mersenne - Twister ") 20 > inits [[3]] <- list ( beta = 0, beta.1 = -10, beta.2 = -10, 21 + sigma = 5, sigma. B = 5, RNG. seed = 12345, RNG. name = " base :: Mersenne - Twister ") 24 > 25 > ## Parameters 26 > pars <- c(" beta ", " beta.1", " beta.2", 27 + " sigma ", " sigma.b", "e.b") 28 > 29 > ## MCMC 30 > model <- jags. model ( file = model. file, 31 + data = list (M = n. block, N = n.data, 32 + X1 = x1, X2 = x2, Y = y), 33 + inits = inits, n. chains = n. chains, 34 + n. adapt = 1000) % 36 > 37 > ## Burn -in 38 > update ( model, n. iter = 1000) 39 ************************************************** 100% 40 > 41 > ## Sampling 42 > post. samp <- coda. samples ( model, n. iter = 5000, thin = 5, 43 + variable. names = pars ) 44 ************************************************** 100% 1 > gelman. diag ( post. samp ) 2 Potential scale reduction factors : 3 4 Point est. Upper C.I. 5 beta

32 6 beta beta e.b [1] e.b [2] e.b [3] e.b [4] e.b [5] e.b [6] e.b [7] e.b [8] sigma sigma.b Multivariate psrf > summary ( post. samp ) Iterations = 2005: Thinning interval = 5 26 Number of chains = 3 27 Sample size per chain = Empirical mean and standard deviation for each variable, 30 plus standard error of the mean : Mean SD Naive SE Time - series SE 33 beta beta beta e. B [1] e. B [2] e. B [3] e. B [4] e. B [5] e. B [6] e. B [7] e. B [8] sigma sigma. B Quantiles for each variable : % 25% 50% 75% 97. 5% 50 beta

33 51 beta beta e. B [1] e. B [2] e. B [3] e. B [4] e. B [5] e.b [6] e. B [7] e. B [8] sigma sigma. B e.b[] 12 Posterior distribution of e.b[] density value 12 e.b[] Nested Indexing 2 Nested Indexing BUGS example3-1_model.txt

34 1 > data <- read. csv (" example3. csv ") 2 > n. block <- 8 # number of blocks 3 > n. row <- nrow ( data ) # number of observations data$block 1 > data$block 2 [1] [27] X1 X2 Y B 6 e.b B 16 1 var 2 M, # Number of blocks 3 N, # Number of observations 4 X1[N], X2[N], # Data 5 Y[N], 6 B[N], # Block 7 e.b[m], # Random effect 8 beta, beta.1, beta.2, # Parameters 9 tau, sigma, 10 tau. B, sigma. B; # Hyperparameters model { 13 for (i in 1:N) { 14 Y[i] ~ dnorm (mu[i], tau ) 15 mu[i] <- beta + beta.1 * X1[i] + 16 beta.2 * X2[i] + e.b[b[i]] 17 } 18 for (i in 1:M) { 19 e.b[i] ~ dnorm (0, tau.b) 20 } ## Priors 23 beta ~ dnorm (0, 1.0 E -6) 24 beta.1 ~ dnorm (0, 1.0E -6) 25 beta.2 ~ dnorm (0, 1.0E -6) 26 tau <- 1 / ( sigma * sigma ) 27 tau.b <- 1 / ( sigma.b * sigma.b) 28 sigma ~ dunif (0, 1.0 E +4) 29 sigma.b ~ dunif (0, 1.0 E +4)

35 30 } jags.model() data 1 model <- jags. model ( file = model.file, 2 data = list (M = n. block, N = n.data, 3 X1 = data$x1, X2 = data$x2, 4 Y = data$y, B = data$block ), 5 inits = inits, n. chains = n. chains, 6 n. adapt = 1000) (centering) (centering) [23] Box beta.1 * X1[i] beta.1 * (X1[i] - X1.bar) (X1.bar X1 ) JAGS glm (JAGS User s Manual[25] 4.6 ) m 1m 36 example4.csv 1 " Plot "," Num "," Light " 2 1,0, ,0, ,3, : 6 36,1,0.412 Plot Num Light (%) 13 log GLM

36 5 4 3 Num Light 13 1 > summary ( glm ( Num ~ Light, family = poisson, data = data )) 2 3 Call : 4 glm ( formula = Num ~ Light, family = poisson, data = data ) 5 6 Deviance Residuals : 7 Min 1 Q Median 3 Q Max Coefficients : 11 Estimate Std. Error z value Pr ( > z ) 12 ( Intercept ) Light ( Dispersion parameter for poisson family taken to be 1) Null deviance : on 35 degrees of freedom 18 Residual deviance : on 34 degrees of freedom 19 AIC : Number of Fisher Scoring iterations : 6 (Residual deviance) (degree of freedom) 2

37 (overdispersion) 0 14 (Zero-inflated) 0 (Zero-Inflated Poisson model; ZIP model) Frequency Number of seedlings 14 ZIP BUGS 1 var 2 N, # Number of observations 3 Y[N], # Number of new seedlings 4 X[N], # Proportion of open canopy 5 lambda [N], # Poisson mean 6 z[ N], # 0: absent, 1: at least latently present 7 p, # Probability of the presence ( at least latently ) 8 beta, # Intercept in the linear model 9 beta. x; # Coefficient of X in the linear model 10 model { 11 for (i in 1:N) { 12 Y[i] ~ dpois ( lambda [i]) 13 lambda [i] <- z[i] * exp ( beta + beta.x * X[i]) 14 z[i] ~ dbern (p) 15 }

38 16 ## Priors 17 p ~ dunif (0, 1) 18 beta ~ dnorm (0, 1.0 E -4) 19 beta.x ~ dnorm (0, 1.0E -4) 20 } z = 0 2. z = 1 0 Poisson(0 λ) z 0 1 p λ log λ = β + β x X rjags burn-in > summary ( post. samp ) 2 3 Iterations = 2010: Thinning interval = 10 5 Number of chains = 3 6 Sample size per chain = Empirical mean and standard deviation for each variable, 9 plus standard error of the mean : Mean SD Naive SE Time - series SE 12 beta beta. x p Quantiles for each variable : % 25% 50% 75% 97. 5% 19 beta beta. x p p % beta beta.x % GLM 0.177

39 4 Models for Ecological Data [3], Hierarchical modelling for the environmental sciences [4], Introduction to WinBUGS for ecologists [12], Bayesian population analysis using WinBUGS [13], Bayesian methods for ccology [23] [17] GLM GLMM 2009 [15] [5] [34] [35] Hamiltonian Monte Carlo Vol.1 [11] JAGS Stan * 15 MCMC [1] Bishop C.M. (2006) Pattern recognition and machine learning. Springer-Verlag, New York. : (2012) /, [2] Brooks S., Gelman A., Jones G.L., Meng X.-L. (2011) Handbook of Markov chain Monte Carlo. Chapman & Hall/CRC, Boca Raton. [3] Clark J.S. (2007) Models for ecological data. Princeton University Press, Princeton. [4] Clark J.S., Gelfand A.E. (2006) Hierarchical modelling for the environmental sciences. Oxford University Press, New York. [5] (2009). 59: [6] (2008) R & WinBUGS.,. [7] Gelman A, Carlin J.B., Stern H.S., Dunson D.B., Vehtari A., Rubin D.B. (2014) Bayesian data analysis, 3rd ed. Chapman & Hall/CRC, Boca Raton. *15

40 [8] Gilks W.R., Richardson S.R., Spiegelhalter D.J. (eds.) (1996) Markov chain Monte Carlo in practice. Chapman & Hall/CRC, Boca Raton. [9] (2003).,. [10] (2005). ( ( ) II, ): [11] ( ) (2015) Vol.1.,. [12] Kéry M. (2010) Introduction to WinBUGS for ecologists: a Bayesian approach to regression, ANOVA, mixed models and related analyais. Academic Press, Waltham. [13] Kéry M, Schaub M. (2011) Bayesian population analysis using WinBUGS. Academic Press, Waltham. [14] Kruschke J. (2014) Doing Bayesian data analysis, 2nd ed.: a tutorial with R, JAGS, and Stan. Academic Press, Waltham. [15] (2009). 59: [16] (2009) [I]. 92: [17] (2012) MCMC.,. [18] (2010).,. [19] Link, W.A., Eaton, M.J. (2012) On thinning of chains in MCMC. Methods in Ecology and Evolution 3: doi: /j X x [20] Lunn D., Spiegelhalter D., Thomas A., Best N. (2009) The BUGS project: Evolution, critique, and future directions. Statistics in Medicine 28: [21] Lunn D., Jackson C, Besk N., Thomas A., Spiegelhalter D. (2012) The BUGS Book. Chapman & Hall/CRC, Boca Raton. [22] Martin A.D., Quinn, K.M. (2006) Applied Bayesian inference in R using MCMCpack. R News 6(1): [23] McCarthy M.A. (2007) Bayesian methods for ecology. Cambridge University Press, New York. : (2009),. [24] Ntzoufras I. (2009) Bayesian modeling using WinBUGS. Wiley, Hoboken. [25] Plummer M. (2015) JAGS version user manual. [26] R Core Team (2015) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna.

41 [27] Robert C.P., Casella G. (2010) Introducing Monte Carlo Methods with R. Springer, New York. : (2012) R,. [28] Spiegelhalter D., Tomas A., Best N., Lunn D. (2003) WinBUGS user manual version [29] Stan Development Team (2015) Stan Modeling Language: User s Guide and Reference Manual, Version [30] (2000).,. [31] (2010) R.,. [32] Thomas A. (2006) The BUGS language. R News 6(1): [33] Thomas A., O Hara B., Ligges U., Sturtz S. (2006) Making BUGS open. R News 6(1): [34] (2008).,. [35] ( ) (2015).,. [36] (2005). ( ( ), ): [37] (2012).,.

Vol. 15 No. 1 JOURNAL OF HARBIN UNIVERSITY OF SCIENCE AND TECHNOLOGY Feb O21 A

Vol. 15 No. 1 JOURNAL OF HARBIN UNIVERSITY OF SCIENCE AND TECHNOLOGY Feb O21 A 5 200 2 Vol 5 No JOURNAL OF HARBIN UNIVERSITY OF SCIENCE AND TECHNOLOGY Feb 200 2 2 50080 2 30024 O2 A 007-2683 200 0-0087- 05 A Goodness-of-fit Test Based on Empirical Likelihood and Application ZHOU