Status

Estimated and checked against book:

Stan code printed in the book or in the rethinking package:

This model is not discussed in my copy of the book:

Libraries

library(tidyverse)
library(tidybayes)
library(rstan)
library(patchwork)
library(rethinking)
options(mc.cores = 4)

Section 16.1: Geometric people

stan_data <- read.csv('data/Howell1.csv', sep = ';') %>%
             mutate(weight = weight / mean(weight),
                    height = height / mean(height)) %>%
             compose_data

stan_program <- "
data {
    int n;
    vector[n] weight;
    vector[n] height;
}
parameters {
    real p;
    real<lower=0> k;
    real<lower=0> sigma;
}
model {
    vector[n] mu;
    for (i in 1:n) {
        mu[i] = log(3.141593 * k * p^2 * height[i]^3);
    }
    weight ~ lognormal(mu, sigma);
    p ~ beta(2, 18);
    k ~ exponential(.5);
    sigma ~ exponential(1);
}
"
m16.1 <- stan(model_code = stan_program, data = stan_data)
m16.1
## Inference for Stan model: f6fb434c6f0d3e19436d6871028d8335.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## p       0.25    0.00 0.06   0.16   0.21   0.24   0.27   0.38   895    1
## k       5.82    0.08 2.68   2.10   3.99   5.24   7.10  12.49  1015    1
## sigma   0.21    0.00 0.01   0.19   0.20   0.21   0.21   0.22  1435    1
## lp__  576.80    0.04 1.22 573.61 576.25 577.10 577.69 578.20  1151    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Aug  2 09:50:54 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Section 16.3: Ordinary differential nut cracking

stan_data <- read.csv('data/Panda_nuts.csv', sep = ';') %>%
             mutate(age = age / max(age)) %>%
             compose_data

stan_program <- "
data {
    int n;
    int nuts_opened[n];
    vector[n] age; 
    vector[n] seconds; 
}
parameters {
    real phi;
    real k;
    real theta;
}
model {
    vector[n] lambda;
    for (i in 1:n) {
        lambda[i] = seconds[i] * phi * (1 - exp(-k * age[i]))^theta;
    }
    nuts_opened ~ poisson(lambda);
    phi ~ lognormal(log(1), .1);
    k ~ lognormal(log(2), .25);
    theta ~ lognormal(log(5), .25);
}
"
m16.4 <- stan(model_code = stan_program, data = stan_data)
m16.4
## Inference for Stan model: b0c1380fa4a6cfc0903d33822ca7fce6.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## phi      0.87    0.00 0.04    0.79    0.84    0.87    0.89    0.95  1456 1.00
## k        5.95    0.02 0.55    4.86    5.58    5.95    6.34    7.02  1001 1.01
## theta    9.78    0.06 1.98    6.41    8.37    9.59   11.04   14.12   966 1.01
## lp__  1993.91    0.03 1.18 1990.82 1993.36 1994.21 1994.78 1995.29  1209 1.01
## 
## Samples were drawn using NUTS(diag_e) at Sun Aug  2 09:51:15 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).