Chapter 9 - Markov Chain Monte Carlo

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. 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.

Questions

9E1. Which of the following is a requirement of the simple Metropolis algorithm?

  1. The parameters must be discrete.
  2. The likelihood function must be Gaussian.
  3. The proposal distribution must be symmetric.
#3

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's distribution of proposed parameter values is adjusted to currect parameter values.
#With Gibbs sampling, high dimension spaces are concentrated. When GS gets stuck, it degenerates towards random walk.

9E3. Which sort of parameters can Hamiltonian Monte Carlo not handle? Can you explain why?

# Hamiltonian Monte Carl cannot handle discrete parameters because there is no slope. It can only be used for continuous parameters

9E4. Explain the difference between the effective number of samples, n_eff as calculated by Stan, and the actual number of samples.

# Effective sample size, n_eff, is an estimate of the number of imdependent draws from the posterior distribution, but not all samples are accepted, n_eff is usually smaller than actual sample size.

9E5. Which value should Rhat approach, when a chain is sampling the posterior distribution correctly?

1 
## [1] 1

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?

X = rnorm(100, 0, 1)
aa = rnorm(1, 0, 1)
bb = rnorm(1, 0, 1)

print(aa, bb)
## [1] -0
simdata = as.double(aa* X + bb)

d = as.data.frame(x = list(y = simdata, X = X))

m_good <- ulam(
    alist(
        y ~ dnorm( mu, sigma),
        mu <- b +  a*X,
        b ~ dnorm( 0 , 1 ),
        a ~ dnorm( 0 , 1 ),
        sigma ~ dexp(1)
    ) , data=d, chains=4, cores=1 , sample=TRUE )
## Trying to compile a simple C file
## Running /usr/local/Cellar/r/4.0.2_1/lib/R/bin/R CMD SHLIB foo.c
## clang -I"/usr/local/Cellar/r/4.0.2_1/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/4.0/site-library/Rcpp/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/4.0/site-library/BH/include" -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/4.0/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include   -fPIC  -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:88:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '7b0ca6606d71a755d3d6fe33ac34aff5' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 4.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.46 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.147129 seconds (Warm-up)
## Chain 1:                0.215694 seconds (Sampling)
## Chain 1:                0.362823 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '7b0ca6606d71a755d3d6fe33ac34aff5' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.262101 seconds (Warm-up)
## Chain 2:                0.164315 seconds (Sampling)
## Chain 2:                0.426416 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '7b0ca6606d71a755d3d6fe33ac34aff5' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.187261 seconds (Warm-up)
## Chain 3:                0.601077 seconds (Sampling)
## Chain 3:                0.788338 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '7b0ca6606d71a755d3d6fe33ac34aff5' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.7e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.352724 seconds (Warm-up)
## Chain 4:                0.595369 seconds (Sampling)
## Chain 4:                0.948093 seconds (Total)
## Chain 4:
## Warning: There were 520 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 14 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 2.58, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
traceplot( m_good)
## [1] 1000
## [1] 1
## [1] 1000

d = as.data.frame(list(y = runif(100)))
m_bad <- ulam(
    alist(
        y ~ dnorm( mu, sigma),
        mu <- b +  a,
        b ~ dnorm( 0 , 1000 ),
        a ~ dnorm( 0 , 1000 ),
        sigma ~ dexp(1)
    ) , data=d, chains=4, cores=1 , sample=TRUE )
## Trying to compile a simple C file
## Running /usr/local/Cellar/r/4.0.2_1/lib/R/bin/R CMD SHLIB foo.c
## clang -I"/usr/local/Cellar/r/4.0.2_1/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/4.0/site-library/Rcpp/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/4.0/site-library/BH/include" -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/4.0/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include   -fPIC  -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:88:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '7d3d1dfeb74b5330df8ad6fcda0f1688' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.29984 seconds (Warm-up)
## Chain 1:                1.52896 seconds (Sampling)
## Chain 1:                2.8288 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '7d3d1dfeb74b5330df8ad6fcda0f1688' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.28931 seconds (Warm-up)
## Chain 2:                1.61719 seconds (Sampling)
## Chain 2:                2.9065 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '7d3d1dfeb74b5330df8ad6fcda0f1688' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.28019 seconds (Warm-up)
## Chain 3:                1.42483 seconds (Sampling)
## Chain 3:                2.70502 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '7d3d1dfeb74b5330df8ad6fcda0f1688' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1.29262 seconds (Warm-up)
## Chain 4:                1.56028 seconds (Sampling)
## Chain 4:                2.8529 seconds (Total)
## Chain 4:
## Warning: There were 1571 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 2.14, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
traceplot( m_bad)
## [1] 1000
## [1] 1
## [1] 1000

9E7. Repeat the problem above, but now for a trace rank plot.

trankplot(m_good)

trankplot(m_bad)

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). Use ulam to estimate the posterior. Does the different prior have any detectible influence on the posterior distribution of sigma? Why or why not?

No significant affect are detected, since the actual value would fall into both distribution.

data(rugged)
d <- rugged
d$log_gdp <- log(d$rgdppc_2000)
dd <- d[ 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)

dd = as.data.frame(list(log_gdp_std = dd$log_gdp_std,
                        rugged_std = dd$rugged_std,
                        cid = dd$cid ))

m9_org <- 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)


m9_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(m9_org , depth=2)
##             mean          sd       5.5%       94.5%
## a[1]   0.8865666 0.015674654  0.8615155  0.91161773
## a[2]   1.0505659 0.009935935  1.0346863  1.06644542
## b[1]   0.1324996 0.074199678  0.0139142  0.25108503
## b[2]  -0.1424674 0.054745890 -0.2299619 -0.05497289
## sigma  0.1094866 0.005934286  0.1000025  0.11897075
precis(m9_unif , depth=2)
##             mean          sd        5.5%       94.5%
## a[1]   0.8865605 0.015680587  0.86149992  0.91162113
## a[2]   1.0505681 0.009939764  1.03468245  1.06645377
## b[1]   0.1325123 0.074226772  0.01388358  0.25114101
## b[2]  -0.1425757 0.054766384 -0.23010295 -0.05504843
## sigma  0.1095292 0.005940061  0.10003588  0.11902261
m9_org %>% extract.samples %>% .$sigma %>% dens(.,xlab="sigma",col="blue")
mtext("posterior")

m9_unif %>% extract.samples %>% .$sigma %>% dens(.,xlab="sigma",col="blue")
mtext("posterior")

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?

After changing, b[cid] to dexp(0.3) we are limiting the effect of betas to be positive, forcing no nagtive effect. and also as beta are assume to be exponentially distribute, the posterior distribution might be having a bigger right tail as well.

m9.2 <- 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] ~ dexp(0.3) ,
        sigma ~ dexp( 1 )
    ) , data=dd , chains=4 , cores=4 )
## Trying to compile a simple C file
## Running /usr/local/Cellar/r/4.0.2_1/lib/R/bin/R CMD SHLIB foo.c
## clang -I"/usr/local/Cellar/r/4.0.2_1/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/4.0/site-library/Rcpp/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/4.0/site-library/BH/include" -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/4.0/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include   -fPIC  -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:88:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
m9_org %>% extract.samples %>% .$sigma %>% dens(.,xlab="sigma",col="blue")
mtext("posterior")

m9.2 %>% extract.samples %>% .$sigma %>% dens(.,xlab="sigma",col="blue")
mtext("posterior")

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?

A : within those values that I tried, there is no big difference.

m9.3_10 <- 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=dd , chains=4 , cores=4 , warmup = 10)
## Trying to compile a simple C file
## Running /usr/local/Cellar/r/4.0.2_1/lib/R/bin/R CMD SHLIB foo.c
## clang -I"/usr/local/Cellar/r/4.0.2_1/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/4.0/site-library/Rcpp/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/4.0/site-library/BH/include" -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/4.0/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include   -fPIC  -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:88:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Warning: There were 1 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
m9.3_100 <- 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=dd , chains=4 , cores=4 , warmup = 100)
## recompiling to avoid crashing R session
## Trying to compile a simple C file
## Running /usr/local/Cellar/r/4.0.2_1/lib/R/bin/R CMD SHLIB foo.c
## clang -I"/usr/local/Cellar/r/4.0.2_1/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/4.0/site-library/Rcpp/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/4.0/site-library/BH/include" -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/4.0/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include   -fPIC  -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:88:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
m9.3_500 <- 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=dd , chains=4 , cores=4 , warmup = 500)
## recompiling to avoid crashing R session
## Trying to compile a simple C file
## Running /usr/local/Cellar/r/4.0.2_1/lib/R/bin/R CMD SHLIB foo.c
## clang -I"/usr/local/Cellar/r/4.0.2_1/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/4.0/site-library/Rcpp/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/4.0/site-library/BH/include" -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/4.0/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include   -fPIC  -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:88:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
precis(m9.3_10, depth = 2)
##             mean         sd        5.5%       94.5%     n_eff    Rhat4
## a[1]   0.8871737 0.02256627  0.86149213  0.91213066 1045.9909 1.001695
## a[2]   1.0509648 0.01290246  1.03508882  1.06665675 1329.3739 1.001382
## b[1]   0.1298607 0.09138849  0.01096412  0.25614217  710.8040 1.002987
## b[2]  -0.1416455 0.06264292 -0.22958022 -0.05097411  985.4112 1.005176
## sigma  0.1126707 0.01868597  0.10244146  0.12258399  442.3459 1.006743
precis(m9.3_100, depth = 2)
##             mean          sd         5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8859330 0.015675124  0.861348185  0.91101727 4125.182 1.0004045
## a[2]   1.0506310 0.010012089  1.034931472  1.06660558 4399.748 0.9992967
## b[1]   0.1296442 0.077185478  0.007725311  0.25511973 1672.435 1.0016397
## b[2]  -0.1409475 0.054263219 -0.228436538 -0.05391883 2460.366 0.9990495
## sigma  0.1115207 0.006242702  0.101998980  0.12197295 2028.682 1.0020798
precis(m9.3_500, depth = 2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8868888 0.016364954  0.86109353  0.91353476 2735.157 0.9986328
## a[2]   1.0506821 0.010208107  1.03384353  1.06680103 3191.299 0.9991357
## b[1]   0.1328722 0.078134485  0.01152166  0.25716951 2189.742 0.9999174
## b[2]  -0.1410870 0.055968558 -0.22761495 -0.05026241 2212.314 0.9991064
## sigma  0.1117332 0.006409402  0.10225555  0.12214929 2606.243 0.9987764

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 )
## Trying to compile a simple C file
## Running /usr/local/Cellar/r/4.0.2_1/lib/R/bin/R CMD SHLIB foo.c
## clang -I"/usr/local/Cellar/r/4.0.2_1/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/4.0/site-library/Rcpp/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/4.0/site-library/BH/include" -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/4.0/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/4.0/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/4.0/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include   -fPIC  -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:88:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /usr/local/lib/R/4.0/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Dense:1:
## /usr/local/lib/R/4.0/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '3bd3f4d287e9cccab124308e5415245c' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.026396 seconds (Warm-up)
## Chain 1:                0.030415 seconds (Sampling)
## Chain 1:                0.056811 seconds (Total)
## Chain 1:
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess
pairs(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?

#Plot a has a normal distribution and plot b has a Cauchy distribution because where it has very long tails