stan_model_string <- "
transformed data {
real sophie_height;
real gwendoline_height;
real sophie_y;
real gwendoline_y;
sophie_height <- 1.75; #m
gwendoline_height <- 1.90; #m
sophie_y <- 39; #pixels above kit
gwendoline_y <- 146; #pixels above kit
# The pixel values were read from this pixture:
# https://i.kinja-img.com/gawker-media/image/upload/1298715223159689320.jpg
# using WebPlotDigitizer: https://arohatgi.info/WebPlotDigitizer/app/
}
parameters {
real<lower=0> kit_height;
# catch all measurement error for all types of unconsidered circumstances, e.g.,
# Kit might be crouching in the picture, the perspective of the picture skewing the lengths, etc.
real<lower=0> height_measurement_error;
# The heights of the heels of the shoes on the actors
real<lower=0> sophie_heel_height;
real<lower=0> gwendoline_heel_height;
real<lower=0> kit_heel_height;
# The size of a pixel, in the image, in meters, bounded between reasonable values.
real<lower=0.0005, upper=0.05> pixel_in_m;
}
model {
# The difference in height between Kit and the two other actors.
real sophie_kit_diff;
real gwen_kit_diff;
# The prior on Kit's height, as I couldn't find any info on the mean and SD of UK
# male height, I used the Irish data from here: https://econ.upf.edu/docs/papers/downloads/1002.pdf
kit_height ~ normal(1.77, 0.08);
# Priors for the height of the heels of the actors, pure guesswork on my part.
sophie_heel_height ~ gamma(7.5, 150); # roughly 5 cm +- 2
gwendoline_heel_height ~ gamma(7.5, 150); # roughly 5 cm +- 2
kit_heel_height ~ gamma(2, 150); # roughly 1 cm +- 1
# Prior on the measurement error, roughly an error SD of 1 cm +- 1 .
height_measurement_error ~ gamma(4, 300);
sophie_kit_diff <- sophie_height + sophie_heel_height - (kit_height + kit_heel_height);
gwen_kit_diff <- gwendoline_height + gwendoline_heel_height - (kit_height + kit_heel_height);
sophie_y ~ normal( sophie_kit_diff / pixel_in_m, height_measurement_error / pixel_in_m);
gwendoline_y ~ normal( gwen_kit_diff / pixel_in_m, height_measurement_error / pixel_in_m);
}"
library(rstan)
fit <- stan(model_code = stan_model_string, iter = 2000)
##
## TRANSLATING MODEL 'stan_model_string' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'stan_model_string' NOW.
##
## SAMPLING FOR MODEL 'stan_model_string' NOW (CHAIN 1).
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## ...
plot(fit)
fit
## Inference for Stan model: stan_model_string.
## 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%
## kit_height 1.73 0.00 0.03 1.67 1.71 1.73
## height_measurement_error 0.01 0.00 0.01 0.00 0.01 0.01
## sophie_heel_height 0.05 0.00 0.02 0.02 0.04 0.05
## gwendoline_heel_height 0.05 0.00 0.02 0.02 0.04 0.05
## kit_heel_height 0.01 0.00 0.01 0.00 0.01 0.01
## pixel_in_m 0.00 0.00 0.00 0.00 0.00 0.00
## lp__ -105.60 0.09 1.95 -110.24 -106.59 -105.24
## 75% 97.5% n_eff Rhat
## kit_height 1.75 1.79 420 1.01
## height_measurement_error 0.02 0.03 602 1.00
## sophie_heel_height 0.06 0.09 417 1.00
## gwendoline_heel_height 0.06 0.10 1042 1.01
## kit_heel_height 0.02 0.03 970 1.00
## pixel_in_m 0.00 0.00 843 1.01
## lp__ -104.19 -102.89 447 1.01
##
## Samples were drawn using NUTS(diag_e) at Tue Jun 16 13:43:59 2015.
## 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).
post_samples <- as.data.frame(fit)
library(BEST)
plotPost(post_samples$kit_height, xlab = "Posterior height of Kit (in m)")
