# ML Estimator parameter of bimodal, polynomial distribution

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.

• I think the situation might become obvious were you to graph the likelihood for a particular value (or values) of the $x_i$ and consider how to find its maximum.
– whuber
Sep 25, 2021 at 21:02
• @whuber: I plotted $L(\mu | X_1)$ for various $X_1$, and $L(\mu | X_1=0, X_2=1.7)$. The last was interesting, because it showed a maximum. For $n=2$, this gives a 4th order polynomial $L(\mu)$, order 3 after $d/d\mu$. No idea how to solve this for larger $n$. I got to $2\sum_i (x_i-\mu) \prod_{j \ne i} (x_j-\mu)^2 = 0$ but don't understand how to progress from here. Sep 25, 2021 at 22:04
• @StubbornAtom: I considered your $\hat \mu$, but all I can think of is plugging that into the argmax-$\prod$ above, and I don't see how that helps me. I guess I don't see where I have to start to turn your $c$-interpolation approach into a solution. Sep 25, 2021 at 22:05
• Hmm. Following whuber's suggestion, I guess the maximum occurs at the endpoints of the domain of the likelihood. Sep 26, 2021 at 19:03
• @StubbornAtom: That doesn't appear to be true, see the edit in my question. There can be one or more internal maxima. Each observation $x_i$ zero to the $f(\mu)$ polynomial, with maxima between the $\mu=x_i$ zeros. If the $x_i$ all lie close together (e.g. all in the left lobe of the distribution), then the maximum is not at any $\frac{\partial f}{\partial \mu}=0$ point, but at the edge of the domain, as you mention. Multiple maxima, or maxima at non-critical-points: I don't know how to tackle that analytically. Sep 26, 2021 at 19:59

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.

• Thank you for your answer! I see some problems: $\sum_i |\hat\mu-x_i|^{-1}$ is a sum of strictly positive values if $x_i, \hat \mu$ are finite, and thus can't be zero, numerical or not. More generally, the log-likelihood approach confuses things a lot because it introduces infinities. Yes, $I(x) > 0$, so one assumes $\mathcal l$ is concave - but that's wrong, since there is pole/negative infinity at every $\hat \mu=x_i$. Thus we can't assert monotony on $\mathcal l$ and -- as described in my comments on the question -- there can be multiple maxima within the lower/upper bound. Sep 27, 2021 at 6:28
• Sorry about that --- I had a mathematical error. I have now edited to correct this.
– Ben
Sep 27, 2021 at 8:46