混合高斯模型

玩了混合高斯模型，先轉幾個參考資料，曾經試過自己寫代碼，結果發現混合高斯模型矩陣運算對我的計算能力要求很高，結果失敗了，上網找了代碼學習一下大牛們的編程思想，事實證明數學寫出來的公式雖然很美，但是現實寫代碼的時候要考慮各種問題~~~

1.http://www.cnblogs.com/cfantaisie/archive/2011/08/20/2147075.html （主要實現代碼）

2.http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/ （對於奇異矩陣的正則化）

3.參考的一本書，Computer Vision Models, Learning, and Inference.（數學推導）

4.斯坦福機器學習課程（數學推導）

利用EM算法：

E-step：

 

 

 

 

M-step:

 

 

 

 

 

 

 

matlab實現代碼：

function prob = MOF_guassPdf(Data,Mu,Sigma)
%

% 根據高斯分布函數計算每組數據的概率密度 Probability Density Function (PDF) 

% 輸入 -----------------------------------------------------------------

%   o Data:  D x N ，N個D維數據

%   o Mu:    D x 1 ，M個Gauss模型的中心初始值

%   o Sigma: D x D ，每個Gauss模型的方差（假設每個方差矩陣都是對角陣，

%                                   即一個數和單位矩陣的乘積）

% Outputs ----------------------------------------------------------------

%   o prob:  N x 1 array representing the probabilities for the 

%            N datapoints.     

[dim,N] = size(Data);

Data = Data' - repmat(Mu',N,1);

prob = sum((Data/Sigma).*Data, 2);

prob = exp(-0.5*prob) / sqrt((2*pi)^dim * (abs(det(Sigma))+realmin));

EM迭代：

function [Alpha, Mu, Sigma] = MOF_EM(Data, Alpha0, Mu0, Sigma0)
%% EM 迭代停止條件
loglik_threshold = 1e-10;

%% 初始化參數

[dim, N] = size(Data); 

M = size(Mu0,2);

loglik_old = -realmax;

nbStep = 0;
                %Data是D*N的矩陣

Mu = Mu0;       %Mu是D*M的矩陣

Sigma = Sigma0; %sigma是D*D*M的三維矩陣

Alpha = Alpha0; %alpha是1*M大小的向量，意思是列代表第M個高斯模型。

Epsilon = 0.0001; 

while (nbStep < 1200)

  nbStep = nbStep+1;

  %% E-步驟 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:M

    % PDF of each point

    Pxi(:,i) = MOF_guassPdf(Data, Mu(:,i), Sigma(:,:,i));          

  end

  

  % 計算后驗概率 beta(i|x)

  Pix_tmp = repmat(Alpha,[N 1]).*Pxi; %Pxi應該是N*M的高斯概率矩陣，意思是第N個數據計算第M個高斯函數的概率

  Pix = Pix_tmp ./ (repmat(sum(Pix_tmp,2),[1 M])+realmin);

  Beta = sum(Pix);

  %% M-步驟 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:M

    % 更新權值

    Alpha(i) = Beta(i) / N;

    % 更新均值

    Mu(:,i) = Data*Pix(:,i) / Beta(i);

    % 更新方差

    Data_tmp1 = Data - repmat(Mu(:,i),1,N);

    Sigma(:,:,i) = (repmat(Pix(:,i)',dim, 1) .* Data_tmp1 * Data_tmp1') / Beta(i);

    %% Add a tiny variance to avoid numerical instability

    Sigma(:,:,i) = Sigma(:,:,i) + 1E-5.*diag(ones(dim,1));

  end

 

%  %% Stopping criterion 1 %%%%%%%%%%%%%%%%%%%%

  for i=1:M

    %Compute the new probability p(x|i)

    Pxi(:,i) = MOF_guassPdf(Data, Mu(:,i), Sigma(i));

  end

  %Compute the log likelihood

  F = Pxi*Alpha';

  F(find(F<realmin)) = realmin;

  loglik = mean(log(F));

  %Stop the process depending on the increase of the log likelihood 

  if abs((loglik/loglik_old)-1) < loglik_threshold

    break;

  end

  loglik_old = loglik;

 

  %% Stopping criterion 2 %%%%%%%%%%%%%%%%%%%%
%{
  v = [sum(abs(Mu - Mu0)), abs(Alpha - Alpha0)];

  s = abs(Sigma-Sigma0);

  v2 = 0;

  for i=1:M

    v2 = v2 + det(s(:,:,i));

  end

 

  if ((sum(v) + v2) < Epsilon)

    break;

  end

  Mu0 = Mu;

  Sigma0 = Sigma;

  Alpha0 = Alpha;
%}
end

測試結果

編寫代碼隨機生成3個高斯分布的數據，參數也隨機生成（注意sigma要半正定對稱）：

這里只畫了一個經過EM算法迭代所得參數的高斯函數圖..（有點丑，不知道怎么將mesh弄透明，求搜索關鍵詞）。