用matlab写一个EM算法的程序

以下是一个简单的EM算法的matlab程序，用于计算两个高斯分布的参数。该程序假设数据是从两个不同的高斯分布中生成的，但是我们不知道哪个高斯分布生成了哪些数据点。

% 生成数据
mu1 = [1; 2];
mu2 = [4; 5];
sigma1 = [1 0.5; 0.5 2];
sigma2 = [2 -1; -1 1];
N = 100;
X = [mvnrnd(mu1, sigma1, N); mvnrnd(mu2, sigma2, N)]';

% 初始化参数
K = 2;
pi = [0.5; 0.5];
mu = [randn(2,1) randn(2,1)];
sigma(:,:,1) = eye(2);
sigma(:,:,2) = eye(2);

% EM算法
for iter = 1:100
% E步骤：计算后验概率
for k = 1:K
P(:,k) = pi(k) * mvnpdf(X', mu(:,k)', sigma(:,:,k));
end
P = P ./ sum(P,2);

% M步骤：更新参数
pi = mean(P,1)';
mu = X * P ./ sum(P,1)';
for k = 1:K
Xshift = X - mu(:,k)';
sigma(:,:,k) = (Xshift' * (Xshift .* repmat(P(:,k),1,2))) ./ sum(P(:,k),1);
end

% 计算似然函数
loglik(iter) = sum(log(sum(P,2)));
end

% 绘制数据点和高斯分布
figure
scatter(X(1,:), X(2,:))
hold on
x1 = linspace(-5, 10);
x2 = linspace(-5, 10);
[X1, X2] = meshgrid(x1, x2);
F1 = mvnpdf([X1(:) X2(:)], mu1', sigma1);
F1 = reshape(F1, length(x2), length(x1));
F2 = mvnpdf([X1(:) X2(:)], mu2', sigma2);
F2 = reshape(F2, length(x2), length(x1));
contour(x1, x2, F1, [0.001 0.01 0.05 0.1 0.2], 'LineWidth', 2)
contour(x1, x2, F2, [0.001 0.01 0.05 0.1 0.2], 'LineWidth', 2)
legend('Data', 'Gaussian 1', 'Gaussian 2')
xlabel('X1')
ylabel('X2')

% 绘制似然函数
figure
plot(loglik, 'LineWidth', 2)
xlabel('Iteration')
ylabel('Log-likelihood')