# MLE of multinomial distribution with missing values

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.

• The reason that $[3,6,2]$ has a higher probability than $[3,6,3]$ is because it is less data. – jerad Jun 19 '13 at 18:27
• @jerad I think you are right, but I am not sure what this means. Does this mean it is incorrect to use MLE here? What would be the proper way of estimating $x_3$? – Bitwise Jun 19 '13 at 18:35
• I can't say anything constructive about that. If you can't observe $x_3$ or $N$, and you have no prior knowledge that constrains the possibilities then all you're left with is making a random guess, which isn't that useful. In any case, if your data sets are rather small (as in your examples) then I'd suggest doing Bayesian inference to estimate $p$. – jerad Jun 19 '13 at 19:04
• @jerad unfortunately my actual numbers are quite large. It seems to me that there should still be a way to fit $x_3$ without prior knowledge. The likelihood decreases as $x_3$ gets very large, and obviously it cannot be negative. BTW if $x_3=0$, there is even less data but the ML is still lower than $x_3=3$, so less data does not necessarily mean higher likelihood. – Bitwise Jun 19 '13 at 19:35
• If you're able to calculate the likelihood of $x_1 = 3$, $x_2=6$, $x_3=2$, you're using the wrong likelihood. – Glen_b Jun 20 '13 at 1:18

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.

• Great answer, thanks! This phenomenon is quite interesting. In my specific application, using Gamma functions makes sense. So you are saying that MLE is reasonable for a large $N$? I am trying to use an EM-type algorithm (for MLE) by alternating between estimating $N$ and estimating $p_i$ (which are more complicated that presented above). – Bitwise Jun 19 '13 at 20:50
• I don't see any reason why MLE is not reasonable for this problem. – whuber Jun 19 '13 at 21:28