Chapter9.Rmd

---
title: "Chapter 9"
author: "Julin N Maloof"
date: "9/26/2019"
output: 
  html_document: 
    keep_md: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, error = TRUE)
```

## Video Notes

* Bayesian is about the posterior but doesn't care how you get there.  It is not about MCMC, that is just one tool to get the posterior.
* King Markov
- Metropolois Archipelago
- Must visit each island in proportion to population density.
- Flip a coin to choose proposed island on left or right. "proposal"
- census population of proposal island and current island
- move to proposal with probability of prop  pop / current pop
- repeat.  this ensures visiting each island in proportion to its population in the long run.
- why would you do this?  useful if you don't know the distribution of population sizes.  Or in this case the distribution of posterior probabilities.  Allow sampling from unknown posterior distribution.
* Metropolis algorithm
- will converge in the long run
- as long as proposals are symmetric
- not very efficient
- Useful to draw samples from posterior distribution
- Island: parameter values
- Population size: proportional to posterior probability
- works for any numbers of dimensions (parameters); continuous or discrete
- Markov chain: history doesn't matter, probability only depends on where you are.
* Why MCMC?
- can't write integrated posterior, or can't use it
- multilevel models, networks, phylogenies, spatial models are are hard to get integrated
- optimization (e.g. quap) not a good strategy in high dimensions -- must have full distributions
- MCMC is not fancy.  It is old and essential.
* Many MCMC strategies
- Metropolis-Hastings (MH): More general
- Gibbs sampling (GS): Efficient version of MH
- Metropolis and Givvs are "guess and check" strategies.  so quality proposals are essential.  If making dumb proposals then you don't move, and don't visit potentially important parts of the distribution.
- Hamiltonian Monte Carlo (HMC) fundamentally different, does not guess and check.
* 

## Problems

### 8E1

8E1. Which of the following is a requirement of the simple Metropolis algorithm?
(1) The parameters must bed iscrete.
(2) The likelihood function must be Gaussian.
(3) The proposal distribution must be symmetric.

3

### 8E2

By using conjugate priors (Whatever those are), allowing more efficient proposals.

### 8E3
_Which sort of parameters can Hamiltonian Monte Carlo not handle? Can you explain why?_

Cannot handle discrete parameters because it glides across the surface.  (Need a continuous surface).

### 8M1
_Re-estimate the terrain ruggedness model from the chapter, but now using a uniform prior and an exponential prior for the standard deviation, sigma. The uniform prior should be dunif(0,10) and the exponential should be dexp(1). Do the different priors have any detectible influence on the posterior distribution?_

get the data and transform it
```{r}
## R code 9.9
library(rethinking)
library(tidyverse)
options(mc.cores = parallel::detectCores())
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 )

## R code 9.11
dat_slim <- list(
  log_gdp_std = dd$log_gdp_std,
  rugged_std = dd$rugged_std,
  cid = as.integer( dd$cid )
)
str(dat_slim)
```

from chapter
```{r}
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 , iter=1000 )
```

```{r}
precis(m9.1, depth = 2)
```

```{r}
pairs(m9.1)
traceplot(m9.1)
trankplot(m9.1)
```

```{r}
m9.1.unif <- 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 ~ dunif( 0,10 )
  ) ,
  data=dat_slim , chains=4 , cores=4 , iter=1000 )
```

```{r}
m9.1.cauchy <- 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 ~ dcauchy( 0, 1 )
  ) ,
  data=dat_slim , chains=4 , cores=4 , iter=1000 )
```

```{r}
precis(m9.1, depth = 2)
```


```{r}
precis(m9.1.unif, depth = 2)
```

```{r}
precis(m9.1.cauchy, depth = 2)
```

```{r}
pairs(m9.1.unif)
```

```{r}
sigmapl <- tibble(exp = extract.samples(m9.1)$sigma,
                  unif = extract.samples(m9.1.unif)$sigma,
                  cauchy = extract.samples(m9.1.cauchy)$sigma) %>%
  gather(key=model, value=sigma)
sigmapl %>% 
  ggplot(aes(x=sigma, fill=model)) +
  geom_density(alpha=.4)
```


No big difference.  Maybe sigma posterior from exp distribution is a bit more peaked

### 8M2

_The Cauchy and exponential priors from the terrain ruggedness model are very weak. They can be made more informative by reducing their scale. Compare the dcauchy and dexp priors for progressively smaller values of the scaling parameter. As these priors become stronger, how does each influence the posterior distribution?_

```{r}
m9.1.exp.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 , iter=1000 )

m9.1.exp.01 <- 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( .01 )
  ) ,
  data=dat_slim , chains=4 , cores=4 , iter=1000 )
```

```{r}
precis(m9.1, depth=2)
```

```{r}
precis(m9.1.exp.1, depth=2)
``` 

```{r}
precis(m9.1.exp.01, depth=2)
```

```{r}
traceplot(m9.1.exp.1)
traceplot(m9.1.exp.01)
```

```{r}
sigmapl <- tibble(exp1 = extract.samples(m9.1)$sigma,
                  exp.1 = extract.samples(m9.1.exp.1)$sigma,
                  exp.01 = extract.samples(m9.1.exp.01)$sigma) %>%
  gather(key=model, value=sigma)
sigmapl %>% 
  ggplot(aes(x=sigma, fill=model)) +
  geom_density(alpha=.4)
```


```{r}
m9.1.cauchy.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 ~ dcauchy(0, .1 )
  ) ,
  data=dat_slim , chains=4 , cores=4 , iter=1000 )

m9.1.cauchy.01 <- 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 ~ dcauchy(0, .01 )
  ) ,
  data=dat_slim , chains=4 , cores=4 , iter=1000 )
```


```{r}
sigmapl <- tibble(cauchy = extract.samples(m9.1.cauchy)$sigma,
                  cauchy.1 = extract.samples(m9.1.cauchy.1)$sigma,
                  cauchy.01 = extract.samples(m9.1.cauchy.01)$sigma) %>%
  gather(key=model, value=sigma)
sigmapl %>% 
  ggplot(aes(x=sigma, fill=model)) +
  geom_density(alpha=.4)
```


```{r}
precis(m9.1.cauchy, depth=2)
```

```{r}
precis(m9.1.cauchy.1, depth=2)
``` 

```{r}
precis(m9.1.cauchy.01, depth=2)
```

No change

### 8H1

_Run the model below and then inspect the posterior distribution and explain what it is accom-
plishing_
```{r}
mp <- ulam(
  alist(
    a ~ dnorm(0,1),
    b ~ dcauchy(0,1)
  ),
  data=list(y=1),
  start=list(a=0,b=0),
  iter=1e4, warmup=100 , chains=4 )
```
```{r}
traceplot(mp)
trankplot(mp)
precis(mp)
```

```{r}
extract.samples(mp)$a %>% dens()
```

```{r}
extract.samples(mp)$b %>% dens()
```


The model simply samples from the normal and cauchy distributions

### 8H2

_Recall the divorce rate example from Chapter 5. Repeat that analysis, using map2stan this time, fitting models m5.1, m5.2, and m5.3. Use compare to compare the models on the basis of WAIC. Explain the results._

```{r}
# load data and copy
library(rethinking)
data(WaffleDivorce)
d <- WaffleDivorce

# standardize variables
d$A <- scale( d$MedianAgeMarriage )
d$D <- scale( d$Divorce )

## R code 5.2
sd( d$MedianAgeMarriage )

## R code 5.3
m5.1 <- quap(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bA * A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

```

```{r}
d$M <- scale( d$Marriage )
m5.2 <- quap(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM * M ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )
```


```{r}
m5.3 <- quap(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )
precis( m5.3 )
```

```{r}
compare(m5.1,m5.2,m5.3)
```


Now with stan

```{r}
dsmall <- list(D=d$D, A=d$A, M=d$M)
```


```{r}
m5.1.stan <- ulam(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bA * A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = dsmall,
  chains=4,
  cores=4,
  log_lik = TRUE)
```

```{r}
m5.2.stan <- ulam(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM * M ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = dsmall ,
  chains=4,
  cores=4,
  log_lik = TRUE)
```

```{r}
m5.3.stan <- ulam(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = dsmall ,
  chains=4,
  cores=4,
  log_lik = TRUE)
```

```{r}
precis(m5.1.stan)
precis(m5.2.stan)
precis(m5.3.stan)
```


```{r}
compare(m5.1.stan, m5.2.stan, m5.3.stan)
```

### optional: 8H6

_Modify the Metropolis algorithm code from the chapter to write your own simple MCMC estimator for globe tossing data and model from Chapter 2._

```{r}
iter <- 1e5
samples <- rep(0,iter)
current <- 0.5
for ( i in 1:iter ) {
  
  # record current position
  samples[i] <- current
  current_lik <- dbinom( 6 , size=9 , prob=current )
  
  # generate proposal
  proposal <- runif(1, min=0,max=1 )
  
  prop_lik <- dbinom(6, size=9, prob = proposal)
  
  # move?
  prob_move <- prop_lik/current_lik
  current <- ifelse( runif(1) < prob_move , proposal , current )
}
```

```{r}
dens(samples, adj=1)
```

plot the chain
```{r}
plot(x=1:iter,y=samples, type="l")
```


## third set of problems

### 8E4
_Explain the difference between the effective number of samples, n_eff as calculated by Stan,
and the actual number of samples._

n_eff takes autocorellation into account.

### 8M3

_Re-estimate one of the Stan models from the chapter, but at different numbers of warmup it- erations. Be sure to use the same number of sampling iterations in each case. Compare the n_eff values. How much warmup is enough?_

Using the marriage model

```{r}
m5.3.stan <- ulam(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = dsmall ,
  chains=4,
  cores=4,
  log_lik = TRUE)
precis(m5.3.stan)
```

0.12 seconds

```{r}
m5.3.stan.w100 <- ulam(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = dsmall ,
  chains=4,
  cores=4,
  iter=600,
  warmup=100,
  log_lik = TRUE)
precis(m5.3.stan.w100)
```

100 seems non-optimal as the sampling is less even.
Total: 0.059

```{r}
m5.3.stan.w10 <- ulam(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = dsmall ,
  chains=4,
  cores=4,
  iter=600,
  warmup=10,
  log_lik = TRUE)
precis(m5.3.stan.w10)
```

But 10 warmup samples is not that much worse.
Total time: 0.05 seconds

```{r}
m5.3.stan.w250 <- ulam(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = dsmall ,
  chains=4,
  cores=4,
  iter=600,
  warmup=250,
  log_lik = TRUE)
precis(m5.3.stan.w250)
```

### 8H3

_Sometimes changing a prior for one parameter has unanticipated effects on other parameters. This is because when a parameter is highly correlated with another parameter in the posterior, the prior influences both parameters. Here’s an example to work and think through. Go back to the leg length example in Chapter 5. Here is the code again, which simulates height and leg lengths for 100 imagined individuals:_

```{r}
## R code 9.27
N <- 100                          # number of individuals
height <- rnorm(N,10,2)           # sim total height of each
leg_prop <- runif(N,0.4,0.5)      # leg as proportion of height
leg_left <- leg_prop*height +     # sim left leg as proportion + error
  rnorm( N , 0 , 0.02 )
leg_right <- leg_prop*height +    # sim right leg as proportion + error
  rnorm( N , 0 , 0.02 )
# combine into data frame
d <- data.frame(height,leg_left,leg_right)
```


```{r}
## R code 9.28
m5.8s <- map2stan(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + bl*leg_left + br*leg_right ,
    a ~ dnorm( 10 , 100 ) ,
    bl ~ dnorm( 2 , 10 ) ,
    br ~ dnorm( 2 , 10 ) ,
    sigma ~ dcauchy( 0 , 1 )
  ) ,
  data=d, chains=4,
  log_lik = TRUE,
  start=list(a=10,bl=0,br=0,sigma=1) )
```

```{r}
trankplot(m5.8s)
plotchains(m5.8s)
precis(m5.8s)
pairs(m5.8s)
```


_Compare the posterior distribution produced by the code above to the posterior distribution pro- duced when you change the prior for br so that it is strictly positive:_


```{r}
## R code 9.29
m5.8s2 <- map2stan(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + bl*leg_left + br*leg_right ,
    a ~ dnorm( 10 , 100 ) ,
    bl ~ dnorm( 2 , 10 ) ,
    br ~ dnorm( 2 , 10 ) & T[0,] ,
    sigma ~ dcauchy( 0 , 1 )
  ) ,
  data=d, chains=4,
  log_lik = TRUE,
  start=list(a=10,bl=0,br=0,sigma=1) )

```


_Note that T[0,] on the right-hand side of the prior for br. What the T[0,] does is truncate the normal distribution so that it has positive probability only above zero. In other words, that prior ensures that the posterior distribution for br will have no probability mass below zero. Compare the two posterior distributions for m5.8s and m5.8s2. What has changed in the pos- terior distribution of both beta parameters? Can you explain the change induced by the change in_

```{r}
trankplot(m5.8s2)
plotchains(m5.8s2)
precis(m5.8s2)
pairs(m5.8s2)
```

So bR and bL are still inversely correlated, as expected, but now bR is always postive (of course) and bL is generally negative.  Interestingly the n_eff is lower.  Probably because in the previous modely that could flip between + and -

### 8H4

_For the two models fit in the previous problem, use DIC or WAIC to compare the effective numbers of parameters for each model. Which model has more effective parameters? Why?_

```{r}
compare(m5.8s, m5.8s2)
```

m5.8s has more effective parameters.  But why?

## Book

```{r, eval=FALSE, echo=FALSE}
## R code 9.1
num_weeks <- 1e5
positions <- rep(0,num_weeks)
current <- 10
for ( i in 1:num_weeks ) {
  # record current position
  positions[i] <- current
  
  # flip coin to generate proposal
  proposal <- current + sample( c(-1,1) , size=1 )
  # now make sure he loops around the archipelago
  if ( proposal < 1 ) proposal <- 10
  if ( proposal > 10 ) proposal <- 1
  
  # move?
  prob_move <- proposal/current
  current <- ifelse( runif(1) < prob_move , proposal , current )
}
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.2
D <- 10
T <- 1e3
Y <- rmvnorm(T,rep(0,D),diag(D))
rad_dist <- function( Y ) sqrt( sum(Y^2) )
Rd <- sapply( 1:T , function(i) rad_dist( Y[i,] ) )
dens( Rd )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.3
# U needs to return neg-log-probability
myU4 <- function( q , a=0 , b=1 , k=0 , d=1 ) {
  muy <- q[1]
  mux <- q[2]
  U <- sum( dnorm(y,muy,1,log=TRUE) ) + sum( dnorm(x,mux,1,log=TRUE) ) +
    dnorm(muy,a,b,log=TRUE) + dnorm(mux,k,d,log=TRUE)
  return( -U )
}
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.4
# gradient function
# need vector of partial derivatives of U with respect to vector q
myU_grad4 <- function( q , a=0 , b=1 , k=0 , d=1 ) {
  muy <- q[1]
  mux <- q[2]
  G1 <- sum( y - muy ) + (a - muy)/b^2 #dU/dmuy
  G2 <- sum( x - mux ) + (k - mux)/d^2 #dU/dmuy
  return( c( -G1 , -G2 ) ) # negative bc energy is neg-log-prob
}
# test data
set.seed(7)
y <- rnorm(50)
x <- rnorm(50)
x <- as.numeric(scale(x))
y <- as.numeric(scale(y))
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.5
library(shape) # for fancy arrows
Q <- list()
Q$q <- c(-0.1,0.2)
pr <- 0.3
plot( NULL , ylab="muy" , xlab="mux" , xlim=c(-pr,pr) , ylim=c(-pr,pr) )
step <- 0.03
L <- 11 # 0.03/28 for U-turns --- 11 for working example
n_samples <- 4
path_col <- col.alpha("black",0.5)
points( Q$q[1] , Q$q[2] , pch=4 , col="black" )
for ( i in 1:n_samples ) {
  Q <- HMC2( myU4 , myU_grad4 , step , L , Q$q )
  if ( n_samples < 10 ) {
    for ( j in 1:L ) {
      K0 <- sum(Q$ptraj[j,]^2)/2 # kinetic energy
      lines( Q$traj[j:(j+1),1] , Q$traj[j:(j+1),2] , col=path_col , lwd=1+2*K0 )
    }
    points( Q$traj[1:L+1,] , pch=16 , col="white" , cex=0.35 )
    Arrows( Q$traj[L,1] , Q$traj[L,2] , Q$traj[L+1,1] , Q$traj[L+1,2] ,
            arr.length=0.35 , arr.adj = 0.7 )
    text( Q$traj[L+1,1] , Q$traj[L+1,2] , i , cex=0.8 , pos=4 , offset=0.4 )
  }
  points( Q$traj[L+1,1] , Q$traj[L+1,2] , pch=ifelse( Q$accept==1 , 16 , 1 ) ,
          col=ifelse( abs(Q$dH)>0.1 , "red" , "black" ) )
}
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.6
HMC2 <- function (U, grad_U, epsilon, L, current_q) {
  q = current_q
  p = rnorm(length(q),0,1) # random flick - p is momentum.
  current_p = p
  # Make a half step for momentum at the beginning
  p = p - epsilon * grad_U(q) / 2
  # initialize bookkeeping - saves trajectory
  qtraj <- matrix(NA,nrow=L+1,ncol=length(q))
  ptraj <- qtraj
  qtraj[1,] <- current_q
  ptraj[1,] <- p
  ```
  
  
  ```{r, eval=FALSE, echo=FALSE}
  ## R code 9.7
  # Alternate full steps for position and momentum
  for ( i in 1:L ) {
    q = q + epsilon * p # Full step for the position
    # Make a full step for the momentum, except at end of trajectory
    if ( i!=L ) {
      p = p - epsilon * grad_U(q)
      ptraj[i+1,] <- p
    }
    qtraj[i+1,] <- q
  }
  ```
  
  
  ```{r, eval=FALSE, echo=FALSE}
  ## R code 9.8
  # Make a half step for momentum at the end
  p = p - epsilon * grad_U(q) / 2
  ptraj[L+1,] <- p
  # Negate momentum at end of trajectory to make the proposal symmetric
  p = -p
  # Evaluate potential and kinetic energies at start and end of trajectory
  current_U = U(current_q)
  current_K = sum(current_p^2) / 2
  proposed_U = U(q)
  proposed_K = sum(p^2) / 2
  # Accept or reject the state at end of trajectory, returning either
  # the position at the end of the trajectory or the initial position
  accept <- 0
  if (runif(1) < exp(current_U-proposed_U+current_K-proposed_K)) {
    new_q <- q  # accept
    accept <- 1
  } else new_q <- current_q  # reject
  return(list( q=new_q, traj=qtraj, ptraj=ptraj, accept=accept ))
}
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.9
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 )

## R code 9.10
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 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.11
dat_slim <- list(
  log_gdp_std = dd$log_gdp_std,
  rugged_std = dd$rugged_std,
  cid = as.integer( dd$cid )
)
str(dat_slim)
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.12
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=1 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.13
precis( m9.1 , depth=2 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.14
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 , iter=1000 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.15
show( m9.1 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.16
precis( m9.1 , 2 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.17
pairs( m9.1 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.18
traceplot( m9.1 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.19
trankplot( m9.1 , n_cols=2 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.20
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=2 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.21
precis( m9.2 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.22
set.seed(11)
m9.3 <- ulam(
  alist(
    y ~ dnorm( mu , sigma ) ,
    mu <- alpha ,
    alpha ~ dnorm( 1 , 10 ) ,
    sigma ~ dexp( 1 )
  ) ,
  data=list(y=y) , chains=2 )
precis( m9.3 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.23
set.seed(41)
y <- rnorm( 100 , mean=0 , sd=1 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.24
m9.4 <- ulam(
  alist(
    y ~ dnorm( mu , sigma ) ,
    mu <- a1 + a2 ,
    a1 ~ dnorm( 0 , 1000 ),
    a2 ~ dnorm( 0 , 1000 ),
    sigma ~ dexp( 1 )
  ) ,
  data=list(y=y) , chains=2 )
precis( m9.4 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.25
m9.5 <- ulam(
  alist(
    y ~ dnorm( mu , sigma ) ,
    mu <- a1 + a2 ,
    a1 ~ dnorm( 0 , 10 ),
    a2 ~ dnorm( 0 , 10 ),
    sigma ~ dexp( 1 )
  ) ,
  data=list(y=y) , chains=2 )
precis( m9.5 )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.26
mp <- map2stan(
  alist(
    a ~ dnorm(0,1),
    b ~ dcauchy(0,1)
  ),
  data=list(y=1),
  start=list(a=0,b=0),
  iter=1e4, warmup=100 , WAIC=FALSE )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.27
N <- 100                          # number of individuals
height <- rnorm(N,10,2)           # sim total height of each
leg_prop <- runif(N,0.4,0.5)      # leg as proportion of height
leg_left <- leg_prop*height +     # sim left leg as proportion + error
  rnorm( N , 0 , 0.02 )
leg_right <- leg_prop*height +    # sim right leg as proportion + error
  rnorm( N , 0 , 0.02 )
# combine into data frame
d <- data.frame(height,leg_left,leg_right)
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.28
m5.8s <- map2stan(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + bl*leg_left + br*leg_right ,
    a ~ dnorm( 10 , 100 ) ,
    bl ~ dnorm( 2 , 10 ) ,
    br ~ dnorm( 2 , 10 ) ,
    sigma ~ dcauchy( 0 , 1 )
  ) ,
  data=d, chains=4,
  start=list(a=10,bl=0,br=0,sigma=1) )
```


```{r, eval=FALSE, echo=FALSE}
## R code 9.29
m5.8s2 <- map2stan(
  alist(
    height ~ dnorm( mu , sigma ) ,
    mu <- a + bl*leg_left + br*leg_right ,
    a ~ dnorm( 10 , 100 ) ,
    bl ~ dnorm( 2 , 10 ) ,
    br ~ dnorm( 2 , 10 ) & T[0,] ,
    sigma ~ dcauchy( 0 , 1 )
  ) ,
  data=d, chains=4,
  start=list(a=10,bl=0,br=0,sigma=1) )

```