# How to find Maximum Likelihood estimates of an *integer* parameter?

If anyone could direct me to the solution of the following problem it would be great.

H.W Question:

$$x_1$$,$$x_2$$,...,$$x_n$$ are Gaussian independent variables with with mean $$\mu$$ and variance $$\sigma^2$$.

$$y$$ is the sum of these $$N$$ random variables -> $$y = \sum_{n=1} ^{N} {x_n}$$

$$N$$ is unknown.

We are interested in estimation of $$N$$ from $$y$$.

a. given $$N_{1}{estimated} = y/\mu$$ determine its bias and variance.

b. given $$N_{2}{estimated} = y^2/\sigma^2$$ determine its bias and variance.

Ignoring the requirement for N to be an integer

c. Is there an efficient estimator (look at both $$\mu = 0$$ and $$\mu != 0$$)

d. Find the maximum likelihood estimate of N from y.

e. Find CRLB of N from y.

f. Does the mean squared error of the estimators $$N_{1}{estimated}, N_{2}{estimated}$$ get CRLB when $$N->\inf$$

Thanks,

• What is the distribution of $Y = \sum_i X_i\,?$ – BruceET Jun 12 at 14:32
• It doesn't say. I suppose it will also be distributed as Gaussian variable since it's a sum of Gaussian variables – Nadav Talmon Jun 13 at 11:37
• If $X_i$ are normal, then $Y = \sum_I Xi$ and $\hat N = Y/\mu$ are normal. What are mean and variance of $\hat N\,?$ That should finish the problem. // In practice, I suppose it makes sense to round $\hat N$ to an integer. That could make a slight difference in mean and variance. You could find out how much difference by simulation. – BruceET Jun 13 at 14:07
• Wouldn't the $Var(N_{estimated})$ will be the $Var(y)/\mu$? Same logic for the mean – Nadav Talmon Jun 13 at 14:24
• Because $N$ is integral, you can't (directly) use Calculus to find the minimum. If this is your obstacle, then please present your work in your question so we can focus on where you actually need help. – whuber Jun 13 at 14:52

You started well by writing down an expression for the likelihood. It is simpler to recognize that $$Y,$$ being the sum of $$N$$ independent Normal$$(\mu,\sigma^2)$$ variables, has a Normal distribution with mean $$N\mu$$ and variance $$N\sigma^2,$$ whence its likelihood is

$$\mathcal{L}(y,N) = \frac{1}{\sqrt{2\pi N\sigma^2}} \exp\left(-\frac{(y-N\mu)^2}{2N\sigma^2}\right).$$

Let's work with its negative logarithm $$\Lambda = -\log \mathcal{L},$$ whose minima correspond to maxima of the likelihood:

$$2\Lambda(N) = \log(2\pi) + \log(\sigma^2) + \log(N) + \frac{(y-N\mu)^2}{N\sigma^2}.$$

We need to find all whole numbers that minimize this expression. Pretend for a moment that $$N$$ could be any positive real number. As such, $$2\Lambda$$ is a continuously differentiable function of $$N$$ with derivative

$$\frac{d}{dN} 2\Lambda(N) = \frac{1}{N} - \frac{(y-N\mu)^2}{\sigma^2N^2} - \frac{2\mu(y-N\mu)}{N\sigma^2}.$$

Equate this to zero to look for critical points, clear the denominators, and do a little algebra to simplify the result, giving

$$\mu^2 N^2 + \sigma^2 N -y^2 = 0\tag{1}$$

with a unique positive solution (when $$\mu\ne 0$$)

$$\hat N = \frac{1}{2\mu^2}\left(-\sigma^2 + \sqrt{\sigma^4 + 4\mu^2 y^2}\right).$$

It's straightforward to check that as $$N$$ approaches $$0$$ or grows large, $$2\Lambda(N)$$ grows large, so we know there's no global minimum near $$N\approx 0$$ nor near $$N\approx \infty.$$ That leaves just the one critical point we found, which therefore must be the global minimum. Moreover, $$2\Lambda$$ must decrease as $$\hat N$$ is approached from below or above. Thus,

The global minima of $$\Lambda$$ must be among the two integers on either side of $$\hat N.$$

This gives an effective procedure to find the Maximum Likelihood estimator: it's either the floor or the ceiling of $$\hat N$$ (or, occasionally, both of them!), so compute $$\hat N$$ and simply choose which of these integers makes $$2\Lambda$$ smallest.

Let's pause to check that this result makes sense. In two situations there is an intuitive solution:

1. When $$\mu$$ is much greater than $$\sigma$$, $$Y$$ is going to be close to $$\mu,$$ whence a decent estimate of $$N$$ would simply be $$|Y/\mu|.$$ In such cases we may approximate the MLE by neglecting $$\sigma^2,$$ giving (as expected) $$\hat N = \frac{1}{2\mu^2}\left(-\sigma^2 + \sqrt{\sigma^4 + 4\mu^2 y^2}\right) \approx \frac{1}{2\mu^2}\sqrt{4\mu^2 y^2} = \left|\frac{y}{\mu}\right|.$$

2. When $$\sigma$$ is much greater than $$\mu,$$ $$Y$$ could be spread all over the place, but on average $$Y^2$$ should be close to $$\sigma^2,$$ whence an intuitive estimate of $$N$$ would simply be $$y^2/\sigma^2.$$ Indeed, neglecting $$\mu$$ in equation $$(1)$$ gives the expected solution $$\hat N \approx \frac{y^2}{\sigma^2}.$$

In both cases, the MLE accords with intuition, indicating we have probably worked it out correctly. The interesting situations, then, occur when $$\mu$$ and $$\sigma$$ are of comparable sizes. Intuition may be of little help here.

To explore this further, I simulated three situations where $$\sigma/\mu$$ is $$1/3,$$ $$1,$$ or $$3.$$ It doesn't matter what $$\mu$$ is (so long as it is nonzero), so I took $$\mu=1.$$ In each situation I generated a random $$Y$$ for the cases $$N=2,4,8,16,$$ doing this independently five thousand times.

These histograms summarize the MLEs of $$N$$. The vertical lines mark the true values of $$N$$.

On average, the MLE appears to be about right. When $$\sigma$$ is relatively small, the MLE tends to be accurate: that's what the narrow histograms in the top row indicate. When $$\sigma \approx |\mu|,$$ the MLE is rather uncertain. When $$\sigma \gg |\mu|,$$ the MLE can often be $$\hat N=1$$ and sometimes can be several times $$N$$ (especially when $$N$$ is small). These observations accord with what was predicted in the preceding intuitive analysis.

The key to the simulation is to implement the MLE. It requires solving $$(1)$$ as well as evaluating $$\Lambda$$ for given values of $$Y,$$ $$\mu,$$ and $$\sigma.$$ The only new idea reflected here is checking the integers on either side of $$\hat N.$$ The last two lines of the function f carry out this calculation, with the help of lambda to evaluate the log likelihood.

lambda <- Vectorize(function(y, N, mu, sigma) {
(log(N) + (y-mu*N)^2 / (N * sigma^2))/2
}, "N") # The negative log likelihood (without additive constant terms)

f <- function(y, mu, sigma) {
if (mu==0) {
N.hat <- y^2 / sigma^2
} else {
N.hat <- (sqrt(sigma^4 + 4*mu^2*y^2) - sigma^2) / (2*mu^2)
}
N.hat <- c(floor(N.hat), ceiling(N.hat))
q <- lambda(y, N.hat, mu, sigma)
N.hat[which.min(q)]
} # The ML estimator

• I could not have asked for a better explanation. Thank you very much, you literally covered everything! – Nadav Talmon Jun 14 at 4:52
• now I need to tell whether an efficient estimator exists(for $\mu != 0$ and $\mu = 0$). I know that if an estimator is unbiased and answers the CRLB than its efficient. I know its unbiased, but taking the second derivative of the L function seems to get me no where. – Nadav Talmon Jun 14 at 6:27
• Ignore the fact that $N$ is integral: that is, allow the estimate to be the global minimum of the negative log likelihood function. Go on from there. – whuber Jun 14 at 12:39
• I took the derivative of the negative log likelihood function as you suggested and attempted to obtain the following expression: $C(N)\times(g(y)-N)$ I managed to do that on $\mu = 0$ but not on $\mu != 0$ this is why we had a unique positive solution on $\mu != 0$? – Nadav Talmon Jun 15 at 9:40
• I don't think so. I find it easier to reparameterize the problem in terms of $\theta=1/N,$ because then the derivative of the log likelihood is a quadratic function of $\theta.$ – whuber Jun 15 at 12:04

The method whuber has used in his excellent answer is a common optimisation "trick" that involves extending the likelihood function to allow real values of $$N$$, and then using the concavity of the log-likelihood to show that the discrete maximising value is one of the discrete values on either side of a continuous optima. This is one commonly used method in discrete MLE problems involving a concave log-likelihood function. Its value lies in the fact that it is usually possible to get a simple closed-form expression for the continuous optima.

For completeness, in this answer I will show you an alternative method, which uses discrete calculus using the forward-difference operator. The log-likelihood function for this problem is the discrete function:

$$\ell_y(N) = -\frac{1}{2} \Bigg[ \ln (2 \pi) + \ln (\sigma^2) + \ln (N) + \frac{(y-N\mu)^2}{N\sigma^2} \Bigg] \quad \quad \quad \text{for } N \in \mathbb{N}.$$

The first forward-difference of the log-likelihood is:

\begin{aligned} \Delta \ell_y(N) &= -\frac{1}{2} \Bigg[ \ln (N+1) - \ln (N) + \frac{(y-N\mu - \mu)^2}{(N+1)\sigma^2} - \frac{(y-N\mu)^2}{N\sigma^2} \Bigg] \\[6pt] &= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+1}{N} \Big) + \frac{N(y-N\mu - \mu)^2 - (N+1)(y-N\mu)^2}{N(N+1)\sigma^2} \Bigg] \\[6pt] &= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+1}{N} \Big) + \frac{[N(y-N\mu)^2 -2N(y-N\mu) \mu + N \mu^2] - [N(y-N\mu)^2 + (y-N\mu)^2]}{N(N+1)\sigma^2} \Bigg] \\[6pt] &= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+1}{N} \Big) - \frac{(y + N \mu)(y-N\mu) - N \mu^2}{N(N+1)\sigma^2} \Bigg]. \\[6pt] \end{aligned}

With a bit of algebra, the second forward-difference can be shown to be:

\begin{aligned} \Delta^2 \ell_y(N) &= -\frac{1}{2} \Bigg[ \ln \Big( \frac{N+2}{N} \Big) + \frac{2 N (N+1) \mu^2 + 2(y + N \mu)(y-N\mu)}{N(N+1)(N+2)\sigma^2} \Bigg] < 0. \\[6pt] \end{aligned}

This shows that the log-likelihood function is concave, so its smallest maximising point $$\hat{N}$$ will be:

\begin{aligned} \hat{N} &= \min \{ N \in \mathbb{N} | \Delta \ell_y(N) \leqslant 0 \} \\[6pt] &= \min \Big\{ N \in \mathbb{N} \Big| \ln \Big( \frac{N+1}{N} \Big) \geqslant \frac{(y + N \mu)(y-N\mu) - N \mu^2}{N(N+1)\sigma^2} \Big\}. \end{aligned}

(The next value will also be a maximising point if and only if $$\Delta \ell_y(\hat{N}) = 0$$.) The MLE (either the smallest, or the whole set) can be programmed as a function via a simple while loop, and this should be able to give you the solution pretty quickly. I will leave the programming part as an exercise.

• I appreciate your time and the thorough explanation. Thank you @Ben ! – Nadav Talmon Jun 15 at 12:52

Comment: Here is a brief simulation in R for $$\mu = 50, \sigma = 3,$$ which should be accurate to 2 or three places, approximating the mean and SD of $$Y.$$ You should be able to find $$E(Y)$$ and $$Var(Y)$$ by elementary analytic methods as indicated in my earlier Comment. If we had $$N = 100$$ then $$E(\hat N)$$ seems unbiased for $$N.$$

N = 100;  mu = 50;  sg = 3
y = replicate( 10^6, sum(rnorm(N, mu, sg))/mu )
mean(y);  sd(y)
[1] 99.99997
[1] 0.6001208
N.est = round(y);  mean(N.est);  sd(N.est)
[1] 99.9998
[1] 0.6649131

• Thank you Bruce! – Nadav Talmon Jun 13 at 14:27
• May I ask one more question? Now i'm asked if there is an efficient estimator regarding what I found, it also states that we now ignore the requirement for N to be an integer. what does it mean that it is not an integer anymore? how would I find the log Likelihood for such a case? – Nadav Talmon Jun 13 at 14:28
• If you would like to understand the potential for bias, don't use large $N:$ try a small value. $N=1$ is especially interesting :-). So is the case $\mu=0.$ – whuber Jun 13 at 14:41