MLE of multinomial distribution with missing values

Question

I am trying to solve a problem and the results I get seem counter-intuitive.

We randomly throw $n$ balls into an area partitioned into 3 bins $b_1,b_2,b_3$. The size of each bins is proportional to the probability the ball will fall in it. Let's call these probabilities $p_1,p_2,p_3$. This can be described by a multinomial distribution.

Now let's say I throw 12 balls, and I know how many landed in each bin ($x_1=3,x_2=6,x_3=3$).

I would like to estimate the size of the bin from the observations. For this I use Maximum Likelihood. It can be shown that the MLE will be $p_1=3/12,p_2=6/12,p_3=3/12$. This is pretty intuitive.

It turns out that the actual likelihood at this point is:

$L(p_1=0.25,p_2=0.5,p_3=0.25|x_1=3,x_2=6,x_3=3)=$

$=\frac{12!}{3!6!3!}0.25^30.6^60.25^6=0.07050$

Now, let's assume I knew in advance that $p_1=p_3$. How would that change my result? It would not - I would still get the same parameter values $p_1=0.25,p_2=0.5,p_3=0.25$.

The twist comes now: let's assume I cannot observe balls that landed in $b_3$. If I know that 12 balls were thrown I am fine, since I can calculate $b_3=n-b_1-b_2=12-3-6=3$. but what happens if I don't know $n$?

I figure that in this case, I would need to estimate $x_3$ (or equivalently $n$) as well. However, if I use MLE, the results start looking weird. Intuitively, I would expect that if I observe $x_1=3,x_2=6$ and I know that $p_1=p_3$, then the MLE will probably be $p_1=0.25,p_2=0.5,p_3=0.25,x_3=3$. However, it is clearly not the maximum, since for example:

$L(p_1=0.24,p_2=0.52,p_3=0.24|x_1=3,x_2=6,x_3=2)=$

$=\frac{11!}{3!6!2!}0.24^30.52^60.24^2=0.07273$

So from this it seems that $x_1=3,x_2=6,x_3=2$ is more likely than $x_1=3,x_2=6,x_3=3$ even if I know that $p_1=p_3$, which seems very counter-intuitive.

My questions are whether my logic is sound, whether my intuition is misleading me and whether this is the correct way to estimate the parameters and missing data.

EDIT:

For future reference, I found this highly relevant paper, which addresses this exact problem. It also discusses the slight skew mentioned in whuber's answer.

Answer 1

One way to square your intuition with ML is to recognize that ML estimates are often biased. The ML estimate of $N$ looks like it's biased a little low.

There are only two parameters, $N$ and $p=p_1$, because $p_3=p_1=p$ and $p_2 = 1-p_1-p_3 = 1-2p$. The log likelihood for observations $(a,b)$ is

$$\log(\Lambda) = \log\binom{N}{a,b,N-a-b} + (N-b)\log(p) + b\log(1-2p)$$

which for any $N$ is maximized at the unique zero of its derivative, $p = (N-b)/(2N)$. It is apparent, too, that as a function of $N$ this is concave downward. Thus to obtain the MLE for $N$ we can scan over $N=a+b, a+b+1, \ldots$ until finding a maximum. (There are more efficient ways, but this works well enough.) For example, with $(a,b)=(3,6)$, the maximum occurs for $\widehat{N}=11$ where $\hat{p} = 5/22$.

Here is a histogram of $10^5$ iid draws of this MLE from a Multinomial$(12; 1/4, 1/2, 1/4)$ distribution:

The shift to a peak at $\widehat{N}=11$ is clear. For larger $N$ the shift still appears to be leftwards by only a small amount. Here is a histogram from a simulation with a Multinomial$(120; 1/4, 1/2, 1/4)$ distribution:

The bias looks like a shift of $1$ or $2$ leftwards (the peak is at $119$ and the mean is $118.96$), but certainly there is not a proportional shift to $11/12 * 120 = 110$.

Although (because $N$ must be integral) some care should be used in applying standard ML results, the mathematical formulation of the likelihood makes sense for non-integral $N$ (via Gamma functions), so you're probably ok using the usual ML-based confidence intervals, etc.