@@ -408,7 +408,7 @@ posteriors in case of separation and with a tighter prior we can
408408shrink coefficients to zero.
409409
410410
411- # Example 8.6: Complete separation: analysis under an informative prior
411+ ### Example 8.6: Complete separation: analysis under an informative prior
412412
413413We now analyze the data of example 8.4. under the more informative
414414prior $\mathcal{N}(\mathbf{0}, \mathbf{I})$. This prior distribution
@@ -910,7 +910,7 @@ We use the same Normal prior as in the Poisson model for the
910910regression effects $\boldsymbol{\beta}$ and a Gamma prior
911911$\mathcal{G}(2 , 0.5)$ for $\alpha$. We first run the full Gibbs sampler for $M=10,000$ iterations after a burn-in of 1000.
912912
913- ``` {r}
913+ ``` {r,negbin }
914914d <- ncol(X)
915915b0=rep(0,d)
916916B0=diag(100,d)
@@ -1229,9 +1229,9 @@ text(4, 0.2, paste0('KS-test: p-value= ', round(ks.a$p.value,4)))
12291229We see very a clear deviation from the identity line and a p-value of 0
12301230for the heterogeneity parameter which indicates that the sampler is erroneous.
12311231
1232- ## Example 8.11
1232+ ## Example 8.11 Verifying the correctness of the partial marginalized Gibbs sampler
12331233
1234- We now analyze the partial marginalized Gibbs sampler
1234+ We now analyze the partial marginalized Gibbs sampler.
12351235
12361236``` {r, echo = -c(1:3)}
12371237if (pdfplots) {
@@ -1263,7 +1263,7 @@ abline(a = 0, b = 1)
12631263text(4, 0.2, paste0('KS-test: p-value= ', round(ks.a$p.value,4)))
12641264```
12651265
1266- Again the points of the QQplot lie close to the
1266+ Again the points of the QQ-plot lie close to the
12671267identity line and the p-value of the Kolmogorov-Smirnov test is larger than
126812680.05 and thus we conclude that the sampler is correct.
12691269
@@ -1296,28 +1296,26 @@ qqplot(alpha_prior,res_check_wrong$alpha_post[thin], xlab = "Prior",
12961296abline(a = 0, b = 1)
12971297text(4, 0.2, paste0('KS-test: p-value= ', round(ks.a$p.value,4)))
12981298```
1299- Also for the partial marginalized sampler with the wrong sampling step for the
1300- heterogeneity parameter we get a Q-Q plot which deviates considerably from the
1301- identity line and a p-value of the Kolmogorov-Smirnov statistics of 0, which
1302- indicates that the sampler is not correct.
1299+ This is not the case for the partial marginalized sampler with the wrong
1300+ sampling step for the heterogeneity parameter: we get a Q-Q plot which deviates
1301+ considerably from the lidentity line and a p-value of the Kolmogorov-Smirnov
1302+ statistics of 0, which indicates that the sampler is not correct.
13031303
1304- ## Example 8.12
1305- In this example we would like to check whether changing the order of the sampling
1304+ ## Example 8.12 : Changing the order of the sampling steps
1305+ In this example we check whether changing the order of the sampling
13061306steps to (c)- (b)-(a) (instead of (a)-(b)-(c) ) still yields a correct Gibbs
13071307sampler.
13081308We first set up the sampler.
13091309
13101310``` {r}
13111311negbin_check_cba <- function(X, e, b0, B0, pri_alpha, c_alpha,
1312- full_gibbs = TRUE, burnin = 1000L,
1313- M = 50000L) {
1314-
1312+ full_gibbs = TRUE, burnin = 1000L, M = 50000L) {
13151313 N <- nrow(X)
13161314 d <- ncol(X)
13171315
13181316 beta_post <- matrix(ncol = d, nrow = M)
13191317 colnames(beta_post) <- colnames(X)
1320- acc_beta <- numeric(length = M)
1318+ acc_beta <- rep(NA_real_, M)
13211319
13221320 alpha_post <- rep(NA_real_, M)
13231321 acc_alpha <- rep(NA_real_, M)
@@ -1331,7 +1329,7 @@ negbin_check_cba <- function(X, e, b0, B0, pri_alpha, c_alpha,
13311329
13321330 # sample new data
13331331 linpred <- X %*% beta
1334- y <- rpois(N, phi * e * exp(linpred))
1332+ y <- rpois(N, phi * e * exp(linpred))
13351333
13361334 # Step c
13371335 phi <- rgamma(N, shape = alpha + y, rate = alpha + e * exp(linpred))
@@ -1371,15 +1369,16 @@ X <- cbind(rep(1, N), c(rep(0,N/2),rep(1,N/2)))
13711369e <- rep(1, N)
13721370
13731371d <- ncol(X)
1374- b0 <- c(0,0.5 )
1372+ b0 <- c(0,1 )
13751373B0 <- diag(0.25, d)
13761374
1377- pri_alpha <- data.frame(shape = 2 , rate = 2 )
1375+ pri_alpha <- data.frame(shape = 1.5 , rate = 3 )
13781376c_alpha <- 0.35
13791377```
1380- We run the full Gibbs sampler under this scheme. To check the correctness of the sampler
1381- we focus on the overdispersion $\frac{\mu_i^2}{\alpha}$$ which we compute from the draws of
1382- the augmented MCMC sampler as well as from draws of the prior distribution.
1378+ We run the full Gibbs sampler under this scheme. To check the correctness of the
1379+ sampler we focus on the overdispersion $\frac{\mu_i^2}{\alpha}$ computed from
1380+ the draws of the augmented MCMC sampler as well as from draws of the prior
1381+ distribution.
13831382
13841383``` {r, echo = -c(1:3)}
13851384if (pdfplots) {
0 commit comments