4.8

d)

Set \(a_\theta = 2\) and \(b_\theta\) = 1. Let \(a_\gamma\)=\(b_\gamma \in \{8,16,32,64, 128\}\). For each of these five values, run a Gibbs sampler of at least 5,000 iterations and obtain \(E[\theta_B - \theta_A| y_1....y_{N_A}, y_1....y_{N_B}]\). Describe the effects of the prior distribution for \(\gamma\) on the results.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.3
ya <- scan("./menchild30bach.dat")
yb <- scan("./menchild30nobach.dat")
yasum <- sum(ya); na <- length(ya)
ybsum <- sum(yb); nb <- length(yb)
atheta <- 2; btheta <- 1

set.seed(8675309)

S <- 10000
gamparam <- c(8, 16, 32, 64, 128)
ev <- c()


sim <- vector(mode = "list",length = 5)
for(i in 1:5){
  
  #starting values
  
  PHI <- matrix(nrow = S, ncol = 5)
  colnames(PHI) <- c("theta", "gamma", "theta_a", "theta_b", "b-a")
  phi <- c(mean(ya), mean(yb)/mean(ya), mean(ya), mean(yb), mean(yb)-mean(ya))
  PHI[1,] <- phi
  
  for(s in 2:S){
    
    #theta conditional
    phi[1] <- rgamma(1, atheta + yasum + ybsum, btheta + na + nb*phi[2])
    
    #gamma conditional
    phi[2] <- rgamma(1, gamparam[i] + ybsum, gamparam[i] + nb * phi[1])
    
    #theta_A
    phi[3] <- phi[1]
    
    #theta_B
    phi[4] <- phi[1]*phi[2]
    
    #theta_B - theta_A
    phi[5] <- phi[4] - phi[3]
    
    PHI[s,] <- phi
    
  }
  sim[[i]] <- PHI
  ev <- c(ev, mean(PHI[, 5]))
}

g = ggplot(data.frame(gamma_parameter = gamparam, theta_difference = ev), aes(x = gamma_parameter, y = theta_difference)) + geom_point() + geom_line() + theme_classic() + scale_x_continuous(breaks = seq(0, 140, 20))
g

As shown by the graph above,\(E[\theta_B - \theta_A| y_1....y_{N_A}, y_1....y_{N_B}]\) decreases as the parameters for \(\gamma\) get larger. That is, as our prior belief that the ratio \(\theta_B / \theta_A = 1\) becomes stronger, the expected difference between the two theta values decreases.