作者：桂。

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

链接： 

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

---

 

前言

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

　　1）背景介绍

　　2）理论推导

　　3）代码实现

　　4）关于拟合的思考

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

 

一、背景介绍

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

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

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

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

 

二、理论推导

假设误差服从高斯分布，故可借助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 guessesN           = length( X );Z           = ones(N,M) * 1/M;                  % indicators vectorP           = zeros(N,M);                       % probabilities vector for each sample and each modelt           = ones(1,M) * 1/M;                  % distribution of the gaussian models in the samplesu           = linspace(min(X),max(X),M);        % mean vectorsig2        = ones(1,M) * var(X) / sqrt(M);     % variance vectorC           = 1/sqrt(2*pi);                     % just a constantIc          = ones(N,1);                        % - enable a row replication by the * operatorIr          = ones(1,M);                        % - enable a column replication by the * operatora          = ones(1,M);                       Q           = zeros(N,M);                       % user variable to determine when we have converged to a steady solutionthresh      = 1e-9;        step        = N;last_step   = 10;         % step/last_stepiter        = 0;min_iter    = 3000;   % main convergence loop, assume gaussians are 1Dwhile ((( abs((step/last_step)-1) > thresh) & (step>(N/5*eps)) ) & (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;endsig     = sqrt( sig2 );
     

给出测试程序：

clc;clear all;close allset(0,'defaultfigurecolor','w')%generate datax = 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);    endend[u,sig,a] = fit_mix_line(x',y',2);yo=[y1,y2];[uo,sigo,ao] = fit_mix_line(x',yo',2);%figuresubplot 211scatter(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 212scatter(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的混合模型到此可以看出：混合模型理论上可以实现各类形状的聚类，而噪声同样可以基于不同的分布假设:

1)Kmeans是对于  中心点（聚类中心） 的分布假设；

2）GMM/LMM是对于 斜率为0的直线（GMM的均值） 的分布假设；

3）更一般地，斜线/二次曲线....以此类推。

 

参考：

李航：《统计学习方法》