Maximum likelihood estimation of mixture parameter between two binomial distributions

Question

I am struggling with the following problem and I hope I can get some hints.

I would like to model the result of sampling balls from 2 urns, each one with a somehow infinite number of black and white balls present at a known rate. Let $p_a, p_b$ determine the fraction of black balls in urn $A, B$ respectively. $p_a$ and $p_b$ are known.

Now let's imagine that a certain number of balls ($n$, known) is drawn from the two urns. This is done many times and constitutes the data ($y$, known).

I would like to estimate the fraction of balls coming from $A$ and $B$, ($\theta$ and $(1-\theta)$), assuming that this fraction stays the same across sampling and given that $p_a, p_b$ are very different.

For instance, if $A$ has 5% of black balls and $B$ has 95% of black balls, I would expect that observing many times something around 10 black balls over 100 drawn would lead to a low $\theta$ estimation.

I think I can write:

$$ y = Bin(n, p) \\ y = Bin(n_a, p_a) + Bin(n_b,p_b) $$ With: $$ n_a + n_b = n \\ n_a = Bin(n, \theta)\\ n_b = n - n_a \\ \theta\in [0,1] $$

So, I can try to write the likelihood, but I am not sure this is correct:

$$ P(y|\theta) = (Bin_a(n_a) + Bin_b(n-n_a)) \times Bin(n,\theta) $$

And I am stuck here. I am familiar with MLE estimation of common distributions, but this goes beyond my formal education and I don't understand whether the problem is ill-posed, or how it should be tackled. Maybe $n_a$ should be marginalized through a summation? Would that be feasible in practice? I hope I stated the problem clearly.

Thank you for your time

Answer 1

Your approach of maximum likelihood estimation is to express the likelihood $L(y; \theta)$ of $Y$ as a function of $\theta$ and then compute the $\theta$ that is the/an argmax of the likelihood.

The distribution of $Y$ is that of a sum of two independent Binomial random variables, which, in the general case of $p_a\ne p_b$ that you are considering, is a bit unwieldy. There actually exists an R package for this situation, sinib. In the belonging paper, you can find references to possible approximations, but the authors state that "Analytical solutions to the density and distribution are usually cumbersome to find and difficult to compute."

Thus, provided that you don't find an analytic approach in the references of this paper, you could use the sinib package to program a function that assigns to each $\theta$ the belonging likelihood $L(y;\theta)$ of your data, and then feed this function to an optimizer like optim.