Code
# prior: Beta(alpha, beta)
= 1
alpha = 2
beta
# observations: binomial draws
= 3
n_successes = 3
n_trials
= function(p) {
gamma_beta_binomial if (p < 0 || p > 1) return(0.0)
dbeta(p, alpha, beta) * dbinom(x = n_successes, size = n_trials, prob = p)
}