EM算法matlab代码

由于EM算法是一种通用的统计算法，其代码的实现方式也因问题而异。以下是一个简单的例子，用于说明如何使用EM算法来估计两个高斯分布的参数。

假设我们有两个高斯分布，其均值和方差未知，我们希望通过观察一些样本来估计它们的参数。我们可以使用EM算法来实现这个目标。

首先，我们生成一些样本：

mu1 = 1;
sigma1 = 2;
mu2 = 4;
sigma2 = 1;
n1 = 100;
n2 = 150;
x1 = mu1 + sigma1 * randn(n1, 1);
x2 = mu2 + sigma2 * randn(n2, 1);
x = [x1; x2];

接下来，我们初始化两个高斯分布的参数：

mu1_hat = 1;
sigma1_hat = 1;
mu2_hat = 3;
sigma2_hat = 1;

然后，我们开始迭代，直到算法收敛为止：

while true
% E step
p1 = normpdf(x, mu1_hat, sigma1_hat);
p2 = normpdf(x, mu2_hat, sigma2_hat);
w1 = p1 ./ (p1 + p2);
w2 = p2 ./ (p1 + p2);

% M step
mu1_hat_new = sum(w1 .* x) / sum(w1);
sigma1_hat_new = sqrt(sum(w1 .* (x - mu1_hat_new).^2) / sum(w1));
mu2_hat_new = sum(w2 .* x) / sum(w2);
sigma2_hat_new = sqrt(sum(w2 .* (x - mu2_hat_new).^2) / sum(w2));

% check for convergence
if abs(mu1_hat_new - mu1_hat) < 1e-6 && abs(sigma1_hat_new - sigma1_hat) < 1e-6 && ...
abs(mu2_hat_new - mu2_hat) < 1e-6 && abs(sigma2_hat_new - sigma2_hat) < 1e-6
break;
end

% update parameters
mu1_hat = mu1_hat_new;
sigma1_hat = sigma1_hat_new;
mu2_hat = mu2_hat_new;
sigma2_hat = sigma2_hat_new;
end

最后，我们可以输出估计出的参数：

fprintf('mu1_hat = %f\n', mu1_hat);
fprintf('sigma1_hat = %f\n', sigma1_hat);
fprintf('mu2_hat = %f\n', mu2_hat);
fprintf('sigma2_hat = %f\n', sigma2_hat);

这段代码实现了一个简单的二元高斯分布的EM算法，以估计其参数。它可以很容易地推广到更多高斯分布或其他类型的分布。