Poisson model:extra program for MH-step

Author: HWagn

vignettes/Chapter08.Rmd

@@ -646,9 +646,47 @@ parms.proposal <- gen.proposal.poisson(y, X, e, b0 = 0, B0 = 100)
 parms.proposal
 ```
 
-Next we set up the independence Metropolis-Hastings algorithm to
-estimate the model parameters.
+To set up the independence Metropolis-Hastings algorithm for 
+the Poisson model, we first write a short program for the MH step
+to sample the regression effects.
+```{r}
+sample_beta<- function(y,X,e, b0, B0, qmean, qvar, beta.old){
+ 
+     beta.proposed <- as.vector(mvtnorm::rmvnorm(1, mean = qmean, sigma = qvar))
+     
+     # Compute log proposal density at proposed and old value
+     lq_proposed <- mvtnorm::dmvnorm(beta.proposed, mean = qmean, sigma = qvar,
+                                log = TRUE)
+     lq_old  <- mvtnorm::dmvnorm(beta.old, mean = qmean, sigma = qvar,
+                            log = TRUE)
 
+     # Compute log prior  of proposed and old value
+     lpri_proposed <- mvtnorm::dmvnorm(beta.proposed, mean = b0, sigma = B0, 
+                                       log = TRUE)
+     lpri_old  <- mvtnorm::dmvnorm(beta.old,  mean = b0, sigma = B0, 
+                                   log = TRUE)
+     # Compute log likelihood of proposed and old value
+     lh_proposed <- dpois(y, e * exp(X %*% beta.proposed), log = TRUE)
+     lh_old  <- dpois(y, e * exp(X %*% beta.old), log = TRUE)
+
+     maxlik <- max(lh_old, lh_proposed)
+     ll <- sum(lh_proposed - maxlik) - sum(lh_old - maxlik)
+    
+     # Compute acceptance probability and accept or not
+     log_acc <- min(0, ll + lpri_proposed - lpri_old + lq_old - lq_proposed)
+
+     if (log(runif(1)) < log_acc) {
+        beta <- beta.proposed
+        acc  <- 1
+     } else {
+        beta <- beta.old
+        acc  <- 0
+     }
+     return(res = list(beta=beta, acc=acc))
+}
+```
+
+We use this program to sample from the posterior.
 ```{r}
 poisson <- function(y, X, e, b0 = 0, B0 = 100, qmean, qvar,
                     burnin = 1000L, M = 10000L) {
@@ -664,43 +702,13 @@ poisson <- function(y, X, e, b0 = 0, B0 = 100, qmean, qvar,
   beta <- as.vector(mvtnorm::rmvnorm(1, mean = qmean, sigma = qvar))
   
   for (m in seq_len(burnin + M)) {
-    beta.old <- beta
-    beta.proposed <- as.vector(mvtnorm::rmvnorm(1, mean = qmean, sigma = qvar))
-    
-    # Compute log proposal density at proposed and old value
-    lq_proposed <- mvtnorm::dmvnorm(beta.proposed, mean = qmean, sigma = qvar,
-                                    log = TRUE)
-    lq_old  <- mvtnorm::dmvnorm(beta.old, mean = qmean, sigma = qvar,
-                                log = TRUE)
-    
-    # Compute log prior of proposed and old value
-    lpri_proposed <- mvtnorm::dmvnorm(beta.proposed, mean = b0, sigma = B0,
-                                      log = TRUE)
-    lpri_old  <- mvtnorm::dmvnorm(beta.old, mean = b0, sigma = B0,
-                                  log = TRUE)
-    
-    # Compute log-likelihood of proposed and old value
-    lh_proposed <- dpois(y, e * exp(X %*% beta.proposed), log = TRUE)
-    lh_old  <- dpois(y, e * exp(X %*% beta.old), log = TRUE)
-    
-    maxlik <- max(lh_old, lh_proposed)
-    ll <- sum(lh_proposed - maxlik) - sum(lh_old - maxlik)
-    
-    # Compute acceptance probability and accept or not
-    log_acc <- min(0, ll + lpri_proposed - lpri_old + lq_old - lq_proposed)
-    
-    if (log(runif(1)) < log_acc) {
-      beta <- beta.proposed
-      accept <- 1
-    } else {
-      beta <- beta.old
-      accept <- 0
-    }
-
-    # Store the beta draws
+   beta.draw <- sample_beta(y,X,e, b0, B0, qmean, qvar,beta)
+   beta<- beta.draw$beta
+   
+   # Store the beta draws
     if (m > burnin) {
         beta.post[m - burnin, ] <- beta
-        acc[m - burnin] <- accept
+        acc[m - burnin] <- beta.draw $acc
     }
   }
   return(list(beta.post = beta.post, accept = mean(acc)))