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's distribution of proposed parameter values is adjusted to currect parameter values. Gibbs sampler uses pairs of priors and likelihoods that have anlytic solutions for the posterior of an individual parameter.

#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 will be no slope. It only uses continuous parameters.

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

# The effective number of samples is an estimate of the number of independent samples from the posterior distribution, in terms of estimating some function like the posterior mean. 

# n_eff is an estimate of effective number of samples, for the purpose of estimating the posterior mean.

# The actual number of samples are 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

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?

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.8865716 0.015674421  0.86152087  0.91162238
## a[2]   1.0505660 0.009935779  1.03468675  1.06644534
## b[1]   0.1325854 0.074198500  0.01400185  0.25116892
## b[2]  -0.1426479 0.054744865 -0.23014080 -0.05515506
## sigma  0.1094849 0.005934046  0.10000113  0.11896863
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 )

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

#traceplot( m9.2 )

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

#trankplot( m9.1 )
#trankplot( m9.2 )

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.8865636 0.015675243  0.86151154  0.91161567
## a[2]   1.0505666 0.009936317  1.03468640  1.06644671
## b[1]   0.1325056 0.074202382  0.01391591  0.25109538
## b[2]  -0.1425758 0.054747837 -0.23007337 -0.05507814
## sigma  0.1094909 0.005934860  0.10000582  0.11897593
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.8865658 0.015680682  0.86150507  0.91162659
## a[2]   1.0505696 0.009939818  1.03468388  1.06645538
## b[1]   0.1325027 0.074227169  0.01387335  0.25113206
## b[2]  -0.1425728 0.054766683 -0.23010055 -0.05504508
## sigma  0.1095299 0.005940146  0.10003636  0.11902336
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.8865654 0.015680696  0.86150466  0.9116262
## a[2]   1.0505694 0.009939827  1.03468363  1.0664552
## b[1]   0.1325017 0.074227239  0.01387223  0.2511312
## b[2]  -0.1425731 0.054766734 -0.23010093 -0.0550453
## sigma  0.1095300 0.005940160  0.10003644  0.1190235
pairs(m8.3_unif)

#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.8865646 0.015678785  0.86150691  0.91162237
## a[2]   1.0505696 0.009938597  1.03468580  1.06645339
## b[1]   0.1325028 0.074218530  0.01388722  0.25111831
## b[2]  -0.1425742 0.054760113 -0.23009139 -0.05505691
## sigma  0.1095163 0.005938302  0.10002572  0.11900682
pairs(m8.3_exp)

# There is no 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 )

#precis( m9.1 , 2 )

#pairs( m9.1 )

#500 warmup is enough.

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 )

#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 around 0 and spread between 2 and -2, while plot b is Cauchy distribution because it contains some extreme value that’s up to over 30 and -40.