How can we determine the parameters of a Beta-Binomial Distribution for given mean and variance? Given a beta-binomial random variate $X$ with $N$ known, how can I choose $\alpha$, $\beta$ such that the distribution's mean matches a chosen quantity $\mu$ and its variance matches a chosen quantity $\sigma^2$?
 A: Start with the mean and variance given on Wikipedia:
$\mu = \frac{n \alpha}{\alpha + \beta}$
$\sigma^2 = \frac{n\alpha\beta(\alpha + \beta + n)}{(\alpha + \beta)^2(\alpha + \beta + 1)}$
We can solve the first equation for $\beta$ and plug into the second to get the answer:
$\beta = \alpha \frac{n-\mu}{\mu}$
$\alpha = \frac{n(\mu(n-\mu) - \sigma^2)}{(1 + \frac{n-\mu}{\mu})(n\sigma^2 - \mu (n - \mu)}$
Here is an R implementation (the mean and var funcs are given below):
# Generate some data
alpha <- 2000
beta <- 100 
N <- 20
mu <- bb.mean(alpha, beta, N)
sig <- sqrt(bb.var(alpha, beta, N)) 

#Estimate the params
alpha <- N * (mu* (N - mu) - sig^2) / ((1 + (N-mu) / mu) * 
       (sig^2*N - mu * (N - mu)))
beta <- alpha * (N - mu) / mu


#Ensure accuracy
> alpha
[1] 200
> beta
[1] 100

Now for the variance and mean functions I used:
#' Get the mean of the beta-binom dist for fixed params
bb.mean <- function(alpha, beta, n) {
    alpha * n / (alpha + beta)
}

#' Get the variance of the beta dist for fixed params
bb.var <- function(alpha, beta, n) {
    num <- n * alpha * beta * (alpha + beta + n)
    denom <- (alpha + beta)^2 * (alpha + beta + 1)
    return(num / denom)
}

