混合高斯模型(GMM)推导及实现

作者:桂。

时间:2017-03-20  06:20:54

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


前言

本文是曲线拟合与分布拟合系列的一部分,主要总结混合高斯模型(Gaussian Mixture Model,GMM),GMM主要基于EM算法(前文已经推导),本文主要包括:

  1)GMM背景介绍;

  2)GMM理论推导;

  3)GMM代码实现;

内容多有借鉴他人,最后一并给出链接。

 

一、GMM背景

  A-高斯模型1

给出单个随机信号(均值为-2,方差为9的高斯分布),可以利用最大似然估计(MLE)求解分布参数:

  B-高斯模型2

对于单个高斯模型2(均值为3,方差为1),同样可以利用MLE求解:

  C-高斯模型3

现在对于一个随机数,每一个点来自混合模型1概率为0.5,来自混合模型2概率为0.5,得到统计信息:

可能已经观察到:只要将信号分为前后两段分别用MLE解高斯模型不就可以?其实这个时候,已经默默地用了一个性质:数据来自模型1或2的概率为0.5,可见一旦该特性确定,混合模型不过是普通的MLE求解问题,可现实情况怎么会这么规律呢,数据来自模型1或2的概率很难通过观察得出。观测数据$Y_1$来自模型1,$Y_2$来自模型2...参差交错。

再分两段看看?如果直接利用MLE求解,这就碰到了与之前分析EM时:硬币第三抛同样的尴尬。先看一下EM解决的效果:

其实硬币第三抛,也是一个混合概率模型:对于任意一个观测点以概率$\pi$选择硬币A,以概率$1-\pi$选择硬币B,对应混合模型为:

$P\left( {{Y_j}|\theta } \right) = {w_1}{P_A} + {w_2}{P_B} = \pi {P_A} + \left( {1 - \pi } \right){P_B}$

同样,对于两个高斯的混合模型(连续分布,故不用分布率,而是概率密度):

推而广之,对于K个高斯的混合模型:

 

二、GMM理论推导

可以看出GMM与抛硬币完全属于一类问题,故采用EM算法求解,按模式识别(2)——EM算法的思路进行求解。

记:观测数据为$Y$={$Y_1,Y_2,...Y_N$},对应隐变量为$Z$={$Z_1,Z_2,...Z_N$}。

写出EM算法中Q函数的表达式:

E-Step:

1)将缺失数据,转化为完全数据

 主要求解:$P\left( {{Z_j}|{Y_j},{\Theta ^{\left( i \right)}}} \right)$,此处的求解与EM算法一文中硬币第三抛的思路一致,只要求出$P\left( {{Z_j} \in {\Upsilon _k}|{Y_j},{\Theta ^{\left( i \right)}}} \right)$即可,${{Z_j} \in {\Upsilon _k}}$表示第$j$个观测点来自第$k$个分模型。同硬币第三抛的求解完全一致,利用全概率公式,容易得到:

为了推导简洁,M-Step时保留隐变量概率的原形式而不再展开。

2)构造准则函数Q

 根据上面给出的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]$为分布$k$对应的参数,$\Theta$  = {$\theta _1$,$\theta _2$,...,$\theta _K$}为参数集合,$N$为样本个数,$K$为混合模型个数。

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

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$为向量或矩阵,对于高斯分布:

关于$\theta_k$利用MLE即可求参,注意对${{{\bf{\Sigma }}_k}}$求偏导时,关于${{{\bf{\Sigma }}^{-1}_k}}$求偏导更方便,得出结果:

至此,完成了参数求导,可见推导前半部分对于任意分布都有效。只是涉及具体求参时,形式不同有差别。

总结一下GMM:

E-Step:

M-Step:

 

三、GMM代码实现

 子程序代码:

function [u,sig,t,iter] = fit_mix_gaussian( X,M )
%
% fit_mix_gaussian - fit parameters for a mixed-gaussian distribution using EM algorithm
%
% format:   [u,sig,t,iter] = fit_mix_gaussian( 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
Q           = zeros(N,M);                       % user variable to determine when we have converged to a steady solution
thresh      = 1e-3;
step        = N;
last_step   = inf;
iter        = 0;
min_iter    = 10;

% main convergence loop, assume gaussians are 1D
while ((( abs((step/last_step)-1) > thresh) & (step>(N*eps)) ) | (iter<min_iter) ) 
    
    % E step
    % ========
    Q   = Z;
    P   = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - 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
    u               = (X')*Z ./ Zm;
    sig2            = sum(((X*Ir - Ic*u).^2).*Z) ./ Zm;
    t               = Zm/N;
end
sig     = sqrt( sig2 );

给出一个示例:

clc;clear all;close all;
set(0,'defaultfigurecolor','w') 
x = [1*randn(100000,1)+3;3*randn(100000,1)-5];
%fitting
x       = x(:);                 % should be column vectors !
N       = length(x);
[u,sig,t,iter] = fit_mix_gaussian( x,2 );
sig = sig.^2;
%Plot
figure;
%Bar
subplot 221
plot(x(randperm(N)),'k');grid on;
xlim([0,N]);
subplot 222
numter = [-15:.2:10];
[histFreq, histXout] = hist(x, numter);
binWidth = histXout(2)-histXout(1);
bar(histXout, histFreq/binWidth/sum(histFreq)); hold on;grid on;
%Fitting plot
subplot 223
y = t(2)*1/sqrt(2*pi*sig(2))*exp(-(numter-u(2)).^2/2/sig(2));
plot(numter,y,'r','linewidth',2);grid on;
hold on;
y = t(1)*1/sqrt(2*pi*sig(1))*exp(-(numter-u(1)).^2/2/sig(1));
plot(numter,y,'g','linewidth',2);grid on;

%Fitting result
subplot 224
bar(histXout, histFreq/binWidth/sum(histFreq)); hold on;grid on;
y = t(2)*1/sqrt(2*pi*sig(2))*exp(-(numter-u(2)).^2/2/sig(2));
plot(numter,y,'r','linewidth',2);grid on;
hold on;
y = t(1)*1/sqrt(2*pi*sig(1))*exp(-(numter-u(1)).^2/2/sig(1));
plot(numter,y,'g','linewidth',2);grid on;

结果便是GMM背景介绍中的图形。

类似的,可以参考混合拉普拉斯分布拟合(LMM),对应效果:

参考:

李航《统计学习方法》.

posted @ 2017-03-21 07:13 桂。 阅读(...) 评论(...) 编辑 收藏