Notes and homework problems for Richard McElreath's Statistical Rethinking.
cavaunpeu / statistical-rethinking Goto Github PK
View Code? Open in Web Editor NEWSolutions for the practice problems
Solutions for the practice problems
Notes and homework problems for Richard McElreath's Statistical Rethinking.
should be: composite.posterior <- composite.unstandardized.posterior / sum(composite.unstandardized.posterior)
statistical-rethinking/chapter-2/homework.R
Line 191 in c8e39e4
2H4 - shouldnt the probabily for correctly identifying a species B panda is 0.65 equate to .65, instead of (1-.65) ?
line 191, replace
composite.posterior <- composite.unstandardized.posterior / sum(composite.posterior)
with
composite.posterior <- composite.unstandardized.posterior / sum(composite.unstandardized.posterior)
I think you have run the models incorrectly by failing to account for damage_norm and min_pressure as independent variables that appear outside of just being interaction terms. Otherwise you hide their influence outside of just being interactions.
Solution to 12M4
# 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
Your code is incomplete. Here is mine if it helps. To do fixed effects from 60 districts is a pain in the ass, that can be more easily done with the greta
package. Below it is the multilevel form alongside the plot with the district-level proportions after logistic
on the posterior a + alpha[district]
.
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
Hi,
Add 8H5 and 8H6 as per below, if you find appropriate:
8H5
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))
8H6
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))
Cheers,
F
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.