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 \Psi; (b) first derivative \Psi'.

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}.