Robust meanΒΆ

The computation of robust mean is based on a re-descending Psi type M-estimator, that uses Hampel’s influence function, see [hampel86] and [huber81]:

(1)\Psi(x) = \begin{cases} x & \textrm{if } 0 \leq |x| \leq a \\ a \operatorname{sign}(x) & \textrm{if } a \leq |x| \leq b \\ \frac{a(c-|x|)}{(c-b)}\operatorname{sign}(x) & \textrm{if } b \leq |x| \leq c \\ 0 & \textrm{if } c \leq |x| \\ \end{cases}

where a = 1.7, b = 3.4 and c = 8.5. The function \Psi(x) and its first derivative \Psi'(x) is plotted on Figure Hampel’s influence function. (a) influence function ; (b) first derivative .

Hampel's influence function

Hampel’s influence function. (a) influence function \Psi; (b) first derivative \Psi'

The starting point for iterations is determined by median as an estimate of location and median of absolute deviations (MAD) for an estimation of scale s.

(2)s = MAD(X)/0.6745

The optimum solution is found by means of the Newton-Raphson optimizing algorithm. In iteration $j$ a new estimation of location \mu_{j} is determined using the following formula:

(3)\mu_{j} = \mu_{j-1} + s\,\frac{\sum \Psi(r_i)}{\sum \Psi'(r_i)}

where r_i is a residual

(4)r_i = \frac{x_i - \mu}{s}

If the input vector has normal distribution, the estimation of variance is:

(5)\sigma^2 = \frac{n}{n-1}\,n\,\frac{\sum\Psi^2(r_i)}{(\sum\Psi'(r_i))^2}

where n is the number of samples. The term \frac{n}{n-1} is Bessel’s correction.

The iterations are stopped when one of the following conditions is satisfied:

  • The number of iterations exceed 100 or
  • the absolute difference |\mu_{j} - \mu_{j-1}| is less than 10^{-4}\sigma, where \sigma is an square root of variance \sigma^2 computed at each step or
  • the absolute difference |\mu_{j} - \mu_{j-1}| is less than 10^{-7}.