This code shows how the posterior is updated with each new observation. The initial prior distribution is a beta-distribution which is rather flexible and has support on [0, 1]. The prior distribution will be updated based on the binomial model (which will actually result in a posterior dfistributino that is itself a beta distribution).
The function f()
will plot the distribution after each
observation and then update the posterior. It takes a vector of binary
outcomes as its first argument, and alpha
and
beta
are the parameters of the prior distribution.
f <- function(y, alpha, beta) {
x <- seq(0, 1, .01)
plot(x, dbeta(x, alpha, beta), type="l",
ylab="", ylim=c(0, 5))
readline(prompt = "Pause. Press <Enter> to continue...")
for (i in y) {
alpha <- alpha+i
beta <- beta + (1-i)
lines(x, dbeta(x, alpha, beta))
readline(prompt = "Pause. Press <Enter> to continue...")
}
}
Now we can generate some data and see how the prior is updated
# Generate data. N=20, probability of success = 1/3
y <- rbinom(20, size=1, p=1/3)
# Starting with a flat prior
f(y, alpha=1, beta=1)
# Starting with a prior that favours high probabilities
f(y, alpha=9, beta=1)
f(y, alpha=.5, beta=.5)
2024