混合高斯模型

玩了混合高斯模型，先转几个参考资料，曾经试过自己写代码，结果发现混合高斯模型矩阵运算对我的计算能力要求很高，结果失败了，上网找了代码学习一下大牛们的编程思想，事实证明数学写出来的公式虽然很美，但是现实写代码的时候要考虑各种问题~~~

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弄透明，求搜索关键词）。