# How to find a maximum likelihood estimator for a discrete parameter?

Let $X=\{x_1,x_2,...x_m\}$ be a random variable with density $f(x)=cx^{2n}$ if $-1\lt x\lt 1$, $0$ otherwise, and $n \in \{ 1,2,3,4\}$. Find the MLE of $n$.

I tried to find the Likelihood function : $(n+1/2)^m\prod_{i=1}^mx_i^{2n}$ but I'm not sure. Suggestions? Thanks

• Given the likelihood function, all that's left to do is to maximize it with respect to n. Do you know how to do that? Hint: What are the possible values of n which have to be considered when maximizing the likelihood function? – Mark L. Stone May 7 '16 at 22:03
• yes, Mark L. Stone, but is correct my likelihood function? – albert May 8 '16 at 15:43

This problem is simpler than it might look: although it might get confusing when one tries to apply routine Calculus methods, it is easy when worked from general principles.

By definition, the likelihood $\mathcal L$ is the probability of the data. Since the data are (implicitly) assumed independent, this is the product of the individual probability densities, each equal to $(n+1/2)(x_i^2)^n$. Consequently, as shown in the question,

$$\mathcal{L}(n) = \prod_{i=1}^m \left((n+1/2)x_i^{2n}\right) = (n+1/2)^m \left(\prod_{i=1}^m x_i^2\right)^n.$$

The only part depending on the data is that product on the right. Since it is almost surely positive for any of the $n$ under consideration, we may write it in terms of its logarithm. A convenient multiple of that is the mean log of the squared data, $$q=\frac{1}{m}\log\prod_{i=1}^m x_i^2.$$ Note that $q$ depends on the data, but not on the unknown parameter $n$: it will be the test statistic on which the Maximum Likelihood estimate is based.

Indeed, let's take logs of both sides, obtaining

$$\log\mathcal{L}(n) = m\log(n+1/2) + mnq.$$

Observe that since $|x_i| \lt 1$ for all $i$, $q \lt 0$. To maximize $\mathcal{L}$ we merely have to select the largest of the four values with $n=1,2,3,4$: namely,

\eqalign{\log \mathcal{L}(1) &= m\log(3/2) + mq,\\ \log \mathcal{L}(2) &= m\log(5/2) + 2mq,\\ \log \mathcal{L}(3) &= m\log(7/2) + 3mq,\\ \log \mathcal{L}(4) &= m\log(9/2) + 4mq.}

From row to row the $m\log(n+1/2)$ terms increase, but they do so more and more slowly; yet the $nmq$ terms decrease the value of $\log\mathcal L(n)$ at a constant rate of $m|q|$. Thus, the largest value of $(1/m)\log\mathcal L$--which corresponds to the largest value of $\mathcal L$ itself--is attained at the point where, in scanning these four values, $|q|$ first exceeds $$\log((n+1)+1/2) - \log(n+1/2) = \log\left(\frac{2n+3}{2n+1}\right).$$

This leads to an extremely simple procedure that can be specified even before collecting the data. Namely,

The increasing sequence $\log(5/3) \gt \log(7/5) \gt \log(9/7)$ partitions the real numbers into four intervals $$I_1=(-\infty, \log\frac{9}{7}],\ I_2=[\log\frac{9}{7}, \log\frac{7}{5}],\ I_3=[\log\frac{7}{5}, \log\frac{5}{3}],\ I_4=[\log\frac{5}{3},\infty).$$ If $|q|\in I_j$, pick $\hat n=j$ as the maximum likelihood estimator.

If $|q|$ lies at the common boundary of two intervals, they give two MLEs for $n$.

For each $n=1,2,3,4$, this procedure was applied to $1000$ independent samples of size $10$, then $1000$ more independent samples of size $500$. The tabulation of the estimates shows decent correspondence between $n$ and $\hat n$ for the small samples--which is as much as one might hope--and near perfect correspondence for the large samples--which is where the Maximum Likelihood method ought to perform well. Here are the results (which can be reproduced with the following R code):

$m = 10 MLE=1 MLE=2 MLE=3 MLE=4 n=1 753 169 15 1 n=2 213 488 241 77 n=3 27 222 340 283 n=4 7 121 404 639$m = 500
MLE=1 MLE=2 MLE=3 MLE=4
n=1  1000     0     0     0
n=2     0   999     0     0
n=3     0     1   998     2
n=4     0     0     2   998


For instance, when the true parameter was $n=2$ and the sample size was only $m=10$, the top table shows the MLE was correct $488$ out of $1000$ times and was off by one (that is, estimating $n$ as either $1$ or $3$) $213+241$ other times. On the whole, the MLE appears a little biased towards the middle values for small $m$ and extremely accurate for large $m$.

#
# Generate random values.
#
rpow <- function(n, p) {
q <- runif(n)
i <- sample.int(2, n, replace=TRUE)
ifelse(i==1, -1, +1) * q^(1/(2*p+1))
}
#
# Confirm rpow works as intended by matching histograms to the PDFs.
#
par(mfrow=c(1,4))
for (n in 1:4) {
hist(rpow(1e4, n), main=paste("n =", n), freq=FALSE)
}
#
# Test the MLE.
#
MLE.pow <- function(x, breaks=log(c(5/3, 7/5, 9/7))) {
# Given data x, return the MLE in {1,2,3,4}.
m <- length(x)
q <- sum(log(x^2)) / m
sum(abs(q) <= breaks) + 1
}
m <- c(10,500) # Specify sample sizes to test
set.seed(17)   # Create a reproducible starting point
sim <- lapply(m, function(m) {
x <- replicate(1e3, sapply(1:4, function(n) MLE.pow(rpow(m, n))))
results <- apply(x, 1, tabulate, nbins=4) # Tabulate the MLEs
rownames(results) <- paste0("n=", 1:4)
colnames(results) <- paste0("MLE=", 1:4)
results
})
names(sim) <- paste("m =", m)
print(sim)

• (+1) I really like it when a forgotten question suddenly becomes filled with useful content! – Alecos Papadopoulos Jun 15 '17 at 18:55

This has a few tricky points, so let's work it out.

First, we note that $2n$ is an even number, so the function that wants to be a density will be non-negative from that respect, as it should, even though the variable $X$ may take negative values.

Second we need to determine the value of $c$ so that we have a proper density. We require

$$\int_{-1}^1 cx^{2n}dx =1 \implies \frac{c}{2n+1}x^{2n+1} \Big|^1_{-1} =1$$

and since $2n+1$ is an odd number, we get

$$\frac{c}{2n+1}[1-(-1)] =1 \implies c= n+1/2$$

Therefore

$$f_X(x) = (n+1/2)x^{2n}, \;\;\; -1<x<1$$

and the likelihood function from an i.i.d. sample is ($I\{\}$ being the indicator function)

$$L = I_{\{-1<x<1\}}\cdot (n+1/2)^m\prod_{i=1}^mx_i^{2n}$$

We usually consider the log-likelihood for various beneficial reasons, but here, if we take logarithms we will be looking at $\ln x_i$ which is not defined if $x_i\leq 0$. But if one attempted to maximize the likelihood directly to avoid taking logs it would again hit on the same problem -the logarithm of the $x$-values would again appear if we considered the derivative with respect to $n$.

But since $x_i^{2n}$ is non-negative, we can write

$$L = I_{\{-1<x<1\}}\cdot (n+1/2)^m\prod_{i=1}^m|x_i|^{2n}$$

$$\implies\ln L = \ln \left(I_{\{-1<x<1\}}\right) + m\ln(n+1/2)+2n\sum_{i=1}^m\ln|x_i|$$

..which leaves us only with the problem of obtaining a realized value $x_i=0$ exactly. Well, from a theoretic point of view we invoke the zero-probability of a continuous r.v. taking a specific value, while from an applied point of view, if our sample contains an exact zero value, we can just discard it.

While this log-likelihood is strictly concave in $n$ and has a straightforward f.o.c

$$\frac{\partial \ln L}{\partial n} = \frac {m}{n+1/2} + 2\sum_{i=1}^m\ln|x_i|$$

$$\implies \hat n = \frac {m}{-2\sum_{i=1}^m\ln|x_i|} - \frac 12$$

...we have just ignored the restrictions on the values of $n$. This is a discrete-optimization problem. Nothing guarantees that the above expression will give an integer, or that it will fall in between $1$ and $4$.

After a little thinking, instead of considering derivatives etc as usual and dance around $\hat n$, it would be simpler to evaluate the log-likelihood at the four possible values of $n$ and see for which it attains the highest value.