% EM测试用于模拟k个正态分布的混合模型的参数估计
% 中国民航大学计算机学院，张志远
% 2010年2月27日
%---------------------------------------------
%生成符合条件的高斯分布的参数
%---------------------------------------------
% 首先指定两个高斯分布的参数
Miu1 = 3;
Miu2 = -2;
Sigma1 = 1;
Sigma2 = 2;
Alpha1 = 0.4;
Alpha2 = 0.6;
% 然后alpha1和alpha2的概率选择两种高斯分布，并生成样本值
N = 500;
X = zeros(1, N);
N1 = floor(N * Alpha1);
N2 = N - N1;
X(1:N1) = randn(1,N1)*Sigma1 + Miu1;
X(N1+1:N) = randn(1,N2)*Sigma2 + Miu2;
hist(X, 50)
%-------------------------------------------------
% EM算法，步骤1， 计算均值
%------------------------------------------------
% 首先假设两个高斯分布的均值为随机的两个值，注意此处的Sigma是Sigma的平方
M = 2;
Miu = randn(1, M);
Sigma = rand(1, M);
Alpha = [0.2, 0.8]
PkMatrix = zeros(N, M);
for step = 1: 100
for i = 1 : N
for k = 1 : M
PkMatrix(i, k) = exp(-1*(X(1,i)-Miu(1,k))^2/(2*Sigma(1,k)))/sqrt(2*pi*Sigma(1,k));
end
end
% SumPkMatrixI 是N*1的矩阵，其每一行的值等于 P1i + P2i + P3i + ... + PMi
% 即SumPkMatrixI 是对PkMatrix按行相加的结果
SumPkMatrixI = sum(PkMatrix, 2);
for i = 1 : N
PkMatrix(i,:) = PkMatrix(i,:)/SumPkMatrixI(i,1);
end
%----------------------------------------------
% EM算法，步骤2， 求最大化的参数Miu, Sigma 和 Alpha
%----------------------------------------------
OldMiu = Miu;
OldSigma = Sigma;
OldAlpha = Alpha;
for k = 1 : M
P1 = 0;
P2 = 0;
for i = 1 : N
P1 = P1 + PkMatrix(i,k);
P2 = P2 + PkMatrix(i,k) * X(1,i);
end
Miu(1,k) = P2/P1;
Alpha(1,k) = P1/N;
%注意必须先计算好Miu才能计算P3，
%因为要计算Sigma用到的P3依赖于的是新的Miu
P3 = 0;
for i = 1 : N
P3 = P3 + PkMatrix(i,k) * (X(1,i) - Miu(1,k))^2;
end
Sigma(1,k) = P3/P1;
end
%-----------------------------------------------------
% 判断假设是否发生了足够的变化
%-------------------------------------------------------
Epsilon = 0.00001;
if sum(abs(Miu - OldMiu)) < Epsilon && ...
sum(abs(Sigma - OldSigma)) < Epsilon && ...
sum(abs(Alpha - OldAlpha)) < Epsilon
break;
end
disp('========================================');
step
Miu
Sigma
Alpha
end
