Unbiased median absolute deviation based on the Harrell-Davis quantile estimator

The median absolute deviation (\(\textrm{MAD}\)) is a robust measure of scale. In the previous post, I showed how to use the unbiased version of the \(\textrm{MAD}\) estimator as a robust alternative to the standard deviation. “Unbiasedness” means that such estimator’s expected value equals the true value of the standard deviation. Unfortunately, there is such thing as the bias–variance tradeoff: when we remove the bias of the \(\textrm{MAD}\) estimator, we increase its variance and mean squared error (\(\textrm{MSE}\)).

In this post, I want to suggest a more efficient unbiased \(\textrm{MAD}\) estimator. It’s also a consistent estimator for the standard deviation, but it has smaller \(\textrm{MSE}\). To build this estimator, we should replace the classic “straightforward” median estimator with the Harrell-Davis quantile estimator and adjust bias-correction factors. Let’s discuss this approach in detail.

Introduction

Let’s consider a sample \(x = \{ x_1, x_2, \ldots, x_n \}\). Its median absolute deviation \(\textrm{MAD}_n\) can be defined as follows:

\[\textrm{MAD}_n = C_n \cdot \textrm{median}(|X - \textrm{median}(X)|) \]

where \(\textrm{median}\) is a median estimator, \(C_n\) is a scale factor (or consistency constant).

Typically, \(\textrm{median}\) assumes the classic “straightforward” median estimator:

• If \(n\) is odd, the median is the middle element of the sorted sample
• If \(n\) is even, the median is the arithmetic average of the two middle elements of the sorted sample

However, we can use other median estimators. Let’s consider \(\textrm{median}_{\textrm{HD}}\) which calculates the median using the Harrell-Davis quantile estimator (see [Harrell1982]):

\[\textrm{median}_{\textrm{HD}}(x) = \sum_{i=1}^n W_i x_i, \quad W_i = I_{i/n} \bigg( \frac{n+1}{2}, \frac{n+1}{2} \bigg) - I_{(i-1)/n} \bigg( \frac{n+1}{2}, \frac{n+1}{2} \bigg) \]

where \(I_t(a, b)\) denotes the regularized incomplete beta function.

When \(n \to \infty\), we have the exact value of \(C_n\) which makes \(\textrm{MAD}_n\) an unbiased estimator for the standard deviation:

\[C_n = C_\infty = \dfrac{1}{\Phi^{-1}(3/4)} \approx 1.4826022185056 \]

For finite values of \(n\), we should use adjusted \(C_n\) values. We already discussed these values for the straightforward median estimator in the previous post (this problem is well-covered in [Park2020]). For \(\textrm{median}_{\textrm{HD}}\), we should use another set of \(C_n\) values. We reuse the approach from [Hayes2014] with the following notation:

\[C_n = \dfrac{1}{\hat{a}_n} \]

Factors for small n

Factors for the small values of \(n\) can be obtained using numerical experiments. Here is the scheme for the performed simulation:

• Enumerate \(n\) from \(2\) to \(100\)
• For each \(n\), generate \(2\cdot 10^8\) samples from the normal distribution of size \(n\) (the Box–Muller transform is used, see [Box1958])
• For each sample, estimate \(\textrm{MAD}_n\) using \(C_n = 1\)
• Calculate the arithmetic average of all \(\textrm{MAD}_n\) values for the same \(n\): the result is the value of \(\hat{a}_n\)

As the result of this simulation, I got the following table with \(\hat{a}_n\) and \(C_n\) values for \(n \in \{ 2..100 \}\).

Here is a visualization of this table:

Factors for huge n

To build the equation for \(n > 100\), we also continue the approach suggested in [Hayes2014] (pages 2208-2209) and use the prediction equation of the following form:

\[\hat{a}_n = \Phi^{-1}(3/4) \Bigg( 1 - \dfrac{\alpha}{n} - \dfrac{\beta}{n^2} \Bigg). \]

We simulated approximations of \(\hat{a}_n\) for some large \(n\) values:

Next, we fitted these values using the multiple linear regression. The dependent variable is \(y = 1 - \hat{a}_n / \Phi^{-1}(3/4)\); the independent variables are \(n^{-1}\) and \(n^{-2}\); the intercept is zero:

\[y = \alpha n^{-1} + \beta n^{-2}. \]

The fitted model was adjusted a little bit to get nice-looking values (keeping the good accuracy level):

\[\alpha = 0.5,\quad \beta = 6.5. \]

These values give pretty accurate estimation of \(\hat{a}_n\) for \(n > 100\) (the accuracy is less than \(10^{-4}\)).

The mean squared error

It’s time to compare the straightforward and the Harrell-Davis-based \(\textrm{MAD}\) estimators in terms of the mean squared error. To estimate the \(\textrm{MSE}\), we also use numerical simulations. Here are the results for some \(n\) values:

As we can see, the difference between estimators for small samples is tangible. For example, for \(n = 3\), we got \(\textrm{MSE} \approx 0.69\) for the straightforward estimator vs. \(\textrm{MSE} \approx 0.27\) for the Harrell-Davis-based estimator. Below you can see density plots of \(\textrm{MSE}\) for the above \(n\) values.

Conclusion

In this post, we discussed the unbiased \(\textrm{MAD}_n\) estimator which is based on the Harrell-Davis quantile estimator instead of the straightforward median estimator. The suggested estimator is much more efficient for small samples because it has smaller \(\textrm{MSE}\). It allows using \(\textrm{MAD}_n\) as a efficient alternative to the standard deviation under the normal distribution with higher accuracy.

References

• [Harrell1982]
  Harrell, F.E. and Davis, C.E., 1982. A new distribution-free quantile estimator. Biometrika, 69(3), pp.635-640.
  https://doi.org/10.2307/2335999
• [Hayes2014]
  Hayes, Kevin. “Finite-sample bias-correction factors for the median absolute deviation.” Communications in Statistics-Simulation and Computation 43, no. 10 (2014): 2205-2212.
  https://doi.org/10.1080/03610918.2012.748913
• [Park2020]
  Park, Chanseok, Haewon Kim, and Min Wang. “Investigation of finite-sample properties of robust location and scale estimators.” Communications in Statistics-Simulation and Computation (2020): 1-27.
  https://doi.org/10.1080/03610918.2019.1699114
• [Box1958]
  Box, George EP. “A note on the generation of random normal deviates.” Ann. Math. Statist. 29 (1958): 610-611.
  https://doi.org/10.1214/aoms/1177706645