```
# Set alpha and beta
alp <- 5
bet <- 7
# How many of these beta-bernoullis do we want to sum?
n <- 10
# How many trials do we want to do?
N <- 100000
# Generate the success probability from the beta
pi <- rbeta(n*N, alp, bet)
# Generate the observation from the binomial
x <- rbinom(n*N, 1, pi)
# Get the desired sum
y <- apply(matrix(x, ncol = n, byrow = T), 1, sum)
# We claim that obs will follow a binomial with n = n, p = alp/(alp + bet)
# This gets us the empirical estimate of the pmf
obs <- table(y)/length(y)
# This gets us the theoretical estimate of the pmf (but only the ones we got outcomes for)
exp <- dbinom(as.numeric(names(table(y))), n, alp/(alp + bet))
# Look at them together. If these are essentially the same that provides evidence
# that it just ends up being binomial
cbind(obs, exp)
```