作者：桂。
時間：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），對應效果：
參考：
李航《統計學習方法》.