# 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: 
  # http://i.kinja-img.com/gawker-media/image/upload/1298715223159689320.jpg
  # 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
  # male height, I used the Irish data from here: http://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)")

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