Stat. 481 Midterm

Instructions: This is an open-book midterm. You can ask questions of the instructor only. Compete the midterm questions in this R Notebook. Change the filename to use your lastname and firstname. The code needed to answer the questions is at the bottom of the notebook.

Question 1

Redo Problem 1, Exercise 6.1, from Homework 4.

1. Change to, “suppose we get a tail.”

Answers:

2. Change to, “flip again and get a tail.”

Answers:

3. Change to, “third time and get a head.”

Answers:

4. Change to, “the order H, T, T instead of T, T, H.”

Answers:

Question 2

1. What is the Bayesian model presented in the Rose Garden Event model? What data is used?

Answers:

2. State the likelihood for the data.

Answers:

3. What are the other parameters, et, theta, and nu?

Answers:

4. What prior are used for the parameters in the model?

Answers:

5. What are the posterior estimates and 95% HDIs for each parameter in the model?

Answers:

parameter | estimate | 95% HDI

eta | | theta | | nu | | Py | |

6. Describe what the Py values are in the JAGS code. Give the posterior estimate and 95% HDI for Py.

Answers:

Question 1 Code

Problem 1

Do Exercise 6.1

source("DBDA2E-utilities.R")  # Load definitions of graphics functions etc.
source("BernBeta.R")          # Load the definition of the BernBeta function

1. Start with a prior distribution that expresses some uncertainty that a coin is fair: beta( theta | 4,4). Flip the coin once; suppose we get a head. What is the posterior distribution?

# Specify the prior:
a = 4
b = 4

Prior = c(a,b)       # Specify Prior as vector with the two shape parameters.

# Specify the data:
N = 1                          # The total number of flips.
z = 1                          # The number of heads.
Data = c(rep(0,N-z),rep(1,z))  # Convert N and z into vector of 0's and 1's.

#openGraph(width=5,height=7)
posterior = BernBeta( priorBetaAB=Prior, Data=Data , plotType="Bars" , 
                      showCentTend="Mode" , showHDI=TRUE , showpD=FALSE )
#saveGraph(file="BernBetaExample",type="png")

2. Use the posterior from the previous flip as the prior for the next flip. Suppose we flip again and get a head. Now what is the new posterior?

# Specify the prior:
a = 5
b = 4

Prior = c(a,b)       # Specify Prior as vector with the two shape parameters.

# Specify the data:
N = 1                          # The total number of flips.
z = 1                          # The number of heads.
Data = c(rep(0,N-z),rep(1,z))  # Convert N and z into vector of 0's and 1's.

#openGraph(width=5,height=7)
posterior = BernBeta( priorBetaAB=Prior, Data=Data , plotType="Bars" , 
                      showCentTend="Mode" , showHDI=TRUE , showpD=FALSE )
#saveGraph(file="BernBetaExample",type="png")

3. Using that posterior as the prior for the next flip, flip a third time and get a tail. Now what is the new posterior? (Hint: Type post = BernBeta(posterior, c(0)).)

# Specify the prior:
a = 6
b = 4

Prior = c(a,b)       # Specify Prior as vector with the two shape parameters.

# Specify the data:
N = 1                          # The total number of flips.
z = 0                          # The number of heads.
Data = c(rep(0,N-z),rep(1,z))  # Convert N and z into vector of 0's and 1's.

#openGraph(width=5,height=7)
posterior = BernBeta( priorBetaAB=Prior, Data=Data , plotType="Bars" , 
                      showCentTend="Mode" , showHDI=TRUE , showpD=FALSE )
#saveGraph(file="BernBetaExample",type="png")

post = BernBeta(posterior, c(0))

4. Do the same three updates but in the order T, H, H instead of H, H, T. Is the final posterior distribution the same for both orderings of the flip results?

# Specify the prior:
a = 4
b = 4

Prior = c(a,b)       # Specify Prior as vector with the two shape parameters.

# Specify the data:
N = 3                          # The total number of flips.
z = 2                          # The number of heads.
Data = c(rep(0,N-z),rep(1,z))  # Convert N and z into vector of 0's and 1's.

#openGraph(width=5,height=7)
posterior = BernBeta( priorBetaAB=Prior, Data=Data , plotType="Bars" , 
                      showCentTend="Mode" , showHDI=TRUE , showpD=FALSE )
#saveGraph(file="BernBetaExample",type="png")

Question 2 Code

Rose Garden Event - Gibbs Sampler - Predictive Distribution

Simualted Data

require(rjags) 

source("DBDA2E-utilities.R") 

Now simulate the sum assuming a particular probability of a positive test.

prob_nu <- 0.05
n <- 300 

y <- rbinom(1, n, prob_nu)
y

JAGS need the data to be in a list of values, vectors, and matrices.

dataList <- list(
  y = y,
  Ntotal = n
)

dataList

Model

modelString <- "
model {
  y ~ dbin( nu, Ntotal )   # Likelihood in code
  nu = eta*pi + (1 - theta)*(1-pi)
  pi ~ dbeta( 1 , 1 )                # prior
  eta ~ dbeta( 910 , 90 )            # prior
  theta ~ dbeta( 950 , 50 )          # prior
  Py ~ dbin( nu, Ntotal )            # Predictive distribution
}
"

Inits

pi_init <- 0.004
eta_init <- 0.91
theta_init <- 0.95

initsList = list( pi=pi_init,
                  eta=eta_init,
                  theta=theta_init) 

Updates

jagsModel <- jags.model( file = textConnection(modelString), 
                         data = dataList, 
                         inits = initsList, 
                         n.chains = 5, 
                         n.adapt = 10000 )

update( jagsModel , n.iter=5000 )

codaSamples = coda.samples( jagsModel,
                            variable.names=c("nu","pi","eta","theta","Py"),
                            n.iter=5000 )

save( codaSamples , file=paste0("Mcmc.Rdata") )

Examine the chains:

Convergence diagnostics:

plot(codaSamples)

diagMCMC( codaObject=codaSamples, parName="nu" )
diagMCMC( codaObject=codaSamples, parName="pi" )
diagMCMC( codaObject=codaSamples, parName="eta" )
diagMCMC( codaObject=codaSamples, parName="theta" )

traceplot(codaSamples)

plotPost( codaSamples[,"nu"], main="nu", xlab=bquote(nu) )
plotPost( codaSamples[,"pi"], main="pi", xlab=bquote(pi) )
plotPost( codaSamples[,"eta"], main="eta", xlab=bquote(eta) )
plotPost( codaSamples[,"theta"], main="theta", xlab=bquote(theta) )

plotPost( codaSamples[,"nu"] , main="nu" , xlab=bquote(nu) , 
          cenTend="median"  )
plotPost( codaSamples[,"pi"] , main="pi" , xlab=bquote(pi) , 
          cenTend="median"  )
plotPost( codaSamples[,"eta"] , main="eta" , xlab=bquote(eta) , 
          cenTend="median"  )
plotPost( codaSamples[,"theta"] , main="theta" , xlab=bquote(theta) , 
          cenTend="median"  )

summary(codaSamples)

Examine the chains:

Convergence diagnostics:

diagMCMC( codaObject=codaSamples , parName="nu" )
diagMCMC( codaObject=codaSamples , parName="pi" )
diagMCMC( codaObject=codaSamples , parName="eta" )
diagMCMC( codaObject=codaSamples , parName="theta" )

Posterior descriptives:

plotPost( codaSamples[,"nu"] , main="nu" , xlab=bquote(nu) )
plotPost( codaSamples[,"pi"] , main="pi" , xlab=bquote(pi) )
plotPost( codaSamples[,"eta"] , main="eta" , xlab=bquote(eta) )
plotPost( codaSamples[,"theta"] , main="theta" , xlab=bquote(theta) )

Re-plot with different annotations:

plotPost( codaSamples[,"nu"] , main="nu" , xlab=bquote(nu) , 
          cenTend="median"  )
plotPost( codaSamples[,"pi"] , main="pi" , xlab=bquote(pi) , 
          cenTend="median"  )
plotPost( codaSamples[,"eta"] , main="eta" , xlab=bquote(eta) , 
          cenTend="median"  )
plotPost( codaSamples[,"theta"] , main="theta" , xlab=bquote(theta) , 
          cenTend="median"  )