# A Bayesian model of the height of Kit Harington using the info given in this article:
# http://jezebel.com/serious-question-how-tall-is-kit-harrington-1711475702
# by Rasmus BÃ¥Ã¥th, @rabaath, http://www.sumsar.net
# Disclaimer: I think there is a 50% chance I made a horrible mistake somewhere...

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:
# using WebPlotDigitizer: http://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
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)")

# Comment
#
# From this analysis it is unclear how tall Kit is, there is much uncertainty
# in the posterior distribution, but according to the analysis
# (which might be quite off) there's a 50% probability he's between
# 1.71 and 1.75 m tall. It is stated in the article that he is NOT 5'8'' (173 cm),
# but according to this analysis it's not an unreasonable height, as the mean of
# the posterior is 173 cm.