混合拉普拉斯

作者：桂。

时间：2017-03-21  07:25:17

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

前言

本文为曲线拟合与分布拟合系列的一部分，主要讲解混合拉普拉斯分布（Laplace Mixture Model，LMM）。拉普拉斯也是常用的统计概率模型之一，网上关于混合高斯模型（GMM）的例子很多，而关于LMM实现的很少。其实混合模型都可以用EM算法推导，只是求闭式解的运算上略有差别，全文包括：

　　1）LMM理论推导；

　　2）LMM代码实现；

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

 

一、LMM理论推导

　　A-模型介绍

对于单个拉普拉斯分布，表达式为：

f(Y)=12be−|Y−μ|b

对于 K个模型的混合分布：

P(Yj|Θ)=∑k=1Kwkf(Yj|μk,bk)

如何拟合呢？下面利用EM分析迭代公式，仅分析Y为一维的情况，其他可类推。（先给出一个结果图）

　　B-EM算法推导

E-Step:

1)求解隐变量，构造完全数据集

同GMM推导类似，利用全概率公式：

2）构造Q函数

基于之前混合高斯模型（GMM）的讨论，EM算法下混合模型的Q函数可以表示为：

Q(Θ,Θ(i))=∑j=1N∑k=1Klog⁡(wk)P(Zj∈Υk|Yj,Θ(i))+∑j=1N∑k=1Klog⁡(fk(Yj|Zj∈Υk,θk))P(Zj∈Υk|Yj,Θ(i))

其中 θk=[μk,bk]为分布 k对应的参数， Θ  = { θ1, θ2,..., θK}为参数集合， N为样本个数， K为混合模型个数。

M-Step:

1）MLE求参

• 首先对 wk进行优化

由于 ∑k=1Mwk=1，利用Lagrange乘子求解：

Jw=∑j=1N∑k=1K[log⁡(wk)P(Zj∈Υk|Yj,Θ(i))]+λ[∑k=1Kwk−1]

求偏导：

∂Jw∂wk=∑J=1N[1wkP(Zj∈Υk|Yj,Θ(i))]+λ=0

 得

• 对各分布内部参数 θk进行优化

给出准则函数：

JΘ=∑j=1N∑k=1Klog⁡(fk(Yj|Zj∈Υk,θk))P(Zj∈Υk|Yj,Θ(i))

仅讨论 Yj为一维数据情况，其他类推。对于拉普拉斯分布：

关于 θk利用MLE即可求参。

首先求解 bk的迭代公式：

由于 μk含有绝对值，因此需要一点小技巧。 JΘ对 μk求偏导，得到：

得到的 μk估计即为：

μk(i+1)=μ^k

在迭代的最终状态，可以认为 i次参数与 i+1次参数近似相等，从而上面的求导结果转化为：

得到参数 μk的迭代公式：

总结一下LMM的求解步骤：

E-Step:

M-Step:

 

二、LMM代码实现

 根据上一篇GMM的代码，简单改几行code，即可得到LMM：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

function  [u,b,t,iter] = fit_mix_laplace( X,M ) % % fit_mix_laplace - fit parameters for a mixed-laplacian distribution using EM algorithm % % format:   [u,b,t,iter] = fit_mix_laplacian( 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 laplacian %           b - fitted standard deviation for each laplacian %           t   - probability of each laplacian in the complete distribution %           iter- number of iterations done by the function % 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 (0.2,1.4,M);         % mean vector b           =  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 Q           =  zeros (N,M);                        % user variable to determine when we have converged to a steady solution thresh      = 1e-7;         step        = N; last_step   = 300;          % step/last_step iter        = 0; min_iter    = 3000;         while  (((  abs ((step/last_step)-1) > thresh) & (step>(N*1e-10)) ) & (iter<min_iter) )      % E step      % ========      Q   = Z;      P   = 1./ (Ic*b) .*  exp ( -(1e-6+ abs (X*Ir - Ic*u))./(Ic*b) );      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      u               =  sum (((X*Ir)./ abs (X*Ir - Ic*u)).*Z) ./ sum (1./ abs (X*Ir - Ic*u).*Z) ;      b               =  sum (( abs (X*Ir - Ic*u)).*Z) ./ Zm ;      t               = Zm/N; end end

给出上文统计分布的拟合程序：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

clc ; clear  all ; %generate random xmin = -10; xmax = 10; Len = 10000000; x =  linspace (xmin,xmax,Len); mu = [3,-4]; b = [0.9 0.4]; w = [0.7 0.3]; fx = w(1)/2/b(1)* exp (- abs (x-mu(1))/b(1))+ w(2)/2/b(2)* exp (- abs (x-mu(2))/b(2)); ymax = 1/b(2); ymin = 0; Y = (ymax-ymin)* rand (1,Len)-ymin; data = x(Y<=fx); %Laplace Mixture Model fitting K = 2; [mu_new,b_new,w_new,iter] = fit_mix_laplace( data',K); %figure subplot  221 hist (data,2000); grid  on; subplot  222 numter = [xmin:.2:xmax]; plot (numter,w_new(1)/2/b_new(1)* exp (- abs (numter-mu_new(1))/b_new(1)), 'r' , 'linewidth' ,2); hold  on; plot (numter,w_new(2)/2/b_new(2)* exp (- abs (numter-mu_new(2))/b_new(2)), 'g' , 'linewidth' ,2); hold  on;   subplot  (2,2,[3,4]) [histFreq, histXout] =  hist (data, numter); binWidth = histXout(2)-histXout(1); %Bar bar (histXout, histFreq/binWidth/ sum (histFreq));  hold  on; grid  on; plot (numter,w_new(1)/2/b_new(1)* exp (- abs (numter-mu_new(1))/b_new(1)), 'r' , 'linewidth' ,2); hold  on; plot (numter,w_new(2)/2/b_new(2)* exp (- abs (numter-mu_new(2))/b_new(2)), 'g' , 'linewidth' ,2); hold  on;

对应结果图（与上文同）：

参考

• Mitianoudis N, Stathaki T. Batch and online underdetermined source separation using Laplacian mixture models[J]. IEEE Transactions on Audio, Speech, and Language Processing, 2007, 15(6): 1818-1832.