How to find maximum likelihood estimates of an integer parameter?

Question

H.W. Question:

$x_1,x_2,\ldots,x_n$ are independent Gaussian variables with mean $\mu$ and variance $\sigma^2$. Define $y = \sum_{n=1}^{N} x_n$ where $N$ is unknown. We are interested in estimation of $N$ from $y$.

a. Given $\hat N_1 = y/\mu$ determine its bias and variance.

b. Given $\hat N_2 = 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 \ne 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 $\hat N_1,\hat N_2$ attain CRLB when $N\to \infty$?

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

Thanks,

Nadav

Answer 1

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

Answer 2

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{equation} \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} \end{equation}$$

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

$$\begin{equation} \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} \end{equation}$$

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

$$\begin{equation} \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} \end{equation}$$

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

Answer 3

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