Merge branch 'main' of github.com:gregorkastner/BayesianLearningCode

Author: bettinagruen

vignettes/Chapter08.Rmd

@@ -17,7 +17,7 @@ knitr::opts_chunk$set(
   fig.height = 6
 )
 knitr::opts_knit$set(global.par = TRUE)
-pdfplots <- FALSE # default: FALSE; set this to TRUE only if you like pdf figures
+pdfplots <- FALSE# default: FALSE; set this to TRUE only if you like pdf figures
 ```
 
 ```{r, include = FALSE}
@@ -816,7 +816,7 @@ Gibbs sampler and the partially marginalised Gibbs sampler using a log random wa
 sample_alpha <- function(y, linpred,e, phi,  pri.alpha,alpha.old, 
                        c.alpha,full.gibbs){
  
-    alpha.proposed <- exp(rnorm(1,log(alpha.old),c.alpha ))
+    alpha.proposed <- exp(rnorm(1,log(alpha.old),c.alpha))
 
    if (full.gibbs) {
       llik_alpha.proposed <- sum(dgamma(phi, shape = alpha.proposed,
@@ -904,7 +904,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=50,000$ iterations after a burn-in of 1000.
+$\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.
 
 ```{r}
 d <- ncol(X)
@@ -914,7 +914,7 @@ B0=diag(100,d)
 pri.alpha <- data.frame(shape = 2, rate = 0.5)
 c.alpha=0.1
 
-M <- 10000L
+M <- 50000L
 
 # Full Gibbs sampler
 set.seed(1234)
@@ -985,11 +985,13 @@ negbin_check_abc <- function(X,e, b0,B0, pri.alpha,c.alpha,
   beta <- as.vector(mvtnorm::rmvnorm(1, mean = b0, sigma = B0))
   alpha <- pri.alpha$shape/pri.alpha$rate
   phi <- rgamma(N, shape = alpha , rate = alpha)
+  linpred=X%*%beta
   
   for (m in seq_len(burnin + M)){
     
     # sample new data
-    y <-  rpois(N, phi*e * exp(X%*%beta))
+    y <-rnbinom(N, alpha, mu=e* exp(linpred))  
+    #rpois(N, phi*e * exp(X%*%beta))
     
     # Step a
     parms.proposal <- gen.proposal.poisson(y, X, e*phi, b0, B0)
@@ -1020,10 +1022,19 @@ negbin_check_abc <- function(X,e, b0,B0, pri.alpha,c.alpha,
 }
 ```
 
-We use a tighter prior for the regression effects.
+We use a tighter prior for the regression effects and a sample size of  N=200 observations.
 
 ```{r}
-B0=diag(1,d)
+N <- 200
+X <- cbind(rep(1,N), rnorm(N), rnorm(N))
+e <- rep(1,N)
+d <- ncol(X)
+
+b0 <- rep(0,d)
+B0 <- diag(0.25,d)
+
+pri.alpha <- data.frame(shape = 4, rate =2)
+
 ```
 
 We  run the sampler and investigate the draws of intercept and heterogeneity
@@ -1035,7 +1046,7 @@ 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)
 
-set.seed(123)
+set.seed(1234)
 res_check_abc <- negbin_check_abc(X,e, b0,B0, pri.alpha,c.alpha,
                              full.gibbs = TRUE, M = M)
 
@@ -1048,7 +1059,7 @@ abline(a = 0, b = 1)
 
 qqplot(alpha.prior,res_check_abc$alpha.post, xlab = "Prior",
        ylab = "Posterior", main = "Heterogeneity parameter", 
-       xlim = c(0, 30), ylim = c(0, 30))
+       xlim = c(0, 12), ylim = c(0, 12))
 abline(a = 0, b = 1)
 ```
 
@@ -1080,7 +1091,8 @@ negbin_check_cba <- function(X,e, b0,B0,pri.alpha,c.alpha,
 
     # sample new data
     linpred <- X%*%beta
-    y <-  rpois(N, phi*e * exp(linpred))
+    y <- rnbinom(N, alpha, mu=e* exp(X%*%beta))  
+      #rpois(N, phi*e * exp(linpred))
 
     # Step c
     phi <- rgamma(N, shape = alpha + y, rate = alpha + e * exp(linpred))
@@ -1104,7 +1116,6 @@ negbin_check_cba <- function(X,e, b0,B0,pri.alpha,c.alpha,
       alpha.post[m - burnin] <- alpha
       acc.alpha[m - burnin] <- alpha.draw$acc
     }
-
   }
   return(res = list(beta.post = beta.post, acc.beta = acc.beta,
                     alpha.post = alpha.post, acc.alpha = acc.alpha))
@@ -1120,7 +1131,7 @@ 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)
 
-set.seed(123)
+set.seed(1234)
 res_check_cba <- negbin_check_cba(X,e, b0,B0,pri.alpha,c.alpha,
                                    full.gibbs = TRUE, M=M)
 
@@ -1130,7 +1141,7 @@ abline(a = 0, b = 1)
 
 qqplot(alpha.prior,res_check_cba$alpha.post, xlab = "Prior",
        ylab = "Posterior", main =  "Heterogeneity parameter",
-       xlim = c(0, 30), ylim = c(0, 30))
+       xlim = c(0, 12), ylim = c(0, 12))
 abline(a = 0, b = 1)
 ```
 
@@ -1144,7 +1155,7 @@ 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)
 
-set.seed(123)
+set.seed(1234)
 # order (a)-(b)-(c)
 res_check_abc <- negbin_check_abc(X,e, b0,B0,pri.alpha,c.alpha,
                              full.gibbs = FALSE, M = M)
@@ -1155,7 +1166,7 @@ abline(a = 0, b = 1)
 
 qqplot(alpha.prior,res_check_abc$alpha.post, xlab = "Prior", 
        ylab = "Posterior", main  = "Heterogeneity parameter", 
-       xlim = c(0, 30), ylim = c(0, 30))
+       xlim = c(0, 12), ylim = c(0, 12))
 abline(a = 0, b = 1)
 ```
 
@@ -1167,7 +1178,7 @@ 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)
 
-set.seed(123)
+set.seed(1234)
 # order (c)- (b)-(a)
 res_check_cba <- negbin_check_cba(X,e, b0,B0,pri.alpha,c.alpha,
                              full.gibbs = FALSE,M = M)
@@ -1178,7 +1189,7 @@ abline(a = 0, b = 1)
 
 qqplot(alpha.prior,res_check_cba$alpha.post, xlab = "Prior", 
        ylab = "Posterior", main = "Heterogeneity parameter",
-       xlim = c(0, 30), ylim = c(0, 30))
+       xlim = c(0, 12), ylim = c(0, 12))
 abline(a = 0, b = 1)
 ```