Assignment #8

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
# The Metropolis algorithm works whenever the probability of proposing a jump to B from A is equal to the probability of proposing A from B, when the proposal distribution is 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?

# It explores posterior distribution more efficiently, so there are fewer rejections compared to base Metropolis algorithm.
# Gibbs sampling leverages information about the analytic form of likelihood and conjugate priors, so it can merge proposal and reject/accept steps.

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

# Hamiltonian Monte Carlo method couldn't handle discrete parameters. Method is based on the physical metaphor of the moving particle. This particle moves in the space of parameters with speed proportional to the gradient of likelihood, so the method should be able to differentiate the space. Discrete parameters aren't differentiable.

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 estimate how many independent samples we should get to collect same information as presented by the current sequence of posterior samples.

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

# Rhat should be close to 1. 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 (does'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?

#Should be uncorrelated, ie wildly going up and down, within reasonable estimates of the parameters. Something like:

plot(y=rnorm(1000,0,1),x=1:1000, type="l")

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

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?

library(rethinking)
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 )

m8.5 <- 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( m8.5 , depth=2 )

##             mean          sd        5.5%       94.5%
## a[1]   0.8865638 0.015675370  0.86151156  0.91161609
## a[2]   1.0505656 0.009936400  1.03468536  1.06644593
## b[1]   0.1325055 0.074202959  0.01391482  0.25109614
## b[2]  -0.1425764 0.054748275 -0.23007475 -0.05507811
## sigma  0.1094918 0.005934983  0.10000653  0.11897703

m8.5_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(m8.5_unif , depth=2)

##             mean          sd        5.5%       94.5%
## a[1]   0.8865702 0.015680391  0.86150987  0.91163046
## a[2]   1.0505661 0.009939628  1.03468062  1.06645151
## b[1]   0.1325022 0.074225809  0.01387504  0.25112939
## b[2]  -0.1425681 0.054765653 -0.23009423 -0.05504204
## sigma  0.1095277 0.005939856  0.10003468  0.11902075

pairs(m8.5)

pairs(m8.5_unif)

There is no visual difference between resulting models, at least difference is indistinguishable from several runs of the same model.

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?

m9.2 <- 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( m9.2 , depth=2 )

##             mean          sd       5.5%       94.5%
## a[1]   0.8865639 0.015678740  0.8615063  0.91162157
## a[2]   1.0505604 0.009938575  1.0346767  1.06644419
## b[1]   0.1325042 0.074218328  0.0138890  0.25111944
## b[2]  -0.1425747 0.054759960 -0.2300917 -0.05505775
## sigma  0.1095160 0.005938260  0.1000255  0.11900644

pairs(m9.2)

# There is no visual difference between resulting models.

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?

formula8m3 <- 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( m9.2 , depth=2 )

##             mean          sd       5.5%       94.5%
## a[1]   0.8865639 0.015678740  0.8615063  0.91162157
## a[2]   1.0505604 0.009938575  1.0346767  1.06644419
## b[1]   0.1325042 0.074218328  0.0138890  0.25111944
## b[2]  -0.1425747 0.054759960 -0.2300917 -0.05505775
## sigma  0.1095160 0.005938260  0.1000255  0.11900644

#too slow to run but results are calculated as below
# m8m3.w1 <- map2stan(formula8m3, data=dd, iter=1001, warmup=1)
# precis(m8m3.w1)#bad results, n_eff=1
# 
# m8m3.w10 <- map2stan(formula8m3, data=dd, iter=1010, warmup=10)
# precis(m8m3.w10)#not so bad results, n_eff=73, troubles with estimates sigma
# plot(m8m3.w10)
# 
# m8m3.w100 <- map2stan(formula8m3, data=dd.trim, iter=1100, warmup=100)
# precis(m8m3.w100)#enough of warmup
# plot(m8m3.w100)

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 /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -D_REENTRANT  -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -include stan/math/prim/mat/fun/Eigen.hpp   -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk -I/usr/local/include  -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/3.6/Resources/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 /Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/3.6/Resources/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.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 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.045428 seconds (Warm-up)
## Chain 1:                0.055077 seconds (Sampling)
## Chain 1:                0.100505 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

# There is no model for y in the formula, so data doesn't influence inference.
# Model will be sampling from the specified prioir distributions.
# We expect to see normal distribution for a and cauchy for b.
plot(mp)

precis(mp)

##        mean        sd      5.5%    94.5%     n_eff    Rhat4
## a 0.1089549  1.031883 -1.485812 1.709142 102.04481 1.021509
## b 3.4479565 32.515644 -6.115737 8.610613  55.51935 1.024127

samples <- extract.samples(mp)
hist(samples$a)#perfect Gausssian

hist(samples$b)#hard to examine, cause has sevral extreme values (as expected from cauchy thick tails)

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?