M估计（稳健估计）的matlab实现

在进行样本点线性拟合时，一个偏离严重的点（这种点通常称为离群点，又称局外点、外点、出格点）产生的误差相比其他多数正常点占主导地位，这些误差会导致拟合过程中产生巨大的偏差。产生这种后果的原因就是误差被平方了。一个很自然的想法就是给那些离群点更小的权重。常用的一个权重函数为 $w(u)=1/(\sigma^{2}+u^{2})$ ，其中 $\sigma$ 为尺度参数，调节权重衰减的快慢， $u$ 为残差。显而易见，该权重函数给残差大的点分配更小的权重，当残差趋于无穷大，权重也趋于0，符合直觉想法。M估计法就是基于这种权重分配思想设计的，其其完整处理流程如fig1所示(参见《计算机视觉——一种现代方法》第二版p216)：

下面给出本人编写的M估计的matlab程序：

function [a,b,c]=M_estimation(x,y,K,R,err)
% M估计法拟合直线，这里拟合直线全部使用完全最小二乘法，即用ax+by+c=0(a^2+b^2=1)拟合
% 输入k为随机抽取样本的次数
% R为单次采样的采样点个数
% err为终止迭代的系数残差模值
% M估计法权重w=1/(sigma^2+u^2)，u为残差
n=length(x);
x=reshape(x,n,1);
y=reshape(y,n,1);
totalErr=zeros(K,1);
abcM=zeros(K,3);
for i=1:Kindex=randperm(n);indexK=index(1:R);xK=x(indexK);yK=y(indexK);w=ones(R,1)/R;[a,b,c,sigma,~]=weightedLS(w,xK,yK);abc_old=[a,b,c]+err*ones(1,3);abc=[a,b,c];errM=a*xK+b*yK+c;w=1./(sigma^2+errM.*errM);w=w/sum(w);while norm(abc-abc_old)>errabc_old=abc;[a,b,c,sigma,thisErr]=weightedLS(w,xK,yK);abc=[a,b,c];w=1./(sigma^2+errM.*errM);w=w/sum(w);endtotalErr(i)=thisErr;abcM(i,:)=abc;
end
index=find(totalErr==median(totalErr));
abc=abcM(index(1),:);
a=abc(1);
b=abc(2);
c=abc(3);

function [a,b,c,sigma,totalErr]=weightedLS(w,x,y)
% 加权最小二乘法
% w为权重向量
% 用ax+by+c=0(a^2+b^2=1)拟合
% sigma为用于M估计法的尺度参数（权重w=sigma^2/(sigma^2+u^2)，u为残差）
% totalErr返回总的加权残差平方和
n=length(x);
a11=n*sum(w.*x.^2)-power(sum(w.*x),2);
a12=n*sum(w.*x.*y)-sum(w.*x)*sum(w.*y);
a22=n*sum(w.*y.^2)-power(sum(w.*y),2);
A=[a11 a12;a12 a22];
[V,D]=eig(A);
Vab1=V(:,1);
Vab2=V(:,2);
%比较两个解的总误差
avgX=sum(x)/n;
avgY=sum(y)/n;
c1=[-avgX,-avgY]*Vab1;
c2=[-avgX,-avgY]*Vab2;
errM1=Vab1(1)*x+Vab1(2)*y+c1;
errM2=Vab2(1)*x+Vab2(2)*y+c2;
totalErr1=sum(w.*errM1.*errM1);
totalErr2=sum(w.*errM2.*errM2);
if totalErr1<totalErr2a=Vab1(1);b=Vab1(2);c=c1;sigma=1.4826*median(abs(errM1));totalErr=totalErr1;
elsea=Vab2(1);b=Vab2(2);c=c2;sigma=1.4826*median(abs(errM2));totalErr=totalErr2;
end

以下是测试脚本：

clear all;close all;clc;
x=reshape(1:10,10,1);
y=5*x+rand(10,1);
y(9)=y(9)+50;
% w=ones(10,1);
% w(9)=0.1;
p=polyfit(x,y,1);
y1=polyval(p,x);
% [a,b,c,sigma,~]=weightedLS(w,x,y);
[a,b,c]=M_estimation(x,y,9,7,10^-3);
y2=(-c-a*x)/b;
figure,plot(x,y1,'g-',x,y2,'r-',x,y,'*');
legend('最小二乘法','M估计法');

由于M估计法中的随机采样过程，程序运行结果并不固定，当采样点数R较小时也可能发生拟合不太好的情况。但是通过增大采样次数K和每次采样点数R可以大幅提升获得良好拟合结果的概率。

以下是随机抽取的两次运行结果：

可以看出，第一次运行结果显著好于原始最小二乘法结果，第二次没有表现出特别的优势，这是因为样本点数量和单次采样数量都太小了，结果受到随机采样的影响较大，但总的来说受异常样本点还是更小一些。当样本个数和单次采样个数较高的时候，程序可以获得更好的结果。

M估计（稳健估计）的matlab实现相关推荐

【特征提取】基于matlab共振峰估计【含Matlab源码 550期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[特征提取]基于matlab共振峰估计[含Matlab源码 550期] 点击上面蓝色字体,直接付费下载,即可. 获取代码方式2: 付费专栏 ...
【OMP信道估计】基于OMP压缩感知的信道估计算法的MATLAB仿真
1.软件版本 MATLAB2021a 2.本算法理论知识 3.核心代码 clc; clear; close all; warning off; addpath 'func\'CYC = 20; for ...
基于倒谱法、自相关法、短时幅度差法的基音频率估计算法（MATLAB及验证）
基音频率检测一.概念何为基音周期?人在发音时,根据声带是否振动可以将语音信号分为清音和浊音两种.浊音携带大量的能量,因此又被称为有声语音,其在时域上有明显的周期性.而清音类似于白噪声,没有明显的周 ...
【盲信道估计】基于matlab的LMS盲信道估计QPSK仿真
目录 1.软件版本 2.核心代码 3.操作步骤与仿真结论 4.参考文献 5.完整源码获得方式 1.软件版本 MATLAB2021a 2.核心代码 % CHANNEL EQUALIZATION USIN ...
从贝叶斯滤波理论到容积卡尔曼滤波算法（CKF）详细推导及编程实现常转弯率模型估计。（matlab）
容积卡尔曼滤波(CKF)是由加拿大学者Arasaratnam和Haykin在2009年提出的.该算法的核心思想是针对非线性高斯系统,通过三阶球面径向容积准则来近似状态的后验均值和协方差,以保证在理论上 ...
【信号处理】CFO估计技术（Matlab代码实现）
目录 1 概述 2 运行结果 3 参考文献 4 Matlab代码实现 1 概述针对存在未知载波频偏(CFO)的线性调制分类,提出一种混合似似然比检测(qHLRT)分类器.并且通过使用多天线接收机来进 ...
基于扩展卡尔曼滤波的SOC估计（附MATLAB代码）
1.卡尔曼滤波原理原理可以参考我之前学习的笔记,使用goodnote完成的. 我认为,对于公式的推导不需要做太多深入的了解,我之前也对公式进行推导的理解,但是没过几天就忘了,只需要掌握住那重要的5个 ...
随机信号功率谱密度函数理论、估计方法及MATLAB代码
文章的内容整理自网络,仅Matlab代码部分进行了部分修正,具体而言: 理论部分来自:现代通信原理2.5:确定信号的能量谱密度.功率谱密度与自相关函数估计和代码部分来自: 随机信号功率谱密度估计 P ...
【滤波估计】基于matlab双卡尔曼滤波SOC和SOH联合估计【含Matlab源码 2335期】
⛄一.双卡尔曼滤波SOC和SOH联合估计 1 引言为实现节能降耗,降低污染,发展节能环保.不依赖化石燃料的电动汽车取代传统燃油车,已成为当今世界汽车行业的重点发展方向.锂离子动力电池准确可靠的状态估 ...
【压缩感知】基于matlab压缩感知理论的窄带信号DOA估计【含Matlab源码 2616期】
⛄一.压缩感知理论阵列信号波达方向(Direction ofArrival,DOA)估计是阵列信号处理领域中主要研究内容之一,广泛应用于军事及民用领域.基于压缩感知理论的稀疏重构算法的阵列信号DOA ...

M估计（稳健估计）的matlab实现

M估计（稳健估计）的matlab实现相关推荐

最新文章

热门文章