# y: Number of infants that were breast feeding at the time of release from hospital
# n: Total number of infants
# y ~ Binomial(n,p)
y <- c( 2, 2, 7, 7,16,14)
n <- c( 6, 5, 9, 9,20,15)
# x: Gestational age of the infant (in weeks) at the time of birth
# logit(p) = log(p/(1-p)) = Beta0 + Beta1 x
x <- c(28,29,30,31,32,33)
make_likelihood <- function(x,y,n,b0,s0,b1,s1) {
function(state) {
p <- 1 / ( 1 + exp(-(state[1]+state[2]*x)) )
sum(dbinom(y,n,p,log=TRUE))
}
}
# b0: Mean of the normal prior distribution for Beta0
# b1: Mean of the normal prior distribution for Beta1
# s0: Std. dev. of the normal prior distribution for Beta0
# s1: Std. dev. of the normal prior distribution for Beta1