基于EM的多直线拟合实现及思考

作者：桂。

时间:2017-03-22  06:13:50

链接：http://www.cnblogs.com/xingshansi/p/6597796.html

声明：欢迎被转载，不过记得注明出处哦~


前言

分布拟合与曲线拟合系列本想简单梳理，却啰嗦的没完没了。本文主要介绍：多直线的拟合，多曲线可以依次类推。全文主要包括：

　　1）背景介绍

　　2）理论推导

　　3）代码实现

　　4）关于拟合的思考

内容多有借鉴他人，最后一并附上链接。

一、背景介绍

对于单个直线，可以借助MLE或者最小二乘进行求参，对于多条直线呢？

假设一堆数据点（$x_j,l_j$），它由两个线性模型产生：

其中$n_{1j}、n_{2j}$分别为对应的随机噪声。

在分析最小二乘与最大似然联系的时候，知道二者可相互转化；另外在分析混合模型（GMM,LMM）时，都是借助最大似然函数。同样，多直线拟合问题是含有隐变量的最小二乘拟合，也就可以转化为最大似然问题，故求解与混合模型（GMM,LMM）方法类似。

二、理论推导

假设误差服从高斯分布，故可借助GMM来解决该问题（误差服从拉普拉斯分布，则借助LMM来解决）。

 　　A-E-Step

1）求解隐变量，转化为完全数据

${{Z_j} \in {\Upsilon _k}}$表示第$j$个观测点来自第$k$个分模型。

2）构造Q函数

$Q\left( {\Theta ,{\Theta ^{\left( i \right)}}} \right) = \sum\limits_{j = 1}^N {\sum\limits_{k = 1}^K {\log \left( {{w_k}} \right)P\left( {{Z_j} \in {\Upsilon _k}|{Y_j},{\Theta ^{\left( i \right)}}} \right)} }  + \sum\limits_{j = 1}^N {\sum\limits_{k = 1}^K {\log \left( {{f_k}\left( {{Y_j}|{Z_j} \in {\Upsilon _k},{\theta _k}} \right)} \right)} } P\left( {{Z_j} \in {\Upsilon _k}|{Y_j},{\Theta ^{\left( i \right)}}} \right)$

其中${{\theta _k}} = [\mu_k,\sigma_k，a_k, b_k]$为分布$k$对应的参数，$\Theta$  = {$\theta _1$,$\theta _2$,...,$\theta _K$}为参数集合，$N$为样本个数，$K$为混合模型个数。

得到$Q$之后，即可针对完全数据进行MLE求参，可以看到每一个分布的概率（即权重w）与该分布的参数在求参时，可分别求解。由于表达式为一般形式，故该性质对所有混合分布模型都适用。所以对于混合模型，套用Q并代入分布具体表达式即可。

　　B-M-Step

1）利用MLE求参

  • 首先对${{w_k}}$进行优化

由于$\sum\limits_{k = 1}^M {{w_k}}  = 1$，利用Lagrange乘子求解：

${J_w} = \sum\limits_{j = 1}^N {\sum\limits_{k = 1}^K {\left[ {\log \left( {{w_k}} \right)P\left( {\left. {{Z_j} \in {\Upsilon _k}} \right|{Y_j},{{\bf{\Theta }}^{\left( i \right)}}} \right)} \right]} }  + \lambda \left[ {\sum\limits_{k = 1}^K {{w_k}}  - 1} \right]$

求偏导：

$\frac{{\partial {J_w}}}{{\partial {w_k}}} = \sum\limits_{J = 1}^N {\left[ {\frac{1}{{{w_k}}}P\left( {{Z_j} \in {\Upsilon _k}|{Y_j},{{\bf{\Theta }}^{\left( i \right)}}} \right)} \right] + } \lambda  = 0$

得

  • 对各分布内部参数$\theta_k$进行优化

给出准则函数：

${J_\Theta } = \sum\limits_{j = 1}^N {\sum\limits_{k = 1}^K {\log \left( {{f_k}\left( {{Y_j}|{Z_j} \in {\Upsilon _k},{\theta _k}} \right)} \right)} } P\left( {{Z_j} \in {\Upsilon _k}|{Y_j},{\Theta ^{\left( i \right)}}} \right)$

对于多直线拟合问题，$Y_j$为拟合残差，假设其服从高斯分布：

可以认为${{l_j} - {a_k}{x_j}}$就是GMM中的$Y_j$，$b_k$就是$\mu_k$。直接套用GMM中的迭代结果：

所不同的是，多了一个对$a_k$的求解，容易得出：

至此，理论推导完成。

三、代码实现

仍然是在之前GMM代码基础上，修改几句指令：

function [u,sig,a,t,iter] = fit_mix_line( X,l,M )
%
% fit_mix_line - fit parameters for a mixed-line using EM algorithm
%
% format:   [u,sig,t,iter] = fit_mix_line( X,M )
%
% input:    X   - input samples, Nx1 vector
%           M   - number of gaussians which are assumed to compose the distribution
%
% output:   u   - fitted mean for each gaussian
%           sig - fitted standard deviation for each gaussian
%           t   - probability of each gaussian in the complete distribution
%           iter- number of iterations done by the function
%

% initialize and initial guesses
N           = length( X );
Z           = ones(N,M) * 1/M;                  % indicators vector
P           = zeros(N,M);                       % probabilities vector for each sample and each model
t           = ones(1,M) * 1/M;                  % distribution of the gaussian models in the samples
u           = linspace(min(X),max(X),M);        % mean vector
sig2        = ones(1,M) * var(X) / sqrt(M);     % variance vector
C           = 1/sqrt(2*pi);                     % just a constant
Ic          = ones(N,1);                        % - enable a row replication by the * operator
Ir          = ones(1,M);                        % - enable a column replication by the * operator
a          = ones(1,M);
Q           = zeros(N,M);                       % user variable to determine when we have converged to a steady solution
thresh      = 1e-9;
step        = N;
last_step   = 10;         % step/last_step
iter        = 0;
min_iter    = 3000;   

% main convergence loop, assume gaussians are 1D
while ((( abs((step/last_step)-1) > thresh) & (step>(N/5*eps)) ) & (iter<min_iter) ) 

    % E step
    % ========
    Q   = Z;
    P   = C ./ (Ic*sqrt(sig2)) .* exp( -(((l*Ir-X*a) - Ic*u).^2)./(2*Ic*sig2) );

    for m = 1:M
        Z(:,m)  = (P(:,m)*t(m))./(P*t(:));
    end

    % estimate convergence step size and update iteration number
    prog_text   = sprintf(repmat( ‘\b‘,1,(iter>0)*12+ceil(log10(iter+1)) ));
    iter        = iter + 1;
    last_step   = step * (1 + eps) + eps;
    step        = sum(sum(abs(Q-Z)));
    fprintf( ‘%s%d iterations\n‘,prog_text,iter );

    % M step
    % ========
    Zm              = sum(Z);               % sum each column
    Zm(find(Zm==0)) = eps;                  % avoid devision by zero
    sig2            = sum((((l*Ir-X*a) - Ic*u).^2).*Z) ./ Zm;
    u               = sum((l*Ir-X*a).*Z) ./ Zm;
    a               = sum((l*Ir - Ic*u).*(X*Ir).*Z) ./ (sum((X*Ir).^2.*Z));
%     a (isnan(a))    = 0.001;
    t               = Zm/N;
end
sig     = sqrt( sig2 );

给出测试程序：

clc;clear all;close all
set(0,‘defaultfigurecolor‘,‘w‘)
%generate data
x = linspace(-40,40,200);
y = zeros(1,length(x));
y1 = zeros(1,length(x)/2);
y2 = zeros(1,length(x)/2);
k1= 0;k2=0;
for i =1 :length(x)
    if mod(i,2)==0
        k1=k1+1;
        y(i) = 5*x(i)-3 + 3*rand;%分别取0.5 和5
        y1(k1)=y(i);
    else
        k2=k2+1;
        y(i)= -7*x(i)+2 + 3*rand;
        y2(k2)=y(i);
    end
end
[u,sig,a] = fit_mix_line(x‘,y‘,2);
yo=[y1,y2];
[uo,sigo,ao] = fit_mix_line(x‘,yo‘,2);
%figure
subplot 211
scatter(x,y,‘k.‘);
hold on;
t = -20:20;
l1 = t*a(1)+u(1);
l2 = t*a(2)+u(2);
plot(t,l1,‘r‘,‘linewidth‘,2);hold on;
plot(t,l2,‘g--‘,‘linewidth‘,2);hold on;
grid on;
subplot 212
scatter(x,yo,‘k.‘);
hold on;
l1 = t*ao(1)+uo(1);
l2 = t*ao(2)+uo(2);
plot(t,l1,‘r‘,‘linewidth‘,2);hold on;
plot(t,l2,‘g--‘,‘linewidth‘,2);hold on;
grid on;

这里分别针对两种多线性进行拟合：

  • 分段多条直线
  • 混合多条直线

理论上二者都适用，但运行却发现二者往往只有一个理想，记录此处，暂时未找出原因。

代码中 y(i) = 5*x(i)-3 + 3*rand;%分别取0.5 和5这一句取0.5时，结果图：

取5时，对应结果图：

理论上应该二者都适用。

四、关于拟合的思考

　　A-以正态分布为例

上面分析的多直线拟合，其实是$ax+b$的形式，由此构造混合分布，对于：

更一般的：

$g$为一般表达式，（如GMM就是$g = ax+b$，且a=0的情况,上文分析的为a不等于0的情况），更一般的$g$理论上可以为任意表达式：

只要将g的具体表达式代入EM求解过程即可。

　　B-其他分布

上文的讨论基于噪声是正态分布，如果是拉普拉斯分布呢？只要将上面更一般表达式提到的外壳换成拉普拉斯分布模型即可。

事实上，EM的混合模型到此可以看出：混合模型理论上可以实现各类形状的聚类，而噪声同样可以基于不同的分布假设。

参考：

李航：《统计学习方法》

时间： 2024-10-20 12:13:33

基于EM的多直线拟合实现及思考的相关文章

OpenCV fitline直线拟合函数学习

下图是OpenCV官方文档中,对直线拟合函数的详细介绍: fitLine()函数用于,对二维或三维空间中的点集进行直线拟合.共有六个参数: param 1:输入的点集,可以是Mat或者vector<>,可以是二维点集或三维点集. 例如: vector<Point> points; param 2:拟合结果,即一条直线.在二维空间中,直线可以定义为 Vec4f line; 在二维平面中,(line[0],line[1])表示直线的方向向量,(line[2],line[3])表示直线上

直线拟合算法(续)

直线拟合算法(续) 曾经写过一篇博客.介绍直线拟合算法. http://blog.csdn.net/liyuanbhu/article/details/50866802 给出的代码事实上有一点小问题,就是 den = 0 时会出现除以 0 的错误. 今天正好也有网友问起这个问题. 我就再写一篇短文来说说怎样解决问题. 首先我们知道: den=D2xy+(λ?Dxx)2??????????????√ 那么 den=0 意味着: Dxy=0λ=Dxx 我们还有关于 λ 的计算式: λ=Dxx+Dyy

直线拟合算法

在计算机视觉的应用中,经常会用到提取一条直线的精确位置这样的工作.这时就要用到直线的拟合算法了. 这里,我也贴一个利用最小二乘法计算最佳拟合直线的代码. 这个代码是我以前学习<机器视觉算法与应用(双语版)>[德] 斯蒂格(Steger C) 著:杨少荣 等 译 的书时写的.所有的公式推导都在书中 3.8.1 ,还算比较有用. 与一元线性回归算法的区别:一元线性回归算法假定 X 是无误差的,只有 Y 有误差. 而这个算法假设每个点的 X Y 坐标的误差都是符合 0 均值的正态分布的. 因此,在计

OpenCV 学习（直线拟合)

OpenCV 学习(直线拟合) Hough 变换可以提取图像中的直线.但是提取的直线的精度不高.而很多场合下,我们需要精确的估计直线的参数,这时就需要进行直线拟合. 直线拟合的方法很多,比如一元线性回归就是一种最简单的直线拟合方法.但是这种方法不适合用于提取图像中的直线.因为这种算法假设每个数据点的X 坐标是准确的,Y 坐标是带有高斯噪声的.可实际上,图像中的每个数据点的XY 坐标都是带有噪声的. 下面就来讲讲适用于提取图像中直线的直线拟合算法. 一个点 (xi,yi) 到直线的距离用 ri 来

最小二乘法直线拟合

最小二乘法的直线拟合 #coding:utf-8 import numpy as np import matplotlib.pyplot as plt dots = np.array([[1,6], [2,5], [3,7], [4,10]]) plt.plot([i[0] for i in dots],      [i[1] for i in dots], 'ro') plt.axis([0, 6, 0, 12]) def nihezhixian(k, x, b):     return k*

ax+by+c=0 型直线拟合算法

所谓直线拟合,通常也叫做线性拟合.一元线性回归.指的是当我们有一批数据(xi,yi),这些数据在平面坐标系下落在一条直线上,或近似的落在一条直线上.我们就要求出这条直线的参数.如果这条直线可以写为: y=kx+b 那么 k=∑(xi?xˉ)(yi?yˉ)∑(xi?xˉ)2 b=yˉ?kxˉ 这个关系式许多教科书上都有详细的推导,无需多说. 今天要说的是另一种情况,当我们的数据有可能落在一条竖直的直线上,也就是k 有可能为∞ 时,应该如何做拟合.这时我们肯定就不能用y=kx+b 了,但是可以将这个

（原）opencv直线拟合fitLine

转载请注明出处 http://www.cnblogs.com/darkknightzh/p/5486234.html 参考网址: http://blog.csdn.net/thefutureisour/article/details/7599537 1 Mat img(60, 60, CV_8UC1, Scalar(0)); 2 std::vector<Point2f> points; 3 points.push_back(Point2f(10.5, 10.2)); 4 points.push

基于MATLAB的多项式数据拟合方法研究-毕业论文

摘要:本论文先介绍了多项式数据拟合的相关背景,以及对整个课题做了一个完整的认识.接下来对拟合模型,多项式数学原理进行了详细的讲解,通过对文献的阅读以及自己的知识积累对原理有了一个系统的认识.介绍多项式曲线拟合的基本理论,对多项式数据拟合原理进行了全方面的理论阐述,同时也阐述了曲线拟合的基本原理及多项式曲线拟合模型的建立.具体记录了多项式曲线拟合的具体步骤,在建立理论的基础上具体实现多项式曲线的MATLAB实现方法的研究,采用MATLAB R2016a的平台对测量的数据进行多项式数据拟合,介绍了M

opencv直线拟合cv::fitLine()

通过2D或者2D点集拟合直线 void fitLine( InputArray points, OutputArray line, int distType,double param, double reps, double aeps ); @param points Input vector of 2D or 3D points, stored in std::vector\<\> or Mat.@param line Output line parameters. In case of 2