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

  • 1
    $\begingroup$ The problem is not ill-posed but do you really need the MLE? It would be much easier to work with the method of moment estimate obtained by equating the sample mean to the theoretical mean given by $\theta n p_a + (1-\theta)n p_b$ and solving for $\theta$. $\endgroup$ Jul 8, 2022 at 10:13
  • $\begingroup$ @JarleTufto I think you are right, and your comment made me also realize that I made a mistake in formalizing the problem. Indeed in my case, the two urns are mixed from the beginning with a fraction $\theta$, and then the samples are taken. I believe I can then model the distribution as $Bin(b, p)$ with $p = \theta p_a + (1-\theta) p_b$. $\endgroup$
    – gianMa
    Jul 8, 2022 at 12:15

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.