基于EM（Expectation Maximization，最大期望）算法的直线拟合matlab实现

EM算法的主要思想就是给出一个初始分组信息，迭代地基于期望最大化进行参数的估计和分组的估计。用在直线拟合上，就是先预设某些点在直线上，然后基于期望最大化估计残差均方根（即估计数据波动范围），然后再估计哪些点属于正常数据，用这些正常数据再进行拟合，然后再估计残差均方根。这个过程不断进行直到收敛。下面给出个人实现的基于EM算法的直线拟合程序：

function [a,b,c]=EM_linefit(x,y,k)
% 用EM方法进行直线拟合
% k为随机初始采样点个数
n=length(x);
index=randperm(n);
chosed=zeros(1,n);
chosed(index(1:k))=1;
chosed=logical(chosed);
flag=true;
while flag
    chosed_old=chosed;
    [a,b,c,sigma]=completeLS(x(chosed),y(chosed));
    errM=abs(a*x+b*y+c);
    chosed=errM<3*sigma; %落在3*sigma波动范围内，认为属于正常数据，用于下一次拟合
    if sum(chosed~=chosed_old)==0 %下次用于拟合的点索引不变，终止迭代
        flag=false;
    end
end
end

function [a,b,c,sigma]=completeLS(x,y)
% 完全最小二乘法
% 用ax+by+c=0(a^2+b^2=1)拟合
% sigma为根据最大似然法估计出的残差均方根
n=length(x);
a11=n*sum(x.^2)-power(sum(x),2);
a12=n*sum(x.*y)-sum(x)*sum(y);
a22=n*sum(y.^2)-power(sum(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(errM1.*errM1);
totalErr2=sum(errM2.*errM2);
if totalErr1<totalErr2
    a=Vab1(1);
    b=Vab1(2);
    c=c1;
    sigma=sqrt(totalErr1/n);
else
    a=Vab2(1);
    b=Vab2(2);
    c=c2;
    sigma=sqrt(totalErr2/n);
end

脚本函数如下：

clear all;close all;clc;
x=reshape(1:10,10,1);
y=5*x+rand(10,1);
y(9)=y(9)+50;
p=polyfit(x,y,1);
y1=polyval(p,x);
[a,b,c]=EM_linefit(x,y,2);
y2=(-c-a*x)/b;
figure,plot(x,y1,'g-',x,y2,'r-',x,y,'*');
legend('最小二乘法','EM法');

随机执行结果如下：

fig1 某次算法执行结果

可以看出，该次执行中程序识别出了异常点并只用正常点进行了拟合，拟合结果优于最小二乘法。不过程序运行结果具有随机性，如果初始选择的样本中含有异常点，有可能会陷入错误拟合无法跳出。可以考虑这样的解决方案：如果选中的都是正常点，那么丢点一个点重新拟合，则拟合误差应该不会有太大的变化，如此便终止迭代。如果丢掉一个点后残差显著变大，则将残差模值最大的数据剔除继续进行迭代。由于这里只是简单介绍下EM算法的应用做个抛砖引玉，所以更完善的解决方案就留给大家了。