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