data { int N; int num[N]; int dead[N]; vector[N] dose; } transformed data { vector[N] centered_dose; real mean_dose; mean_dose = mean(dose); centered_dose = dose - mean_dose; } parameters { real alpha_star; real beta; } transformed parameters { vector[N] m; m = alpha_star + beta * centered_dose; } model { alpha_star ~ normal(0.0, 1.0E4); beta ~ normal(0.0, 1.0E4); dead ~ binomial_logit(num, m); } generated quantities { real alpha; real LD50; real p[N]; // real llike[N]; // real rhat[N]; for (i in 1:N) { p[i] = inv_logit(m[i]); // llike[i] = dead[i]*log(p[i]) + (num[i]-dead[i])*log(1-p[i]); // rhat[i] = p[i]*num[i]; // fitted values } alpha = alpha_star - beta*mean_dose; LD50 = -1.0*alpha/beta; }