change in proposal poisson and minor changes

Author: HWagn

vignettes/Chapter08.Rmd

@@ -604,7 +604,7 @@ gen_proposal_poisson <- function(y, X, e, b0 = 0, B0 = 100, t.max = 20) {
   N <- length(y)
   d <- ncol(X)
 
-  beta_new <- matrix(c(log(mean(y/e)), rep(0, d - 1)), nrow = d)
+  beta_new <- matrix(c(log(mean(y/e)),rnorm(d - 1)/10), nrow = d)
 
   XtXi <- lapply(seq_len(N), function(i) tcrossprod(X[i,]))
   B0.inv <- solve(B0)
@@ -723,7 +723,7 @@ poisson <- function(y, X, e, b0 = 0, B0 = 100, burnin = 1000L, M = 10000L) {
 }
 ```
 
-We perform MCMC and report the results.
+We perform MCMC for !0000 iterations after a burnin of 1000 and report the results.
 
 ```{r}
 set.seed(1234)
@@ -776,7 +776,7 @@ knitr::kable(res_poisson2)
 With 2.8 deadly or seriously injured children per 10000 at risk, the
 estimated baseline risk is very similar to that from Model 1. Also the
 estimated intervention effect is very similar in both models,
-indicating a reduction of the risk by a factor of 0.69 in Model 2
+indicating a reduction of the risk by a factor of 0.692 in Model 2
 (compared to 0.697 in Model 1). The monthly effects have rather wide
 95\% HPD intervals that cover 0 for all months except for July and
 August. For these two holiday months they are clearly negative,
@@ -908,7 +908,7 @@ negbin<- function(y,X,e, b0,B0, pri_alpha,c_alpha,
 
 We use the same Normal prior as in the Poisson model for the
 regression effects $\boldsymbol{\beta}$ and a Gamma prior
-$\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.
+$\mathcal{G}(2 , 0.5)$ for $\alpha$. We first run the full Gibbs sampler for $M=50,000$ iterations after a burn-in of 1000.
 
 ```{r,negbin}
 d <- ncol(X)
@@ -966,14 +966,14 @@ sampler.
 
 
 ## Section 8.2.3: Evaluating MCMC samplers
+
 ### Example 8.10 Verifying the correctness of the full conditional MCMC samper
 
 To check the MCMC algorithm for correctness, we extend the sampler by adding sampling the data from the prior as a further sampling step.
 
 ```{r}
 negbin_check <- function(X, e, b0, B0, pri_alpha, c_alpha,
-                             full_gibbs = TRUE, burnin = 1000L, 
-                             M = 50000L) {
+                          full_gibbs = TRUE, burnin = 1000L,   M = 50000L) {
     
     N <- nrow(X)
     d <- ncol(X)
@@ -1338,13 +1338,13 @@ negbin_check_cba <- function(X, e, b0, B0, pri_alpha, c_alpha,
     alpha_draw <- sample_alpha(y, e * exp(linpred), phi, pri_alpha, alpha,
                                c_alpha, full_gibbs)
     alpha <- alpha_draw$alpha
-    
+  
     # Step a
     parms_proposal <- gen_proposal_poisson(y, X, e * phi, b0, B0)
     beta_draw <- sample_beta(y, X, e * phi, b0, B0, parms_proposal$mean, 
                              parms_proposal$var, beta)
     beta <- beta_draw$beta
-    
+  
     # Save the draws
     if (m > burnin) {
       beta_post[m - burnin, ] <- beta
@@ -1369,14 +1369,15 @@ X <- cbind(rep(1, N), c(rep(0,N/2),rep(1,N/2)))
 e <- rep(1, N)
 
 d <- ncol(X)
-b0 <- c(0,1)
-B0 <- diag(0.25, d)
+b0 <- c(0.5,2)
+B0 <- diag(0.2, d)
 
-pri_alpha <- data.frame(shape = 1.5, rate = 3)
+pri_alpha <- data.frame(shape = 5, rate = 10)
 c_alpha <- 0.35
 ```
 We run the full Gibbs sampler under this scheme. To check the correctness of the
-sampler we focus on the overdispersion  $\frac{\mu_i^2}{\alpha}$  computed from 
+sampler we focus on the overdispersion  $\frac{\mu_i^2}{\alpha}$ for $x_i=0$ and
+$x-i=1$ which we compute from 
 the draws of the augmented MCMC sampler as well as from draws of the prior
 distribution.
 
@@ -1386,11 +1387,11 @@ if (pdfplots) {
 }
 par(mfrow = c(1, 2), mar = c(2.5, 2.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
 
-M=20000
+M=30000
 set.seed(1234)
 res_check_full <- negbin_check_cba(X, e, b0, B0, pri_alpha, c_alpha,
                                   full_gibbs = TRUE, M = M)
-h=200
+h=300
 thin=seq(from=1, to=M,by=h)
 
 mu1=exp(res_check_full$beta_post[,1])
@@ -1404,21 +1405,26 @@ print(coda::effectiveSize(ov2))
 set.seed(1234)
 beta_prior <- t(mvtnorm::rmvnorm(M/h, mean = b0, sigma = B0))
 alpha_prior <- rgamma(M/h, shape = pri_alpha$shape,rate = pri_alpha$rate)
-ov1_pri<- exp(beta_prior[1,])^2/alpha_prior
-ov2_pri<- exp(beta_prior[1,]+beta_prior[2,])^2/alpha_prior
+ov1_prior<- exp(beta_prior[1,])^2/alpha_prior
+ov2_prior<- exp(beta_prior[1,]+beta_prior[2,])^2/alpha_prior
+print(coda::effectiveSize(ov2_prior))
+print(coda::effectiveSize(ov1_prior))
 
-ks1<- ks.test(ov1_pri,ov1)
-qqplot(ov1_pri, ov1,xlab = "Prior",xlim=c(0,50), ylim=c(0,50),
+ks1<- ks.test(ov1_prior,ov1)
+qqplot(log(ov1_prior), log(ov1),xlab = "Prior",xlim=c(0,6), ylim=c(0,6),
        ylab = "Posterior", main = "Overdispersion for X=0")
 abline(a = 0, b = 1)
-text(30,0.1, paste0('KS-test: p-value= ', round(ks1$p.value,4)))
+text(3,0.1, paste0('KS-test: p-value= ', round(ks1$p.value,4)))
 
-ks2<- ks.test(ov2_pri,ov2)
-qqplot(ov2_pri,ov2, xlab = "Prior",xlim=c(0,300), ylim=c(0,300),
+ks2<- ks.test(ov2_prior,ov2)
+qqplot(log(ov2_prior),log(ov2), xlab = "Prior",xlim=c(0,10), ylim=c(0,10),
        ylab = "Posterior", main = "Overdispersion for X=1")
 abline(a = 0, b = 1)
-text(200,0.1, paste0('KS-test: p-value= ', round(ks2$p.value,4)))
+text(5,0.1, paste0('KS-test: p-value= ', round(ks2$p.value,4)))
 ```
+Both p_values are larger than 0.05, which is as expected for the full-conditional
+Gibbs sampler.
+
 Now we check the partially marginalised Gibbs sampler.
 ```{r, echo = -c(1:3)}
 if (pdfplots) {
@@ -1431,26 +1437,29 @@ res_check_partial <- negbin_check_cba(X, e, b0, B0, pri_alpha, c_alpha,
                                   full_gibbs = FALSE, M = M)
 
 mu1=exp(res_check_partial$beta_post[,1])
-ov1<-(mu1^2/res_check_partial$alpha_post)[thin]
+ov1<-((mu1^2)/res_check_partial$alpha_post)[thin]
 print(coda::effectiveSize(ov1))
 
 mu2=exp(res_check_partial$beta_post[,1]+res_check_partial$beta_post[,2])
-ov2<-(mu2^2/res_check_partial$alpha_post)[thin]
+ov2<-((mu2^2)/res_check_partial$alpha_post)[thin]
 print(coda::effectiveSize(ov2))
 
-ks1<- ks.test(ov1_pri,ov1)
-qqplot(ov1_pri, ov1,xlab = "Prior",xlim=c(0,30), ylim=c(0,30),
+ks1<- ks.test(ov1_prior,ov1)
+qqplot(log(ov1_prior), log(ov1),xlab = "Prior",xlim=c(0,5), ylim=c(0,5),
        ylab = "Posterior", main = "Overdispersion for X=0")
 abline(a = 0, b = 1)
-text(20,0.1, paste0('KS-test: p-value= ', round(ks1$p.value,4)))
+text(3,0.1, paste0('KS-test: p-value= ', round(ks1$p.value,4)))
 
-ks2<- ks.test(ov2_pri,ov2)
-qqplot(ov2_pri,ov2, xlab = "Prior",xlim=c(0,300), ylim=c(0,300),
+ks2<- ks.test(ov2_prior,ov2)
+qqplot(log(ov2_prior),log(ov2), xlab = "Prior",xlim=c(0,10), ylim=c(0,10),
        ylab = "Posterior", main = "Overdispersion for X=1")
 abline(a = 0, b = 1)
-text(200,0.1, paste0('KS-test: p-value= ', round(ks2$p.value,4))) 
+text(5,0.1, paste0('KS-test: p-value= ', round(ks2$p.value,4))) 
 
 ```
+Also for the partially marginalised Gibbs sampler both p-values are larger than
+0.05 and hence we fail to detect that this sampler is wrong. 
+
 
 # Section 8.3: Beyond i.i.d. Gaussian error distributions