cavaunpeu / statistical-rethinking Goto Github PK

composite.posterior <- composite.unstandardized.posterior / sum(composite.posterior)

composite.posterior <- composite.unstandardized.posterior / sum(composite.unstandardized.posterior)

# Sun Feb 24 17:47:16 2019 ------------------------------
library(rethinking)
data("chimpanzees")
head(chimpanzees)
chimpanzees$recipient <- NULL
chimpanzees$block_id <- chimpanzees$block

# New cross-clasified model
mod1 <- map2stan(alist(
      pulled_left ~ dbinom(1, p),
      logit(p) <- a_actor[actor] + a_block[block_id] +
            (bP + bPC * condition)*prosoc_left,
      a_actor[actor] ~ dnorm(mu_actor, sigma_actor),
      a_block[block_id] ~ dnorm(mu_block, sigma_block),
      c(bP, bPC, mu_actor, mu_block) ~ dnorm(0, 10),
      c(sigma_actor, sigma_block) ~ dcauchy(0, 1)),
      data = chimpanzees, warmup = 1000, iter = 5000, chains = 4, cores = 3)

# Chapter's model for reference (incl. a and var intercept priors with mean == 0)
mod2 <- map2stan(alist(
      pulled_left ~ dbinom(1, p),
      logit(p) <- a + a_actor[actor] + a_block[block_id] +
            (bP + bPC * condition)*prosoc_left,
      a_actor[actor] ~ dnorm(0, sigma_actor),
      a_block[block_id] ~ dnorm(0, sigma_block),
      c(a, bP, bPC) ~ dnorm(0, 10),
      c(sigma_actor, sigma_block) ~ dcauchy(0, 1)),
      data = chimpanzees, warmup = 1000, iter = 5000, chains = 4, cores = 3)

compare(mod1, mod2)
precis(mod1)
precis(mod2)
# WAIC clearly favours the chapter model
# the posterior dists in the new model have much larger sds (~6),
# and this is reflected in n_eff from the HMC samples

data("bangladesh")
bangladesh$district <- as.integer(factor(bangladesh$district))

mod <- map2stan(alist(
      use.contraception ~ dbinom(1, p),
      logit(p) <- a + alpha[district],
      alpha[district] ~ dnorm(0, sigma),
      a ~ dnorm(0, 10),
      sigma ~ dcauchy(0, 1)),
      data = bangladesh)
samp <- extract.samples(mod)
samp$alpha <- logistic(apply(samp$alpha, 2, function(k){k + samp$a}))
plot(1:60, colMeans(samp$alpha), ylim = c(0, 1),
     xlab = "District ID", ylab = "Expected proportion contraception")
cis <- apply(samp$alpha, 2, HPDI)
segments(x0 = 1:60, x1 = 1:60, y0 = cis[1,], y1 = cis[2,])
# There is high variability in the intervals

num_weeks <- 1e5
positions <- rep(0,num_weeks)
current <- 10

# Make up populations for the ten islands
probs_islands <- c(4, 23, 41, 13, 7, 21, 33, 2, 31, 9)
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 <- probs_islands[proposal]/probs_islands[current]
      current <- ifelse( runif(1) < prob_move , proposal , current )
}
barplot(table(positions))

trials <- 1e5
positions <- rep(NA, trials)
current <- "W"
probs_waterland <- c("W" = .7, "L" = .3)

for ( i in 1:trials ) {
      # record current position
      positions[i] <- current
      
      # if current is L, proposal is W and vice versa
      proposal <- ifelse(current == "W", "L", "W")
      
      # move? prob_move will always either be .7/.3 or .3/.7
      prob_move <- probs_waterland[proposal]/probs_waterland[current]
      current <- ifelse(runif(1) < prob_move, proposal, current)
}
barplot(table(positions))