The Stan project for statistical computation has a great collection of curated case studies which anybody can contribute to, maybe even me, I was thinking. But I don’t have time to worry about that right now because I’m on vacation, being on the yearly visit to my old family home in the north of Sweden. What I do worry about is that my son will be stung by a bumblebee. His name is Torsten, he’s almost two years old, and he loves running around barefoot on the big lawn. Which has its fair share of bumblebees. Maybe I should put shoes on him so he wont step on one, but what are the chances, really.
Well, what are the chances? I guess if I only had Data on the bumblebee density of the lawn.
Data on the size of Torsten’s feet and how many steps he takes when running around.
A reasonable Bayesian model, maybe implemented in Stan.
I could figure that out. “How hard can it be?”, I thought. And so I made an attempt.
Getting the data
To get some data on bumblebee density I marked out a 1 m² square on a representative part of the lawn. During the course of the day, now and then, I counted up how many bumblebees sat in the square.
Most of the time I saw zero bumblebees, but 1 m² is not that large of an area. Let’s put the data into R:
bumblebees parameters {
real mu_bees;
real mu_steps;
real precision_steps;
}
model {
# Since we have contrained the parameters to be positive we get implicit
# half-cauchy distributions even if we declare them to be 'full'-cauchy.
mu_bees ~ cauchy(0, (4.0 + 4.0) / (10.0 * 10.0) );
mu_steps ~ cauchy(0, 60);
precision_steps ~ cauchy(0, 1);
bumblebees ~ poisson(mu_bees);
toddler_steps ~ neg_binomial_2(mu_steps, precision_steps);
}
The final step is to predict how many bumblebees Torsten will step on during, say, one hour of active play. We do this in the generated quantities code block in Stan. The code below will step through 60 minutes (a.k.a. one hour) and for each minute: (1) Sample pred_steps from the negative binomial, (2) calculate the area (m²) covered by these steps, and (3) sample the number of bees in this area from the Poisson giving a predicted stings_by_minute. Finally we sum these 60 minutes worth of strings into stings_by_hour.
generated quantities {
int stings_by_minute[60];
int stings_by_hour;
int pred_steps;
real stepped_area;
for(minute in 1:60) {
pred_steps = neg_binomial_2_rng(mu_steps, precision_steps);
stepped_area = pred_steps * foot_m2;
stings_by_minute[minute] = poisson_rng(mu_bees * stepped_area);
}
stings_by_hour = sum(stings_by_minute);
}
Looking at the result
We have data, we have a model, and now all we need to do is fit it. After we’ve put the whole Stan model into bee_model_code as a string we do it like this using the rstan interface:
# First we put all the data into a list
data_list