ML Estimator parameter of bimodal, polynomial distribution

Question

Let $X_1, \ldots X_n$ i.i.d with density of (using the indicator function $\mathbf 1$) $$f(x|\mu) = \frac 3 2 (x-\mu)^2 \cdot \mathbf 1_{[\mu-1, \mu+1]}(x) = \left\{ \begin{array}{ll} \frac{3}{2}(x-\mu)^2 & \quad \mu-1 \le x \le \mu+1 \\ 0 & \quad \textrm{otherwise} \end{array} \right.$$

The moment method estimator is easy to find, $\hat \mu_{MM} = \bar X_n$. I also want to find the maximum-likelihood estimator $\hat \mu_{ML}$, but I'm a bit stumped.

My thought process:

• Disregarding the $x \in [\mu - 1, \mu+1]$ constraint gives non-sensical results, suggesting that $\hat \mu_{ML} \rightarrow \pm \infty$
• The $|x-\mu| \le 1$ enforces that $\hat \mu_{ML} \in [\max_i(X_i)-1, \min_i(X_i)+1]$. That is a tight bound with increasing $n$, especially, as the distribution is weighted so heavily towards the edges. The problem is that "a tight bound" is not the estimator itself.
• Combining the two approaches, I find $$\hat \mu_{ML} = \operatorname{argmax}_{\mu} \left (\mathbf 1_{[\max_i(X_i)-1, \min_i(X_i)+1]}(\mu) \cdot \prod_i (X_i - \mu)^2\right)$$ Sadly, I can't simplify this further, or leverage it to calculate a value for $\hat \mu_{ML}$.
• Clearly, $\hat \mu^* = \frac 1 2 (\max_i(X_i) + \min_i(X_i))$ is a fairly reasonable estimator, but it does not seem to be the ML estimator, since it disregards the information contained in the $\prod$-term above.

Am I missing something, or is there really no closed-form expression for $\hat \mu_{ML}$? Any pointers would be appreciated.

Edit: Following @whuber's suggestion, I plotted $L(\mu | X_1 = 0, X_2 = 0.15, X_3=1.7, X_4=1.8)$:

As you can see, the maximum is not at the edges, but in the interior. It can be at the edges for different choices of the observed $X_i$. The same plot without the $|X_i - \mu| \le 1$ restriction with a slightly wider $\mu$ axis:

It's an 8th order polynomial for n=4. $\hat \mu$ would correspond to the maximum around ~0.9, but there are other maxima (one between each consecutive pairs of $X_i$).

Differentiating w.r.t. $\mu$ gives a sum-of-($2n-1$)-order polynomials, that I don't see how to solve. I'm a bit lost here, to be honest.

Answer 1

This is quite a tricky problem, and I think it will require a difficult optimisation of a function that will have multiple local maxima in general. Given an observed sample vector $\mathbf{x} = (x_1,...,x_n)$, your likelihood function is:

$$\begin{align} L_\mathbf{x}(\mu) &= \prod_{i=1}^n f(x_i|\mu) \\[6pt] &= \text{const.} \times \prod_{i=1}^n (\mu-x_i)^2 \cdot \mathbb{I}(\mu-1 \leqslant x_i \leqslant \mu+1) \\[6pt] &= \text{const.} \times \mathbb{I}(\max x_i-1 \leqslant \mu \leqslant \min x_i+1) \prod_{i=1}^n (\mu-x_i)^2. \\[6pt] \end{align}$$

If $\max x_i - \min x_i > 2$ you have a problem with your model specification --- we leave this case aside. Assuming your sample vector satisfies the requirement $\max x_i - \min x_i \leqslant 2$, the corresponding log-likelihood function over this range (removing the previous constant) is:

$$\begin{align} \ell_\mathbf{x}(\mu) &= 2 \sum_{i=1}^n \log|\mu-x_i|, \\[6pt] \end{align}$$

and at all values $\mu \notin \{ x_1,...,x_n \}$ the score and functions are given by:

$$\begin{align} s_\mathbf{x}(\mu) \equiv \frac{d \ell_\mathbb{x}}{d \mu}(\mu) &= \sum_{i=1}^n \frac{2}{\mu-x_i}, \\[6pt] I_\mathbf{x}(\mu) \equiv - \frac{d^2 \ell_\mathbb{x}}{d \mu^2}(\mu) &= \sum_{i=1}^n \frac{2}{(\mu-x_i)^2} >0. \\[6pt] \end{align}$$

We can also find the score values at the extremes of the allowable parameter range, which are:

$$\begin{align} s_{\mathbf{x}, \text{lower}} \equiv s_\mathbf{x}(\max x_i - 1) &= \sum_{i=1}^n \frac{2}{\max x_i - x_i - 1}, \\[6pt] s_{\mathbf{x}, \text{upper}} \equiv s_\mathbf{x}(\min x_i + 1) &= \sum_{i=1}^n \frac{2}{\min x_i - x_i + 1}. \\[6pt] \end{align}$$

The log-likelihood function has poles (with negitive infinte value) at all points $\mu \in \{ x_1,...,x_n \}$ so you need to examine the function on a piecewise basis in the regions between the points:

$$\max x_i-1 \leqslant x_{(1)} \leqslant \cdots \leqslant x_{(n)} \leqslant \min x_i+1.$$

Between each of its poles (occurring at the data values) the function is strictly concave so it will have a number of local maxima. Between the endpoints of the allowable parameter range and the nearest poles it is also strictly concave and may be monotonic. It is simple to compute the values $s_{\mathbf{x}, \text{lower}}$ and $s_{\mathbf{x}, \text{upper}}$ to see if the log-likelihood is monotonic in the end regions, which allows you to determine whether the local maximising point in that region is at a boundary. For all other regions the local maximum will be in the interior at the unique critical point.

Now, so far as I can see, the only way I can see to find the MLE is to compute each of the local maxima in the pieces separated by the values above and then compare them to see which is the biggest. In each region (except possibly the end regions) the MLE will be in the interior of the region at the unique critical point defined by the score equation $s_\mathbf{x}(\hat{\mu}) = 0$. This can be written more clearly as:

$$0 = \frac{1}{\hat{\mu}-x_{(1)}} + \cdots + \frac{1}{\hat{\mu}-x_{(n)}}.$$

This equation gives the local maximising point in each region; it will need to be solved numerically with the range constraint for each region. Once you compute all the local maxima you can compare them to see which is the biggest, which will give you the MLE.