Update Chapter08

Author: HWagn

R/.gitignore

@@ -0,0 +1 @@
+rtruncnorm.R

data/.gitignore

@@ -0,0 +1 @@
+labor.rda

vignettes/Chapter08_b.Rmd

@@ -0,0 +1,97 @@
+---
+title: "Chapter 8: Bayesian Learning Beyond Standard Regression Analysis"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteEncoding{UTF-8}
+  %\VignetteIndexEntry{Chapter 8: Bayesian Learning Beyond Standard Regression Analysis}
+  %\VignetteEngine{knitr::rmarkdown}
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.width = 8,
+  fig.height = 6
+)
+knitr::opts_knit$set(global.par = TRUE)
+pdfplots <- FALSE # default: FALSE; set this to TRUE only if you like pdf figures
+```
+
+```{r, include = FALSE}
+par(mgp = c(1.6, .6, 0), mar = c(2.6, 2.6, 2.6, .4), lwd = 1)
+```
+
+# Section 8.1
+We illustrate probit regression analysis for the labor market data.
+```{r}
+library("BayesianLearningCode")
+data("labor", package = "BayesianLearningCode")
+```
+
+We model the binary variable unemployment and use as covariates the variables female
+(binary), wcollar (binary) and age18 (quantitative, centered at 18 years).
+
+```{r}
+y <- labor$unemp98
+N <- length(y)  # number of observations 
+
+X <- cbind("intercept" = rep(1, N), "female"=labor$female, "wcollar"=labor$wcollar,
+            "age18"=labor$age-18) # regressor matrix
+
+d <- dim(X)[2] # number regression effects 
+p <- d - 1 # number of regression effects without intercept
+```
+We specify  the  prior on the regression effects as a rather flat multivariate
+Normal.
+
+```{r}
+# define prior parameters 
+B0.inv <- diag(rep(1 / 10000, d), nrow = d) 
+b0 <- rep(0, d) 
+B0inv.b0 <- iB0%*%b0
+```
+The regression coefficients are estimated using  data augmentation and Gibbs 
+sampling.
+```{r}
+set.seed(1)
+
+#define burnin and M
+burnin <- 1000
+M <- 10000
+
+
+# define quantities for the Gibbs sampler
+XX <- crossprod(X)
+BN <- solve(B0.inv + XX)
+
+ind0=(y==0) # indicators for zeros and ones
+ind1=(y==1)
+
+# starting values
+beta <- rep(0,k)
+z <- rep(NA_real_,n)
+
+
+for (m in 1:(burnin + M)) {
+    
+  # Draw z conditional on y and beta
+    u<- runif(N)
+    eta <- X %*% beta
+    pi<- pnorm(eta) 
+   
+    z[ind0] <- eta[ind0] + qnorm(u[ind0]*(1-pi[ind0]))
+    z[ind1] <- eta[ind1] + qnorm (1-u[ind1]*pi[ind1])
+  
+    # sample beta from the full conditional 
+    bN <- BN %*% (B0inv.b0 + t(X) %*% z)
+    beta <- t(mvtnorm::rmvnorm(1, mean = bN, sigma = BN))
+    
+    # Store the beta draws
+    betas[m, ] <- beta
+    
+}
+```
+