changes negbin

Author: HWagn

vignettes/Chapter08.Rmd

@@ -606,7 +606,10 @@ gen_proposal_poisson <- function(y, X, e, b0 = 0, B0 = 100, t.max = 20) {
   d <- ncol(X)
 
   beta_new <- matrix(c(log(mean(y/e)),rnorm(d - 1)/10), nrow = d)
-
+  if(is.nan(beta_new[1])){
+    beta.new[1]<- log(mean(y)/mean(e))
+  }
+  
   XtXi <- lapply(seq_len(N), function(i) tcrossprod(X[i,]))
   B0.inv <- solve(B0)
   for (t in seq_len(t.max)) {
@@ -919,7 +922,7 @@ b0=rep(0,d)
 B0=diag(100,d)
 
 pri_alpha <- data.frame(shape = 2, rate = 0.5)
-c_alpha=0.1
+c_alpha=0.3
                                                     
 M <- 50000L / mcmcspeedup
                                                     
@@ -970,7 +973,7 @@ sampler.
 
 ## Section 8.2.3: Evaluating MCMC samplers
 
-### Example 8.10 Verifying the correctness of the full conditional MCMC samper
+### Example 8.10 Verifying the correctness of the full conditional MCMC sampler
 
 To check the MCMC algorithm for correctness, we extend the sampler by adding sampling the data from the prior as a further sampling step.
 
@@ -1392,11 +1395,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 <- 300000 / (10 * mcmcspeedup)
+M <- 1000000 / (10 * mcmcspeedup)
 set.seed(1234)
 res_check_full <- negbin_check_cba(X, e, b0, B0, pri_alpha, c_alpha,
                                   full_gibbs = TRUE, M = M)
-h=300
+h=1000
 thin=seq(from=1, to=M,by=h)
 
 mu1=exp(res_check_full$beta_post[,1])
@@ -1461,10 +1464,9 @@ 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(5,0.1, paste0('KS-test: p-value= ', round(ks2$p.value,4))) 
-
 ```
 
-We next run the ivalid partially marginalised Gibbs sampler.
+We next run the invalid partially marginalised Gibbs sampler.
 
 ```{r, echo = -c(1:3)}
 if (pdfplots) {
@@ -1495,7 +1497,6 @@ 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(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. 
@@ -1506,7 +1507,16 @@ versus the heterogeneity parameter $\alpha$
 if (pdfplots) {
   pdf("8-2_7.pdf", width = 8, height = 4)
 }
-par(mfrow = c(2,3), mar = c(2.5, 2.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
+par(mfrow = c(3,3), mar = c(2.5, 2.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
+
+plot(res_check_full$beta_post[thin,1],res_check_full$beta_post[thin,2],
+      xlab="beta1", ylab="beta2",main="full",xlim=c(-1,2),
+     ylim=c(0.5,3.5))
+plot(res_partial$beta_post[thin,1],res_partial$beta_post[thin,2],
+           xlab="beta1", ylab="beta2",main="correct partial",xlim=c(-1,2), ylim=c(0.5,3.5))
+plot(res_check_partial$beta_post[thin,1],res_check_partial$beta_post[thin,2],
+           xlab="beta1", ylab="beta2",main="incorrect partial",xlim=c(-1,2),
+     ylim=c(0.5,3.5))
 
 plot(res_check_full$beta_post[thin,1],res_check_full$alpha_post[thin],
       xlab="beta1", ylab="alpha",main="full",xlim=c(-1,2), ylim=c(0,1.4))
@@ -1516,11 +1526,11 @@ plot(res_check_partial$beta_post[thin,1],res_check_partial$alpha_post[thin],
            xlab="beta1", ylab="alpha",main="incorrect partial",xlim=c(-1,2), ylim=c(0,1.4))
 
 plot(res_check_full$beta_post[thin,2],res_check_full$alpha_post[thin],
-        xlab="beta2", ylab="alpha",xlim=c(0.5,3.5), ylim=c(0,1.4))
+        xlab="beta2", ylab="alpha",main="full",xlim=c(0.5,3.5), ylim=c(0,1.4))
 plot(res_partial$beta_post[thin,2],res_partial$alpha_post[thin],
-      xlab="beta2", ylab="alpha",xlim=c(0.5,3.5), ylim=c(0,1.4))
+      xlab="beta2", ylab="alpha",main="correct partial",xlim=c(0.5,3.5), ylim=c(0,1.4))
 plot(res_check_partial$beta_post[thin,2],res_check_partial$alpha_post[thin],
-        xlab="beta2", ylab="alpha",xlim=c(0.5,3.5), ylim=c(0,1.4))
+        xlab="beta2", ylab="alpha",main="incorrect partial",xlim=c(0.5,3.5), ylim=c(0,1.4))
 ```
 # Section 8.3: Beyond i.i.d. Gaussian error distributions