作者：桂。

時(shí)間:2017-03-22 ?06:13:50

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

聲明：歡迎被轉(zhuǎn)載，不過記得注明出處哦~

---

?

前言

分布擬合與曲線擬合系列本想簡單梳理，卻啰嗦的沒完沒了。本文主要介紹：多直線的擬合，多曲線可以依次類推。全文主要包括：

　　1）背景介紹

　　2）理論推導(dǎo)

　　3）代碼實(shí)現(xiàn)

　　4）關(guān)于擬合的思考

內(nèi)容多有借鑒他人，最后一并附上鏈接。

?

一、背景介紹

對于單個(gè)直線，可以借助MLE或者最小二乘進(jìn)行求參，對于多條直線呢？

假設(shè)一堆數(shù)據(jù)點(diǎn)（$x_j,l_j$），它由兩個(gè)線性模型產(chǎn)生：

其中$n_{1j}、n_{2j}$分別為對應(yīng)的隨機(jī)噪聲。

在分析最小二乘與最大似然聯(lián)系的時(shí)候，知道二者可相互轉(zhuǎn)化；另外在分析混合模型（GMM,LMM）時(shí)，都是借助最大似然函數(shù)。同樣，多直線擬合問題是含有隱變量的最小二乘擬合，也就可以轉(zhuǎn)化為最大似然問題，故求解與混合模型（GMM,LMM）方法類似。

?

二、理論推導(dǎo)

假設(shè)誤差服從高斯分布，故可借助GMM來解決該問題（誤差服從拉普拉斯分布，則借助LMM來解決）。

?　　A-E-Step

1）求解隱變量，轉(zhuǎn)化為完全數(shù)據(jù)

${{Z_j} \in {\Upsilon _k}}$表示第$j$個(gè)觀測點(diǎn)來自第$k$個(gè)分模型。

2）構(gòu)造Q函數(shù)

$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$對應(yīng)的參數(shù)，$\Theta$ ?= {$\theta _1$,$\theta _2$,...,$\theta _K$}為參數(shù)集合，$N$為樣本個(gè)數(shù)，$K$為混合模型個(gè)數(shù)。

得到$Q$之后，即可針對完全數(shù)據(jù)進(jìn)行MLE求參，可以看到每一個(gè)分布的概率（即權(quán)重w）與該分布的參數(shù)在求參時(shí)，可分別求解。由于表達(dá)式為一般形式，故該性質(zhì)對所有混合分布模型都適用。所以對于混合模型，套用Q并代入分布具體表達(dá)式即可。

　　B-M-Step

1）利用MLE求參

• 首先對${{w_k}}$進(jìn)行優(yōu)化

由于$\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]$

求偏導(dǎo)：

$\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$

?得

• 對各分布內(nèi)部參數(shù)$\theta_k$進(jìn)行優(yōu)化

給出準(zhǔn)則函數(shù)：

${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$為擬合殘差，假設(shè)其服從高斯分布：

可以認(rèn)為${{l_j} - {a_k}{x_j}}$就是GMM中的$Y_j$，$b_k$就是$\mu_k$。直接套用GMM中的迭代結(jié)果：

所不同的是，多了一個(gè)對$a_k$的求解，容易得出：

至此，理論推導(dǎo)完成。

?

三、代碼實(shí)現(xiàn)

仍然是在之前GMM代碼基礎(chǔ)上，修改幾句指令：

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:MZ(:,m) = (P(:,m)*t(m))./(P*t(:));end% estimate convergence step size and update iteration numberprog_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 columnZm(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)==0k1=k1+1;y(i) = 5*x(i)-3 + 3*rand;%分別取0.5 和5y1(k1)=y(i);elsek2=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;

這里分別針對兩種多線性進(jìn)行擬合：

• 分段多條直線
• 混合多條直線

理論上二者都適用，但運(yùn)行卻發(fā)現(xiàn)二者往往只有一個(gè)理想，記錄此處，暫時(shí)未找出原因。

代碼中?y(i) = 5*x(i)-3 + 3*rand;%分別取0.5 和5這一句取0.5時(shí)，結(jié)果圖：

取5時(shí)，對應(yīng)結(jié)果圖：

理論上應(yīng)該二者都適用。

?

四、關(guān)于擬合的思考

　　A-以正態(tài)分布為例

上面分析的多直線擬合，其實(shí)是$ax+b$的形式，由此構(gòu)造混合分布，對于：

更一般的：

$g$為一般表達(dá)式，（如GMM就是$g = ax+b$，且a=0的情況,上文分析的為a不等于0的情況），更一般的$g$理論上可以為任意表達(dá)式：

只要將g的具體表達(dá)式代入EM求解過程即可。

　　B-其他分布

上文的討論基于噪聲是正態(tài)分布，如果是拉普拉斯分布呢？只要將上面更一般表達(dá)式提到的外殼換成拉普拉斯分布模型即可。

事實(shí)上，EM的混合模型到此可以看出：混合模型理論上可以實(shí)現(xiàn)各類形狀的聚類，而噪聲同樣可以基于不同的分布假設(shè):

1)Kmeans是對于 ?中心點(diǎn)（聚類中心） 的分布假設(shè)；

2）GMM/LMM是對于 斜率為0的直線（GMM的均值） 的分布假設(shè)；

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

?

參考：

李航：《統(tǒng)計(jì)學(xué)習(xí)方法》

轉(zhuǎn)載于:https://www.cnblogs.com/xingshansi/p/6597796.html

總結(jié)

以上是生活随笔為你收集整理的基于EM的多直线拟合实现及思考的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。