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 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 an efficient way of reducing a parameter vector with large dimension to smaller ones (Ex. Vectors with single parameter). Each iteration will result in each subvector randomly sampled using the subvector's posterior density and also dependent on the the current values of the other subvectors(Duchateau & Janssen, 2007, p.234).

## Limitation: If the parameter vector is very large and the subdivision will result in very small groups, Gibbs sampling will take a very long time to converge.

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

## Hamiltonian Monte Carlo cannot handle discrete parameters because there is no slope. We can only use continuous parameters.

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

## n_eff is the number of independent samples.
## Samples here means the number of iterations of the Markov chains, not data points. Markov chains are typically autocorrelated, making sequential samples not independent. 

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

## Rhat should approach 1 from the top.

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

##             mean          sd        5.5%       94.5%
## a[1]   0.8865640 0.015674842  0.86151253  0.91161538
## a[2]   1.0505706 0.009936056  1.03469086  1.06645034
## b[1]   0.1325115 0.074200548  0.01392474  0.25109835
## b[2]  -0.1425748 0.054746448 -0.23007023 -0.05507943
## sigma  0.1094880 0.005934470  0.10000357  0.11897243

dat_slim <- list(
    log_gdp_std = dd$log_gdp_std,
    rugged_std = dd$rugged_std,
    cid = as.integer( dd$cid )
)
str(dat_slim)

## List of 3
##  $ log_gdp_std: num [1:170] 0.88 0.965 1.166 1.104 0.915 ...
##  $ rugged_std : num [1:170] 0.138 0.553 0.124 0.125 0.433 ...
##  $ cid        : int [1:170] 1 2 2 2 2 2 2 2 2 1 ...

m9.1 <- 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 )
show( m9.1 )

## Hamiltonian Monte Carlo approximation
## 2000 samples from 4 chains
## 
## Sampling durations (seconds):
##         warmup sample total
## chain:1   0.05   0.04  0.09
## chain:2   0.04   0.03  0.07
## chain:3   0.04   0.03  0.07
## chain:4   0.06   0.03  0.09
## 
## Formula:
## 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)

traceplot( m9.1 )
y <- c(-1,1)
set.seed(11)
m9.2 <- ulam(
    alist(
        y ~ dnorm( mu , sigma ) ,
         mu <- alpha ,
    alpha ~ dnorm( 0 , 1000 ) ,
    sigma ~ dexp( 0.0001 )
) , data=list(y=y) , chains=3 )

## 
## SAMPLING FOR MODEL 'd26c527083e7eda89b17a8c2eccd6019' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 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.055 seconds (Warm-up)
## Chain 1:                0.033 seconds (Sampling)
## Chain 1:                0.088 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'd26c527083e7eda89b17a8c2eccd6019' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 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.068 seconds (Warm-up)
## Chain 2:                0.013 seconds (Sampling)
## Chain 2:                0.081 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'd26c527083e7eda89b17a8c2eccd6019' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 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.078 seconds (Warm-up)
## Chain 3:                0.016 seconds (Sampling)
## Chain 3:                0.094 seconds (Total)
## Chain 3:

## Warning: There were 42 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## Warning: Examine the pairs() plot to diagnose sampling problems

## Warning: The largest R-hat is 1.11, 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

show( m9.2 )

## Hamiltonian Monte Carlo approximation
## 1500 samples from 3 chains
## 
## Sampling durations (seconds):
##         warmup sample total
## chain:1   0.06   0.03  0.09
## chain:2   0.07   0.01  0.08
## chain:3   0.08   0.02  0.09
## 
## Formula:
## y ~ dnorm(mu, sigma)
## mu <- alpha
## alpha ~ dnorm(0, 1000)
## sigma ~ dexp(1e-04)

traceplot( m9.2 )

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

trankplot( m9.1 )

## Warning in if (class(x) == "numeric") x <- array(x, dim = c(length(x), 1)): the
## condition has length > 1 and only the first element will be used

## Warning in if (class(x) == "numeric") x <- array(x, dim = c(length(x), 1)): the
## condition has length > 1 and only the first element will be used

## Warning in if (class(x) == "numeric") x <- array(x, dim = c(length(x), 1)): the
## condition has length > 1 and only the first element will be used

## Warning in if (class(x) == "numeric") x <- array(x, dim = c(length(x), 1)): the
## condition has length > 1 and only the first element will be used

## Warning in if (class(x) == "numeric") x <- array(x, dim = c(length(x), 1)): the
## condition has length > 1 and only the first element will be used

trankplot( m9.2 )

## Warning in if (class(x) == "numeric") x <- array(x, dim = c(length(x), 1)): the
## condition has length > 1 and only the first element will be used

## Warning in if (class(x) == "numeric") x <- array(x, dim = c(length(x), 1)): the
## condition has length > 1 and only the first element will be used

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?

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

##             mean          sd        5.5%       94.5%
## a[1]   0.8865660 0.015675078  0.86151419  0.91161779
## a[2]   1.0505679 0.009936208  1.03468791  1.06644787
## b[1]   0.1325350 0.074201585  0.01394649  0.25112342
## b[2]  -0.1425568 0.054747270 -0.23005354 -0.05506012
## sigma  0.1094897 0.005934696  0.10000487  0.11897445

pairs(m8.3)

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

##             mean          sd        5.5%       94.5%
## a[1]   0.8865646 0.015680645  0.86150390  0.91162530
## a[2]   1.0505685 0.009939796  1.03468276  1.06645419
## b[1]   0.1325028 0.074227013  0.01387368  0.25113189
## b[2]  -0.1425733 0.054766564 -0.23010089 -0.05504579
## sigma  0.1095296 0.005940112  0.10003617  0.11902306

pairs(m8.3_unif)

## From the results, we can see it does not have detectible influence on the posterior distribution of sigma.

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?

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

##             mean          sd        5.5%       94.5%
## a[1]   0.8865649 0.015679232  0.86150647  0.91162335
## a[2]   1.0505696 0.009938884  1.03468537  1.06645388
## b[1]   0.1325026 0.074220566  0.01388385  0.25112145
## b[2]  -0.1425740 0.054761661 -0.23009371 -0.05505429
## sigma  0.1095195 0.005938737  0.10002822  0.11901072

pairs(m8.3_exp)

## There is not differences in the posterior distribution. 

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?

m9.1 <- 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 )

## recompiling to avoid crashing R session

precis( m9.1 , 2 )

##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8867544 0.015997741  0.86215748  0.91227352 3450.006 0.9989706
## a[2]   1.0505851 0.010190587  1.03446373  1.06672717 2534.798 1.0000621
## b[1]   0.1330407 0.075296395  0.01564941  0.25222800 2595.367 0.9988701
## b[2]  -0.1413581 0.056017108 -0.22998829 -0.04893445 2559.966 0.9997827
## sigma  0.1115539 0.005802838  0.10263806  0.12113434 2527.737 0.9998120

pairs( m9.1 )

## Warmup of 500 is enough in this case.

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 )

## 
## SAMPLING FOR MODEL 'bcf56ee89f6cf2a4224a4139ff01c7d4' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 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.014 seconds (Warm-up)
## Chain 1:                0.024 seconds (Sampling)
## Chain 1:                0.038 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

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? ## From the plot we can see plot a should be a normal distribution as the prior is aroung 0 and spread in between 2 and -2. Plot b is Cauchy distribution which contains some extreme value go up to over 30 and -40.