Robust estimator of parameters of a mixture distribution

I have a random variable taking one of three values in a arithmetic progression, $a$, $a+b$, $a+2b$, with probabilities $1/4,1/2,1/4$, where $b\approx 0.5a$ to $2a$. It is perturbed by a little noise (say Gaussian with deviation equal to a tenth of $a$), and there can be a small proportion of outliers, say $10\%$.

I would like to build a robust estimator of $a$ and $b$ from a sample of $N$ values, where $N$ is of the order of $50$.

I could for instance try and find three modes in the histogram, then classify the values by assigning them to the closest mode and estimate the class averages. I also assume that the median would be a good estimator of $a+b$. (And the distribution of $|v-\overline{a+b}|$ would be bimodal.)

Can you advise an approach for $b$ ?

Following the comments by @whuber, I can recast the question as estimating the means of a three membered Gaussian mixture, knowing that these means are equidistant, the variances are equal and there is a good separation between the members (but a few outliers).

• (1) This is an arithmetic progression, not a geometric one. (2) The distribution is continuous: it's not discrete. A more standard description would be that it is a finite mixture distribution. Knowing what to call these elements of your problem may help in identifying related problems and searching for solutions. (3) Do you perhaps mean that the Gaussian SD is a tenth of $b$? The size of $a$ isn't relevant because $a$ appears only as a location parameter, whereas the size of $b$ matters because it is a scale parameter.
– whuber
Commented Apr 6, 2018 at 14:59
• @whuber: (1) yep, typo. (2) this is why I said quasi-discrete. You are right, this is a mixture. (3) No, of $a$. But in the given problem, $a$ and $b$ are of the same order of magnitude and the scale is given collectively by $a$ and $b$. Commented Apr 6, 2018 at 15:58
• @whuber: (1) yep, typo. (2) this is why I said quasi-discrete. You are right, this is a mixture. (3) No, of $a$. But in the given problem, $a$ and $b$ are of the same order of magnitude and the scale is given by both. Notice that I changed the range of the ratio, which was incorrect. Commented Apr 6, 2018 at 16:05

A cluster analysis, followed by finding the medians of the clusters, would likely work well. However, there are demonstrably robust, effective, simple procedures that could be recommended because (1) they are relatively easy to study and analyze; (2) they capitalize on the assumed arithmetic progression and the assumption that about half the data are near the middle mode; and (3) they can be tuned to achieve a desired amount of robustness (within limits).

One procedure begins by estimating $a+b$ as the median of the data. You then expect the absolute residuals to be clustered, in approximately equal proportions, near $0$ and $b$. By specifying a suitable quantile $q$ between $0.5$ and $1.0$ you can obtain a preliminary estimate $b^{*}$ of $b$: for instance, the $0.75$ quantile should be a good estimate. With this preliminary estimate in hand, complete a simple cluster analysis by partitioning the absolute residuals into those less than $b^{*}$ and those greater than $b^{*}$. The median of the latter group is a robust estimate of $b.$ The estimate of $a$ is obtained by subtracting the estimate of $b$ from the estimate of $a+b.$ (See the code at the end of this post.)

There's a trade-off: if $q$ isn't high enough, there may be an appreciable chance that this quantile will lie in the cluster near $0$ and thereby grossly underestimate $b$. If $q$ is too high, the quantile will be sensitive to outliers and perhaps grossly overestimate $b.$ Adjusting $q$ to achieve a desired balance between these behaviors allows you to tune the procedure.

Here's a quick analysis. The chance that $b^{*}$ is too small in a dataset of size $n$ is approximately the chance that more than $qn$ of the data come from the component located at $a+b$. Ignoring outliers--most of which are likely not to be close to $a+b,$ anyway--this has a Binomial$(n, 1/2)$ distribution. For instance, with $n=50$ and $q=0.75$, the chance is only $0.00015.$ That may be appreciable enough (it happened once to me in a set of $200$ simulated datasets) that you might want to increase $q$ to somewhere between $0.8$ and $0.9.$ That will, however, reduce the breakdown point: with $q=0.9,$ for instance, more than five outliers can radically change $b^{*}.$

I performed a simulation with $n=50$ and $q=0.875,$ using perturbations equal to $0.1b$ times a Cauchy distribution. It's a fairly severe test because this distribution produces outliers at a great rate: it has about a $1$ in $8$ chance of changing a value by more than $b/2,$ thereby moving it closer to another mode (or far from all the modes).

Here are histograms of the first $20$ iterations. I fixed $a=b=1$ and have marked the locations $a=1,a+b=2,a+2b=3$ with vertical black dashes. The estimates are shown with red lines. They are in good agreement with the true values for these iterations. (The choices of $a$ and $b$ don't matter: all that's relevant is the dispersion of the random noise, expressed on a scale where $b$ is one unit.)

The outliers are evident in iterations 5, 12, and 14 because they are so extreme. You can also see plenty of data falling vaguely between two of the modes, making them difficult to classify.

The next figure is a scatterplot to show the relationship between actual values of the modes (that is, $1,2,3$) and the estimated values from all 2000 iterations.

Evidently none of the estimates was bad and they tend to be unbiased. The vast majority were well within $0.1$ of the true values. (Over $92\%$ of the estimates of $a$ were between $0.9$ and $1.1$ while almost $99\%$ of the estimates of $b$ were between $0.9$ and $1.1$.)

Finally, let's look at how the estimates of $a$ and $b$ are related, by drawing their scatterplot:

As you might expect, they are negatively correlated. (The curve is a local nonlinear smooth, so it's remarkable that it is nearly linear throughout the full range of estimates.) Because the cloud is concentrated at $(1,1)=(a,b),$ the estimates are unbiased. (Their averages in this simulation were $1.009$ and $0.993,$ respectively, with a correlation coefficient of $-60\%.$)

Finally, the R code to create these estimates is fast and simple. It returns estimates of all three modes. It is guaranteed that the middle estimate is exactly halfway between the two extreme estimates. (This differs from a general procedure to estimate the locations of a mixture of three components. It also differs insofar as it relies on the assumption that the weight of the middle component is one-half the total.)

estimator <- function(x, q=0.875) {      # x is an array of data
ab.hat <- median(x)                    # Estimate of a+b
z <- abs(x - ab.hat)                   # Absolute residuals
b.hat <- quantile(z, q)                # Preliminary estimate of b
threshold <- b.hat / 2                 # Cutoff between 0 and b
b.hat <- median(z[z >= threshold])     # Final estimate of b
c(ab.hat-b.hat, ab.hat, ab.hat+b.hat)
}


Quick proposal. Feel free to improve it.

Start from your partial solution, that is, let $m_2$ be your overall median. Since you know that approximately half of the point should be in the middle group, removing half of the points, the ones that are closest to your overall median, could be a good idea.

At this point you will be left with the other half, which should be divided into two groups of approximately the same size. Taking the median of all point that are less than $m_2$, name it $m_1$, and then the median of all points that are bigger than $m_2$, name it $m_3$, is an interesting option.

Finally, let $(a,b) = \mathrm{argmin} \sum_{i=1}^3 \{a+b(i-1) - m_i\}^2$.

That solution might not be the best if your groups do not possess similar variances, or huge variances so that they overlap a lot.

• I like the idea of removing the middle group. As there won't be exactly $n/2$ elements in it, they can appear as outliers in the lateral groups, or "steal" a bit of these, but taking the median of the lateral groups can cope. In fact, these estimators are the elements of rank $n/8,n/2$ and $7n/8$. I plan to use them as initial values for a k-medioids process. Commented Apr 6, 2018 at 17:15