前言

数据探测和稳健估计均可以探测出粗差，但是数据探测仅仅在数据中仅一个观测量含有粗差时好用，如果想要检验多个粗差，需要逐步剔除已检验粗差，再继续检验。而稳健估计在处理包含有多个粗差时有一定优势。

稳健估计

稳健估计分为三大类：M估计，R估计，L估计，M估计最常用，它是一种广义的极大似然估计，有分为选权迭代法和P范数最小法。

M估计的准则

其中ρ为M估计的函数，p为权
M估计的稳健性与ρ的选择有关，当ρ=v**2时，即为最小二乘估计，它并不具有抗差性。
选取不同的ρ函数，会得出不同的M估计，因此M估计不是一个估计，而是一类估计。

选权迭代法

流程：
（1）建立数学模型：V=BX-L
（2）按照最小二乘法求解参数估计值及其残差：X=inv(N)*w，V=BX-L
（3）由V求解观测值的等价权矩阵P_，应用抗差最小二乘（P_代替P）进行迭代计算，直到ΔX<=η，η为迭代停止条件
（4）最后计算出X,V,权函数ρ

几种常用的选权迭代法

Huber，Hampel，丹麦法，IGG，一次范数最小法，P范数最小法，相关等价权法。
其中P范数最小法
ρ函数：

权因子：

python实现

输入数据

import numpy as npnp.set_printoptions(suppress=True)  # 抑制科学记数法
# B,L,P矩阵B = np.array([0, 1, 0, -1, 1, 0, -1, 0, -1, 1]).reshape(5, 2)
L = np.array([0, -1, 0, 4, 8])
P = np.array([2, 2, 2, 2, 1])
P = np.diag(P)

计算X,V

def cal_canshu(B, P, L):N = np.dot(B.T, P).dot(B)W = np.dot(B.T, P).dot(L)inv_N = np.linalg.inv(N)X = inv_N.dot(W)V = B.dot(X) - Lreturn X, V
# 由初始权阵计算X,V
X, V = cal_canshu(B, P, L)
# print(X,V)

计算第一次的w

# 选权迭代：P范最小法m = 1.5
k = 0.0002
yita = 0.000001  # 终止循环条件
yita_array = np.array([yita] * 2)
w = []  # 权因子矩阵
for v in V:w_v = (abs(v) ** (2 - m) + k) ** -1w.append(w_v)P_ = P * w  # 更新权阵# print(P_)

计算更新之后X,V

# 计算更新权阵之后的X,V
X_, V_ = cal_canshu(B, P_, L)

选权迭代

delta = abs(X_ - X)
delta[0] = round(delta[0], 7)
delta[1] = round(delta[1], 7)  # 保留7位小数while all(delta >= yita_array):  # 这里注意如果要比较矩阵大小，需选择all（所有因子）/any（一个及以上）进行比较w_v = []for i in V_:w_a = (abs(i) ** (2 - m) + k) ** (-1)w_v.append(w_a)# print(w_v)X = X_P_ = P_ * w_v  # 更新权阵X_, V_ = cal_canshu(B, P_, L)delta = abs(X_ - X)delta[0] = round(delta[0], 7)delta[1] = round(delta[1], 7)  # 保留7位小数

查看ρ函数

ρ = []
for i in V_:ρ.append([abs(i) ** m])
print(ρ, X_)

结果

[[0.9999999999960162], [4.328242386676566e-18], [7.999999908958161], [5.28663897358714e-12], [5.196152501558078]] [-3.99999997  1.        ]

从步骤可以看出，该方法关键是确定等价权，权函数是一个在平差过程中随改正数变化的量，经过多次迭代，含有粗差的异常观测的权函数接近0
**权函数不是权阵，异常值的权会变大(L2,L4)

稳健估计,P范数最小法相关推荐

l20范数最小化求解系数方程_贪婪组稀疏方法(Greedy group sparsity)
l20范数最小化求解系数方程_贪婪组稀疏方法(Greedy group sparsity) 本文章部分参考Fast group sparse classification l20范数最小化求解系数方程 ...
三维网格去噪算法（L0范数最小化，包含二维图像去噪）
参考文章(技术来源):Mesh denoising via L0 minimization 上面参考文章提出了一种基于L0范数最小化的三角网格去噪算法.该思想由二维图像平滑引申而来,所以先从基于L0范 ...
压缩感知的尽头: 原子范数最小化
文章目录前言问题建模 Toeplitz 矩阵的范德蒙德分解 DOA估计的一般框架 ℓ0\ell_0ℓ0-原子范数 ℓ0\ell_0ℓ0-原子范数与范德蒙德分解原子范数多维原子范数证明 ...
稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现
文章目录一.前言二.问题重述三.构造 ℓ1\ell_1ℓ1 范数四.ℓ1\ell_1ℓ1 范数最小化问题转换为线性规划问题五.基于linprog的基追踪Python代码六.运行测试七 ...
最大最小法及α-β剪枝算法图解
(网上讲的都不是很好理解,贡献一下之前听慕课做的笔记,适合初学者比较简洁明了.) 要想理解α-β剪枝算法,必须从最大最小法的博弈问题讲起!注意不要跳过第一节往下看. 最大最小法场景:双方博弈前提: ...
l1范数最小化快速算法【文献阅读】
1:解决的问题模型如下: 或者约束条件可以适当的松弛,即为如下模型: 当然约束条件取范数,数据获取的比较准确,结果逼近的效果更好,防止过拟合.如果取范数,则是获取的数据,受到污染比较严重. ...
l1范数最小化快速算法
1:解决的问题模型如下: 或者约束条件可以适当的松弛,即为如下模型: 当然约束条件取
图像处理-基于图像梯度L0范数最小化(L0smooth)的保护边缘平滑滤波
算法程序备注: (1)下面是对一幅自然图像进行处理的结果: 可以看到图像有非常明显的变化,图像分成了一块一块,这是图像平滑后的结果,因为保护了边界,因此明显的边界仍然存在,但是不可避免的细节部分被磨平 ...
欠定方程组的最小范数解
欠定方程组的最小范数解(超曲面和原点的距离,原点球心球面和超曲面的切点) 今天小白问了我一个问题,我觉得颇有意思,做个记录. 问题描述小白:插入一个问题,这个切点怎么求?就是固定超平面Ax=b 怎么 ...
原子范数 Atomic norm最小化：简单的Matlab例程
前言基于压缩感知的尽头: 原子范数最小化中的原子范数最小化算法, 笔者做了一些matlab的仿真, 作为简单的例程,希望帮助大家进一步理解算法和自定义的拓展. 由于凸问题的求解需要使用 CVX, ...

稳健估计,P范数最小法

前言