Hampel 估计稳健均值——理论和Python实现

Posted 2021-06-29 zhuo木鸟

文章目录

• Hampel 方法①
• • 理论步骤
  • Python 代码
• Hampel 方法②
• • 理论步骤
  • Python 代码

本文的参考标准是：ISO 13528 的 C.5.3.2 和 C.5.3.3 条款。

设样本为 $x_1,x_2,\cdots ,x_n$，$n$ 为样本容量。

Hampel 方法①

Hampel 方法有两种实现方法，第一种方法的结果和具体参数选择有关。

理论步骤

1. 首先求出样本的中位数为 $x^{\ast }$；
2. 求出样本的稳健标准差 $s^{\ast }$，可用方法有：MADe、Qn/Q、或算法 A。
3. 对每一个样本点，求 $q_i$：
   $$q_i=|\frac{x_i-x^{\ast }}{s^{\ast }}|$$
4. 根据公式计算权重：
   $$w_i=\{\begin{array}{cc}0& |q|\ge c\\ \frac{a(c-q)}{q(c-b)}& b<|q|\le c\\ a/q& a<|q|\le b\\ 1& |q|\le a\end{array}$$
   其中，参数一般选择为：$a=1.5,b=3.0,c=4.5$，且 $a\le b\le c$，如果提高 $a,b,c$ 的范围，则可以提高该算法的效率；如果缩小范围，则可以提高算法抵抗离群值，和小“峰” 的能力。
5. 根据下述公式，重新计算 $x^{\ast }$：
   $$x^{\ast }=\frac{{\sum }_{i=1}^n{w}_i\cdot x_i}{{\sum }_{i=1}^n{w}_i}$$
   重复 3 到 5 环节，直到 $x^{\ast }$ 不变为止。所谓不变的条件（收敛条件）是：上一次的 $x^{\ast }$ 与当前的 $x^{\ast }$ 的差值不超过 $0.01s^{\ast }/\sqrt{n}$。

由于该 Hampel 方法的输出随着所选的参数，如 $a,b,c$ 以及 $s^{\ast }$ 的计算方法的不同而不同，所以一般要多多尝试，并选择最接近中位数的值，作为最终的输出结果。

Python 代码

def weight_fun(q, a=1.5, b=3.0, c=4.5):
    if a>b or b>c:
        # 若 a 大于 b 或者 b 大于 c，则会陷入死循环：
        raise(Exception('在调整权重函数时，参数必须满足：a<=b<=c'))
    
    if np.abs(q) > c:
        return 0
    elif np.abs(q) > b and np.abs(q) <= c:
        w = a*(c-q)/(q*(c-b))
        return w
    elif np.abs(q) > a and np.abs(q) <= b:
        w = a/q
        return w
    elif np.abs(q) <= a:
        return 1

def perform_hampel_reweighting_estimate(x, s_star=None, epochs=100, **tunning):
    '''
    根据标准 ISO 13528：2015 的 C.5.3.2 的需要迭代更新权重的 Hampel 算法求解稳健均值

    Parameters
    ----------
    x : Series
        样本数据.
    s_star : float
        标准差的稳健估计，参考标准，可用的估计有：MADe、Qn、Q 或者算法 A 等。
        若不提供，则默认使用 Q 作为稳健标准差的估计。
    epochs : int
        最大迭代步数，若超过该限度，则迭代自动终止
    weight_fun : function handler
        函数引用，默认值为 C.5.3.2 中步骤 v 规定的权重函数
    
    Returns
    -------
    x_star : float
        稳健平均值.

    '''
    # 迭代步数
    epoch = 0
    a = tunning.get('a', 1.5)
    b = tunning.get('b', 3.0)
    c = tunning.get('c', 4.5)

    if isinstance(x, pd.Series):
        x = x.astype('float').values
    
    # 样本个数 
    p = len(x)
    
    x_star = np.median(x)
    if s_star is None:
        # 若没有提供 s*，则默认采用 Q 估计作为 s*
        s_star = perform_Q_estimate(x)

    x_star_curr = x_star
    while True:
        # 步数加 1
        epoch += 1
        x_star_last = x_star_curr
        q_arr = (x-x_star)/s_star
        
        # 求权重
        w_arr = [weight_fun(wi, a=a, b=b, c=c) for wi in q_arr]
        w_arr = np.array(w_arr)
        
        # 按照权重更新 x
        numerator = sum(w_arr*x)
        denominator = sum(w_arr)
        x_star = numerator/denominator
        x_star_curr = x_star
        
        if np.abs(x_star_last - x_star_curr) <= 0.01*s_star/np.sqrt(p):
            return x_star_curr
        # 若达到最大迭代次数
        if epoch >= epochs:
            return x_star_curr

Hampel 方法②

理论步骤

这种 Hampel 方法求解的输出是唯一的，其原理在于解决一个方程，如下：
$$\sum _{i=1}^n\Psi (\frac{x_i-x^{\ast }}{s^{\ast }})=0$$
其中函数 $\Psi$ 如下：
Ψ ( q ) = { 0 0 ≥ q ≤ 4.5 − 4.5 − q − 4.5 < q ≤ − 3 − 1.5 − 3 < q ≤ − 1.5 q − 1.5 < q ≤ 1.5 1.5 1.5 < q ≤ 3 4.5 − q 3 < q ≤ 4.5 0 q > 4.5 \\Psi (q) =\\left\\{\\begin{array}{c} 0 & 0 \\geq q\\leq 4.5 \\\\ -4.5-q & -4.5 < q \\leq -3 \\\\ -1.5 & -3 < q \\leq -1.5\\\\ q & -1.5 < q \\leq 1.5 \\\\ 1.5 & 1.5 < q \\leq 3 \\\\ 4.5 - q & 3 < q \\leq 4.5 \\\\ 0 & q > 4.5 \\end{array}\\right. Ψ(q)=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪Trimmed 稳健均值估计与 中位数-中位数配对偏差法估计标准差——理论与 Python 实现

简单稳健估计法——原理与 Python 实现

简单稳健估计法——原理与 Python 实现

Q 方法求稳健标准差——理论和 Python实现

Q 方法求稳健标准差——理论和 Python实现

算法 A 求稳健平均值和稳健标准差