hotbath

hotbath.ipynb

model.stan

data {
  int<lower=0> Y[4];
}

parameters {
  real mu;
  real<lower=0> sigma;
}

transformed parameters {
  vector<lower=0>[4] theta;
  theta[1] = normal_cdf(39.5 | mu, sigma);
  theta[2] = normal_cdf(40.5 | mu, sigma) - normal_cdf(39.5 | mu, sigma);
  theta[3] = normal_cdf(41.5 | mu, sigma) - normal_cdf(40.5 | mu, sigma);
  theta[4] = 1 - normal_cdf(41.5 | mu, sigma);
}

model {
  Y ~ multinomial(theta);
}

generated quantities {
  int<lower=0> Y_pred[4];
  Y_pred = multinomial_rng(theta, sum(Y));
}