Beta-Binomial parameter estimation The MLE or method of moments estimation of parameters of a beta-binomial distribution makes use of (c, y) -- total number and positive counts. However, if we only have one such pair, then $\frac{\alpha}{\beta}$, of the distribution, is supervised but $\alpha+\beta$ tends to $\infty$. Now if I break the observation: (c, y) into several randomly subsetted counts $(c_i, y_i)$, where $c_i=\sum_{k\in S}1, y_i=\sum_{k\in S}\rho$, S is randomly chosen subset of examples and $\rho \in \{0, 1\}$ is observation per point. We still get the same estimate for the mean: $\frac{\alpha}{\alpha+\beta}$, but this time $\alpha+\beta$ would also be bounded.
Is there something wrong with the argument? Which of these estimates is correct? In general, what qualifies as observations for a beta-binomial and what could be wrong with such random subsetting of observation?
Edit
Problem setting:
I have several {0, 1} observations from only one Bernoulli arm. I would like to estimate the mean of this arm while keeping track of its uncertainty, so modelling the accuracy of the arm through beta distribution is a natural choice. Since I observe a certain number of positive and total count (y, c) for this arm, the probability of the observations would be $\theta \sim Beta(\alpha, \beta)$ and $(c, y) \sim Binomial(\theta)$.
 A: 
Is there something wrong with the argument? Which of these estimates is correct? In general, what qualifies as observations for a beta-binomial and what could be wrong with such random subsetting of observation?

Which is correct depends on how you choose to model the data! There are two distinct models that are masquerading as one when it comes to topics like these.

*

*You have a sequence of Bernoulli trials (same probability of success, independent) and you're interested in the total. If you model the probability of success as a scalar (one true $p$) then the estimation of $p$ does not depend on the way you partition the trials- any way you do so, you end up with a likelihood of a binomial distribution with all counts grouped together. If, on the other hand, you place a prior distribution on the probability of success, say, $p \sim \text{Beta}(\alpha, \beta)$, then the posterior distribution is  $\text{Beta}(\alpha + k, \beta + n - k)$, and your posterior mean is as you gave above. Here, the beta distribution models your uncertainty about the true parameter, but it does not vary between trials.


*If instead you have a sequence of Beta-Binomial trials, where some number of binary outcomes are observed for each, then the generative process for the data is different- for each trial, a success probability is sampled from a beta distribution, and all the binary outcomes within that trial share the same success probability. Crucially, the beta distribution here reflects the fact that the success probability is itself a random variable. More concretely, if you have 10 Bernoulli trials and group them all together, you would be generating one success probability from a beta distribution, then using that to sample an outcome from 10 trials via a binomial distribution. If you instead split the Bernoulli trials into two groups of 5 each, you would be generating two success probabilities and drawing from a binomial over 5 trials twice, once with each.
It's sneaky because in either case, the posterior predictive distribution is a Beta-Binomial, but the posterior distribution in 1. is a beta distribution over the success probability, while in 2. the prior distribution is an ugly mess with no known closed form for the normalization constant over the alpha and beta parameters.
The reason that you're seeing $\alpha + \beta$ blow up when you use maximum likelihood on a single partition is straightforward- if you had a single binomial trial, the MLE is the proportion of successes. When you use a beta binomial, you're in essence averaging the binomial likelihood over a distribution on p. If p is any less or more than the proportion of successes, the likelihood is smaller than at the MLE. So, the best beta-binomial distribution is one that is always exactly the proportion of successes. However, that cannot be the case for any finite values of the beta parameters. With the ratio of the two parameters fixed but letting their sum go to infinity, the beta distribution approaches a point mass at the mean.
What does this mean for your problem? It depends on the specifics.

*

*do you have any compelling reason to partition the trials- are some under differing conditions or is there any reason to believe that the success probability varies between partitions? If so, partitioning them in a way that the context supports and using a beta-binomial is the way to go.

*do you have reason to believe that the success probability is fixed across all the trials grouped together and you want to find a posterior distribution on the probability of success encoding the uncertainty? Then you should use a beta prior distribution and update the parameters with success/failure counts appropriately, but you should not try to optimize over the beta distribution parameters with a single group of trials. Either they are prior parameters and fixed, or you're going to have the parameters go to infinity indicating that a binomial distribution is sufficient and you can skip the beta distribution altogether.

This is really only an issue when optimizing over the beta parameters on a set of trials with identical numbers of successes and failures per trial. If you're doing this estimation for data where most groups have multiple binomial trials with potentially varying success probabilities within each group but some only have one trial, there are a few ways to handle it. The simplest is to build a hierarchical model, where you have multiple beta-binomial trials sharing parameters per group, and each group has its parameters drawn from a common prior. More precisely, a common way to do this is
$$ m_i \sim Beta(a', b') $$
$$ s_i \sim Gamma(c, d) $$
$$ \alpha_i = s_i m_i $$
$$ \beta_i = s_i (1-m_i) $$
Where $i$ ranges over the groups, and I've reparameterized the per-group beta distributions in terms of a "precision" $s = \alpha + \beta$ and a mean parameter $m$. This way, the gamma prior on $s$ prevents the convergence to a point mass and the mean probability of success $m$ is decoupled from the precision. This is no more complex than directly optimizing $\alpha$ and $\beta$ numerically, and it has the benefit of giving a straightforward way to  handle cases with a single trial without resorting to using special cases. You'll need to use an MCMC method to get an idea of the posterior on the $s$ and $ m$ parameters as this model is not completely conjugate but you'd need to do that even if you used a conjugate prior on the beta parameters since the normalization constant is intractable.
