This chapter has been an informal introduction to Markov chain Monte Carlo (MCMC) estimation. The goal has been to introduce the purpose and approach MCMC algorithms. The major algorithms introduced were the Metropolis, Gibbs sampling, and Hamiltonian Monte Carlo algorithms. Each has its advantages and disadvantages. The ulam function in the rethinking package was introduced. It uses the Stan (mc-stan.org) Hamiltonian Monte Carlo engine to fit models as they are defined in this book. General advice about diagnosing poor MCMC fits was introduced by the use of a couple of pathological examples.
Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them. Problems are labeled Easy (E), Medium (M), and Hard(H).
Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points.
9E1. Which of the following is a requirement of the simple Metropolis algorithm?
##3. The proposal distribution must be symmetric.9E2. Gibbs sampling is more efficient than the Metropolis algorithm. How does it achieve this extra efficiency? Are there any limitations to the Gibbs sampling strategy?
#Gibbs sampling is more efficient than the Metropolis algorithm since through the Gibbs sampling method you can get a good estimate of the posterior with many fewer samples than a comparable metropolis approach. The improvement arises from adaptive proposals in which the distribution of proposed parameter values adjusts itself intelligently depending upon the parameter values at the moment.
#Gibbs sampling leverages information about the analytic form of likelihood and conjugate priors, so it can merge proposal and reject/accept steps
#Some of the limitations in using Gibbs sampling strategy, is the Long convergence time especially with the high dimensional data.
#Convergence time also depends on the shape of the distribution. And Difficulty in finding the posterior for each variable.9E3. Which sort of parameters can Hamiltonian Monte Carlo not handle? Can you explain why?
# Hamiltonian Monte Carlo cannot handle discrete parameters.
## Because it cannot glide through discrete parameters without slopes.
#Considering the particles in space theory, since the energy of the total system needs to conserved as the model runs a physics simulation, if and when the total energy changes, the numerical approximation is bad, this might reject the proposal. 9E4. Explain the difference between the effective number of samples, n_eff as calculated by Stan, and the actual number of samples.
# N_effective aims to estimate the number of 'ideal' samples.
#Ideal samples are entirely uncorrelated. Due to way MCMC works each next sample is actually correlated with the previous one to some extent.
# n_eff as calculated by Stan: It is an estimated of effective number of samples, for the purpose of estimating the posterior mean.
# The actual number of samples: samples we use for accurate inference.9E5. Which value should Rhat approach, when a chain is sampling the posterior distribution correctly?
#Rhat should approach 1 when a chain is sampling the posterior distribution correctly.
#It reflects the fact that inner variance and outer variance between chains is roughly the same, so we could expect that inference is not broken (doesn't depend on the chain).9E6. Sketch a good trace plot for a Markov chain, one that is effectively sampling from the posterior distribution. What is good about its shape? Then sketch a trace plot for a malfunctioning Markov chain. What about its shape indicates malfunction?
# library(rethinking)
# data(rugged)
# df <- rugged
# df$log_gdp <- log(df$rgdppc_2000)
# dd <- df[ complete.cases(d$rgdppc_2000) , ]
# dd$log_gdp_std <- dd$log_gdp / mean(dd$log_gdp)
# dd$rugged_std <- dd$rugged / max(dd$rugged)
# dd$cid <- ifelse( dd$cont_africa == 1 , 1 , 2)
# m9e6 <- quap(
# alist(
# log_gdp_std ~ dnorm(mu, sigma),
# mu <- a[cid] + b[cid]*( rugged_std - 0.215 ),
# a[cid] ~ dnorm(1, 0.1 ),
# b[cid] ~ dnorm(0, 0.3 ),
# sigma ~ dexp(1)),
# data = dd)
# precis(m9e6 , depth = 2)
# dat <- list(
# log_gdp_std = dd$log_gdp_std,
# rugged_std = dd$rugged_std,
# cid = as.integer( dd$cid ))
# str(dat)
# HMC <- ulam(
# alist(
# log_gdp_std ~ dnorm( mu , sigma ) ,
# mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) ,
# a[cid] ~ dnorm( 1 , 0.1 ) ,
# b[cid] ~ dnorm( 0 , 0.3 ) ,
# sigma ~ dexp( 1 )
# ) , data=dat , chains=4 , cores=4 )
#show( HMC)
# traceplot( HMC )
#
# y <- c(-1,1)
# set.seed(11)
# HMC2 <- ulam(
# alist(
# y ~ dnorm( mu , sigma ) ,
# mu <- alpha ,
# alpha ~ dnorm( 0 , 1000 ) ,
# sigma ~ dexp( 0.0001 )
# ) , data=list(y=y) , chains=3 )
#show(HMC2)
#traceplot(HMC2)9E7. Repeat the problem above, but now for a trace rank plot.
#trankplot(HMC)
#trankplot(HMC2)9M1. Re-estimate the terrain ruggedness model from the chapter, but now using a uniform prior for the standard deviation, sigma. The uniform prior should be dunif(0,1). Visualize the priors. Use ulam to estimate the posterior. Visualize the posteriors for both models. Does the different prior have any detectible influence on the posterior distribution of sigma? Why or why not?
# data(rugged)
# df <- rugged
# df$log_gdp <- log(df$rgdppc_2000)
# dd <- df[ complete.cases(df$rgdppc_2000) , ]
#
# dd$log_gdp_std <- dd$log_gdp/ mean(dd$log_gdp)
# dd$rugged_std<- dd$rugged/max(dd$rugged)
#
# dd$cid<-ifelse(dd$cont_africa==1,1,2)
#
# HMC <- quap(
# alist(
# log_gdp_std ~ dnorm( mu , sigma ) ,
# mu <- a[cid] + b[cid]* (rugged_std-0.215) ,
# a[cid] ~ dnorm(1,0.1),
# b[cid] ~ dnorm(0,0.3),
# sigma ~ dexp(1)
# ) ,
# data=dd)
#
# precis(HMC , depth=2)
# pairs(HMC)
# HMC_unif <- quap(
# alist(
# log_gdp_std ~ dnorm( mu , sigma ) ,
# mu <- a[cid] + b[cid]* (rugged_std-0.215) ,
# a[cid] ~ dnorm(1,0.1),
# b[cid] ~ dnorm(0,0.3),
# sigma ~ dunif(0,1)
# ) ,
# data=dd)
#
# precis(HMC_unif , depth=2)
#pairs(HMC_unif)9M2. Modify the terrain ruggedness model again. This time, change the prior for b[cid] to dexp(0.3). What does this do to the posterior distribution? Can you explain it?
# HMC_exp <- quap(
# alist(
# log_gdp_std ~ dnorm( mu , sigma ) ,
# mu <- a[cid] + b[cid]* (rugged_std-0.215) ,
# a[cid] ~ dnorm(1,0.1),
# b[cid] ~ dnorm(0,0.3),
# sigma ~ dexp(0.3)
# ) ,
# data=dd)
#
# precis(HMC_exp , depth=2)
#pairs(HMC_exp)9M3. Re-estimate one of the Stan models from the chapter, but at different numbers of warmup iterations. Be sure to use the same number of sampling iterations in each case. Compare the n_eff values. How much warmup is enough?
# HMC_stan <- ulam(
# alist(
# log_gdp_std ~ dnorm( mu , sigma ) ,
# mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) ,
# a[cid] ~ dnorm( 1 , 0.1 ) ,
# b[cid] ~ dnorm( 0 , 0.3 ) ,
# sigma ~ dexp( 1 )
# ) , data=dat_slim , chains=4 , cores=4, warmup=5, iter=1000 )
#precis( HMC_stan , depth=2 )
#pairs(HMC_stan)
#traceplot(HMC_stan)9H1. Run the model below and then inspect the posterior distribution and explain what it is accomplishing.
# mp <- ulam(
# alist(
# a ~ dnorm(0,1),
# b ~ dcauchy(0,1)
# ), data=list(y=1) , chains=1 )
# precis( mp , 2 )
# pairs(mp)
# traceplot(mp)Compare the samples for the parameters a and b. Can you explain the different trace plots? If you are unfamiliar with the Cauchy distribution, you should look it up. The key feature to attend to is that it has no expected value. Can you connect this fact to the trace plot?