双全估计算均值和方差——理论与 Python 实现

Posted 2021-06-30 zhuo木鸟

文章目录

• 均值估计（location estimate）
• • 步骤
  • Python 实现
• 方差估计（scale estimate）
• • 步骤
  • Python 代码

本文主要参考标准：16269-4:2010 的 5.2.3 条款和 5.3.3 条款

设样本为 $x_1,x_2,\cdots ,x_n$。

均值估计（location estimate）

步骤

双权（biweight）均值估计主要用在样本采样与一个不对称分布，并存在离群值的情景。他是通过迭代方法求解均值的，迭代过程如下：
$$T_n=\frac{{\sum }_{|u_i|<1}{x}_i(1-u_i^2{)}^2}{{\sum }_{|u_i|<1}(1-u_i^2{)}^2}$$
其中：
$$u_i=\frac{x_i-T_n}{c\cdot M_{ad}}$$
其中 $M_{ad}=\text{median}(|x_i-\text{median}(x)|)$。$T_n$ 的初始值可以取为样本中位数。迭代停止的条件是：$T_n$ 和 $T_{n+1}$ 基本不变为止。此时 $T_n$ 就是稳健均值。

注：系数 $c$ 一般取为 6。意味着样本的总体为正态分布时，那些落在 $\text{median}(x)\pm 4\times \sigma$ 之外的点，参与计算 $T_n$ 是，其权重为 0 。$\sigma$ 为标准差。

Python 实现

def perform_beweight_location_estimate(x, c=6, epochs=100):
    if isinstance(x, pd.Series):
        x = x.astype('float').values
    
    x_med = np.median(x)
    T_last = x_med
    epoch = 1
    while True:
        mad = np.median(np.abs(x-x_med))
        u_arr = (x-T_last)/(c*mad)
        idx = np.abs(u_arr) < 1
        
        coeff = (1-u_arr[idx]**2)**2
        
        Tn = sum(x[idx]*coeff)/sum(coeff)
        T_curr = Tn
        if np.abs(T_curr - T_last) <= 1e-5:
            return T_curr
        epoch = epoch + 1
        if epoch >= epochs:
            print('迭代达到最大步数，自动停止')
            return T_curr
        T_last = T_curr

方差估计（scale estimate）

首先给出 beweight 估计的修正系数 $s_{bi}$ 与样本容量 n 的表格如下（来源 16269-4 的表格 D.1）：

步骤

根据 beweight scale estimate 算法，只需要一步到位计算即可：
$$S_{bi}=s_{bi}\times \frac{n}{\sqrt{n-1}}\frac{\sqrt{{\sum }_{|u_i|<1}(x_i-M{)}^2(1-u_i^2{)}^2}}{|{\sum }_{|u_i|<1}(1-u_i^2)(1-5u_i^2)|}$$
其中 M 为样本中位数，$u_i$ 的计算公式同上。此时常熟 $c$ 一般取 9，意味着样本来源为正态分布时，那些落在 $M\pm 6\times \sigma$ 的点，在计算 $S_{bi}$ 时的权重为 0 。

Python 代码

将 $s_{bi}$ 的系数表格存储在一个与 Python 代码在同一文件夹中的 excel 表格，取名为 correction_factors.xlsx：

def perform_beweight_scale_estimate(x, c=9):
    if isinstance(x, pd.Series):
        x = x.astype('float').values
    
    x_med = np.median(x)
    mad = np.median(np.abs(x-x_med))
    u_arr = (x-x_med)/(c*mad)
    
    # 样本数量
    n = len(x)
    path = r'./correction_factors.xlsx'
    factor_data = pd.read_excel(path, header=0)
    n_list = factor_data['n']
    if n not in n_list:
        if n <= 150:
            n_list_interpolate = np.arange(2, 151)
            sbi_list = factor_data['sbi'].values
            f = interpolate.interp1d(n_list, sbi_list, kind='cubic')
            sbi_list = f(n_list_interpolate)
            sbi = sbi_list[n-2]
        elif 150 < n and n <= 200:
            sbi = 0.9914
        elif 200 < n and n<=300:
            sbi = 0.9912
        else:
            sbi = 0.9910
    else:
        sbi = factor_data.loc[factor_data['n']==n]['sbi']
    
    sbi = float(sbi)
    idx = np.abs(u_arr) < 1
    numerator = np.sqrt(sum((x[idx]-x_med)**2*(1-u_arr[idx]**2)**4))
    denominator = np.abs(A15 算法估计均值方差——原理和 Python 实现

 A15 算法估计均值方差——原理和 Python 实现
 Trimmed 稳健均值估计与 中位数-中位数配对偏差法估计标准差——理论与 Python 实现
 Trimmed 稳健均值估计与 中位数-中位数配对偏差法估计标准差——理论与 Python 实现
 学机器学习要学一些什么？机器学习和深度学习项目实战分享
 Hampel 估计稳健均值——理论和Python实现