參考文章
http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/
《統計學習方法》 李航
EM算法是一種疊代算法，1977年由Dempster等人總結出，用于含有隐變量的機率模型參數的極大似然估計，或極大後驗機率估計。EM算法的每次疊代由兩步組成：E步，求期望(expectation);M步，求極大(maximization).是以這一算法稱為期望極大算法(expectation maximization algorithm),簡稱EM算法。
首先準備一些預備知識，如：二維高斯混合函數的形式，協方差的定義
高斯混合函數
對于一維混合高斯函數，如(1)式所示。(所謂的混合高斯，就是将一些高斯函數利用一些權值系數組合起來) （1）
二維高斯混合函數，如(2)式所示。
（2）
可以看出二維高斯混合模型，其中的d為資料的次元，在這裡等于2；∑ 代表協方差矩陣，兩行兩列；這裡的x以及μ都是二維的資料，用矩陣表示就是一行兩列。這裡提到了協方差矩陣，那麼什麼是協方差矩陣，與方差有什麼差別？
協方差矩陣
了解協方差矩陣的關鍵在于牢記它計算的是不同次元之間的協方差，而不是不同樣本之間。拿到一個樣本矩陣，我們最先明确的是一行是一個樣本還是一個次元（重要）。
方差表示的是每個樣本值與全體樣本值的平均數之差的平方值的平均數，反映了一組資料的離散程度。
既然有了方差還要協方差幹啥？顯然，方差隻能衡量一維資料的離散程度，對于二維資料，就用到了協方差，如下式所示。其中X，Y代表一組資料的兩個次元。COV(X,Y)的值大于零，表示兩個次元的資料成正相關，小于零則表示成負相關，等于零表示兩個次元的資料互相獨立。顯然當X=Y時，COV(X,X)計算的就是方差。由協方差的定義可以得到COV(X,Y)=COV(Y,X)。
二維資料的協方差矩陣表示為 ，顯然協方差矩陣為對角矩陣。
三維資料的協方差矩陣表示為，同理知道N維資料的協方差矩陣。
EM算法
假設一組資料x1,x2,...,xn是來自上述二維高斯混合模型，這裡我們利用EM算法确定各個高斯部件的參數，充分拟合給定的資料。算法目的具體的理論知識與算法步驟參考《統計學習方法》李航一書中對于EM算法的講解(李航的這一本書，絕對值得學習)。
我主要寫一下自己的了解，也為自己以後複習用。
EM算法其實就是機率模型參數的極大似然估計法，隻不過存在一些隐變量，進而不能像一般的機率模型參數估計一樣直接求似然函數的最大值，因為這裡的似然函數存在隐變量，因而不存在解析解，就不能直接對 似然函數求偏導來求極大值，隻能疊代求解。
基本思想：
E步，假設參數已知，去估計隐藏變量的期望；
M步，利用E步求得的隐藏變量的期望，根據最大似然估計求得參數，然後新得到的參數值重新被用于E步。。。直至收斂到局部最優解。
需要注意的一點是：理論中，E步最重要的是求Q函數（李航 9.29式），即完全資料的對數似然函數的期望，M步要做的是用Q函數對各參數求偏導，進而确定滿足Q函數最大的各參數的值。而在程式中，因為M步Q函數對參數求偏導後，各參數的值由Q函數中的Gamma(j,k)決定，是以E步隻需計算Gamma(j,k)即可，Gamma(j,k)代表第j個樣本屬于第k個高斯模型的機率，這樣M步各參數的值也就确定了，是以切記執行個體理論與程式實作的差別，好像是程式隻是利用了理論的精簡結果。
學習完理論之後，最好寫個程式練習一下，加深自己的了解。下面的程式(matlab編寫)，是這樣一個過程：
1.首先産生兩組資料點。每一組資料點都是符合給定均值Mu和協方差矩陣∑的。程式和圖如下所示，
%産生2個二維正态資料
MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];
X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
scatter(X(:,1),X(:,2),10,'*');
2.對于上述産生的資料，利用EM算法預測這組資料對應的均值μ和協方差矩陣∑
對于各參數初值的設定：兩個均值初始值各設定為兩組資料點中的随機一個樣本點，權值系數初始值設定為1/k,k是高斯混合模型中模型的個數，協方差矩陣初值設定為所有資料的協方差。
程式：GMM.m
%程式中Psi對應二維高斯函數，其中的(x-μ)與其轉置的順序與上述提到的高斯函數不同，這樣是為了保證矩陣可乘性，不影響結果.
%Gamma 為隐變量的值，Gamma(i,j)代表第i個樣本屬于第j個模型的機率。
%Mu為期望，Sigma為協方差矩陣%Pi為各模型的權值系數%産生2個二維正态資料
%産生2個二維正态資料
MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];
%生成1000行2列(預設)均值為mu标準差為sigma的正态分布随機數
X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
%顯示
scatter(X(:,1),X(:,2),10,'.');
%====================
K=2;
[N,D]=size(X);
Gamma=zeros(N,K);
Psi=zeros(N,K);
Mu=zeros(K,D);
LM=zeros(K,D);
Sigma =zeros(D, D, K);
Pi=zeros(1,K);
%選擇随機的兩個樣本點作為期望疊代初值
Mu(1,:)=X(randint(1,1,[1 N/2]),:);
Mu(2,:)=X(randint(1,1,[N/2 N]),:);
%所有資料的協方差作為協方差初值
for k=1:K
Pi(k)=1/K;
Sigma(:, :, k)=cov(X);
end
LMu=Mu;
LSigma=Sigma;
LPi=Pi;
while true
%Estimation Step
for k = 1:K
Y = X - repmat(Mu(k,:),N,1);
Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y'))); %Psi第一列代表第一個高斯分布對所有資料的取值
end
Gamma_SUM=zeros(D,D);
for j = 1:N
for k=1:K
Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi'); %Psi的第一行分别代表兩個高斯分布對第一個資料的取值
end
end
%Maximization Step
for k = 1:K
%update Mu
Mu_SUM= zeros(1,D);
for j=1:N
Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);
end
Mu(k,:)= Mu_SUM/sum(Gamma(:,k));
%update Sigma
Sigma_SUM= zeros(D,D);
for j = 1:N
Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));
end
Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));
%update Pi
Pi_SUM=0;
for j=1:N
Pi_SUM=Pi_SUM+Gamma(j,k);
end
Pi(1,k)=Pi_SUM/N;
end
R_Mu=sum(sum(abs(LMu- Mu)));
R_Sigma=sum(sum(abs(LSigma- Sigma)));
R_Pi=sum(sum(abs(LPi- Pi)));
R=R_Mu+R_Sigma+R_Pi;
if ( R<1e-10)
disp('期望');
disp(Mu);
disp('協方差矩陣');
disp(Sigma);
disp('權值系數');
disp(Pi);
obj=gmdistribution(Mu,Sigma);
figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);
break;
end
LMu=Mu;
LSigma=Sigma;
LPi=Pi;
end
%=====================
%%%%%%%%運作結果%%%%%%%
%
% 期望
% 1.0334 2.0479
% -0.9703 -0.9879
%
% 協方差矩陣
%
% (:,:,1) =
%
% 0.9604 -0.0177
% -0.0177 0.4463
%
%
% (:,:,2) =
%
% 0.9976 0.0667
% 0.0667 1.0249
%
% 權值系數
% 0.4935 0.5065
%%%%%%%%%%%%%%%%%%%%%%%%%
參考文章
http://blog.163.com/baolong_zhu/blog/static/196311091201421185531966/
《統計學習方法》 李航
EM算法是一種疊代算法，1977年由Dempster等人總結出，用于含有隐變量的機率模型參數的極大似然估計，或極大後驗機率估計。EM算法的每次疊代由兩步組成：E步，求期望(expectation);M步，求極大(maximization).是以這一算法稱為期望極大算法(expectation maximization algorithm),簡稱EM算法。
首先準備一些預備知識，如：二維高斯混合函數的形式，協方差的定義
高斯混合函數
對于一維混合高斯函數，如(1)式所示。(所謂的混合高斯，就是将一些高斯函數利用一些權值系數組合起來) （1）
二維高斯混合函數，如(2)式所示。
（2）
可以看出二維高斯混合模型，其中的d為資料的次元，在這裡等于2；∑ 代表協方差矩陣，兩行兩列；這裡的x以及μ都是二維的資料，用矩陣表示就是一行兩列。這裡提到了協方差矩陣，那麼什麼是協方差矩陣，與方差有什麼差別？
協方差矩陣
了解協方差矩陣的關鍵在于牢記它計算的是不同次元之間的協方差，而不是不同樣本之間。拿到一個樣本矩陣，我們最先明确的是一行是一個樣本還是一個次元（重要）。
方差表示的是每個樣本值與全體樣本值的平均數之差的平方值的平均數，反映了一組資料的離散程度。
既然有了方差還要協方差幹啥？顯然，方差隻能衡量一維資料的離散程度，對于二維資料，就用到了協方差，如下式所示。其中X，Y代表一組資料的兩個次元。COV(X,Y)的值大于零，表示兩個次元的資料成正相關，小于零則表示成負相關，等于零表示兩個次元的資料互相獨立。顯然當X=Y時，COV(X,X)計算的就是方差。由協方差的定義可以得到COV(X,Y)=COV(Y,X)。
二維資料的協方差矩陣表示為 ，顯然協方差矩陣為對角矩陣。
三維資料的協方差矩陣表示為，同理知道N維資料的協方差矩陣。
EM算法
假設一組資料x1,x2,...,xn是來自上述二維高斯混合模型，這裡我們利用EM算法确定各個高斯部件的參數，充分拟合給定的資料。算法目的具體的理論知識與算法步驟參考《統計學習方法》李航一書中對于EM算法的講解(李航的這一本書，絕對值得學習)。
我主要寫一下自己的了解，也為自己以後複習用。
EM算法其實就是機率模型參數的極大似然估計法，隻不過存在一些隐變量，進而不能像一般的機率模型參數估計一樣直接求似然函數的最大值，因為這裡的似然函數存在隐變量，因而不存在解析解，就不能直接對 似然函數求偏導來求極大值，隻能疊代求解。
基本思想：
E步，假設參數已知，去估計隐藏變量的期望；
M步，利用E步求得的隐藏變量的期望，根據最大似然估計求得參數，然後新得到的參數值重新被用于E步。。。直至收斂到局部最優解。
需要注意的一點是：理論中，E步最重要的是求Q函數（李航 9.29式），即完全資料的對數似然函數的期望，M步要做的是用Q函數對各參數求偏導，進而确定滿足Q函數最大的各參數的值。而在程式中，因為M步Q函數對參數求偏導後，各參數的值由Q函數中的Gamma(j,k)決定，是以E步隻需計算Gamma(j,k)即可，Gamma(j,k)代表第j個樣本屬于第k個高斯模型的機率，這樣M步各參數的值也就确定了，是以切記執行個體理論與程式實作的差別，好像是程式隻是利用了理論的精簡結果。
學習完理論之後，最好寫個程式練習一下，加深自己的了解。下面的程式(matlab編寫)，是這樣一個過程：
1.首先産生兩組資料點。每一組資料點都是符合給定均值Mu和協方差矩陣∑的。程式和圖如下所示，
%産生2個二維正态資料
MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];
X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
scatter(X(:,1),X(:,2),10,'*');
2.對于上述産生的資料，利用EM算法預測這組資料對應的均值μ和協方差矩陣∑
對于各參數初值的設定：兩個均值初始值各設定為兩組資料點中的随機一個樣本點，權值系數初始值設定為1/k,k是高斯混合模型中模型的個數，協方差矩陣初值設定為所有資料的協方差。
程式：GMM.m
%程式中Psi對應二維高斯函數，其中的(x-μ)與其轉置的順序與上述提到的高斯函數不同，這樣是為了保證矩陣可乘性，不影響結果.
%Gamma 為隐變量的值，Gamma(i,j)代表第i個樣本屬于第j個模型的機率。
%Mu為期望，Sigma為協方差矩陣%Pi為各模型的權值系數%産生2個二維正态資料
%産生2個二維正态資料
MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];
%生成1000行2列(預設)均值為mu标準差為sigma的正态分布随機數
X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];
%顯示
scatter(X(:,1),X(:,2),10,'.');
%====================
K=2;
[N,D]=size(X);
Gamma=zeros(N,K);
Psi=zeros(N,K);
Mu=zeros(K,D);
LM=zeros(K,D);
Sigma =zeros(D, D, K);
Pi=zeros(1,K);
%選擇随機的兩個樣本點作為期望疊代初值
Mu(1,:)=X(randint(1,1,[1 N/2]),:);
Mu(2,:)=X(randint(1,1,[N/2 N]),:);
%所有資料的協方差作為協方差初值
for k=1:K
Pi(k)=1/K;
Sigma(:, :, k)=cov(X);
end
LMu=Mu;
LSigma=Sigma;
LPi=Pi;
while true
%Estimation Step
for k = 1:K
Y = X - repmat(Mu(k,:),N,1);
Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y'))); %Psi第一列代表第一個高斯分布對所有資料的取值
end
Gamma_SUM=zeros(D,D);
for j = 1:N
for k=1:K
Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi'); %Psi的第一行分别代表兩個高斯分布對第一個資料的取值
end
end
%Maximization Step
for k = 1:K
%update Mu
Mu_SUM= zeros(1,D);
for j=1:N
Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);
end
Mu(k,:)= Mu_SUM/sum(Gamma(:,k));
%update Sigma
Sigma_SUM= zeros(D,D);
for j = 1:N
Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));
end
Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));
%update Pi
Pi_SUM=0;
for j=1:N
Pi_SUM=Pi_SUM+Gamma(j,k);
end
Pi(1,k)=Pi_SUM/N;
end
R_Mu=sum(sum(abs(LMu- Mu)));
R_Sigma=sum(sum(abs(LSigma- Sigma)));
R_Pi=sum(sum(abs(LPi- Pi)));
R=R_Mu+R_Sigma+R_Pi;
if ( R<1e-10)
disp('期望');
disp(Mu);
disp('協方差矩陣');
disp(Sigma);
disp('權值系數');
disp(Pi);
obj=gmdistribution(Mu,Sigma);
figure,h = ezmesh(@(x,y)pdf(obj,[x,y]),[-8 6], [-8 6]);
break;
end
LMu=Mu;
LSigma=Sigma;
LPi=Pi;
end
%=====================
%%%%%%%%運作結果%%%%%%%
%
% 期望
% 1.0334 2.0479
% -0.9703 -0.9879
%
% 協方差矩陣
%
% (:,:,1) =
%
% 0.9604 -0.0177
% -0.0177 0.4463
%
%
% (:,:,2) =
%
% 0.9976 0.0667
% 0.0667 1.0249
%
% 權值系數
% 0.4935 0.5065
%%%%%%%%%%%%%%%%%%%%%%%%%
---------------------
作者：lpkinging
原文：https://blog.csdn.net/u011582757/article/details/70003351