# Status

Estimated and checked against book:

• m16.1
• m16.4

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

• m16.2
• m16.5

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

• m16.3

# 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).