EM算法 matlab

以下是一个使用EM算法进行高斯混合模型（GMM）参数估计的Matlab代码示例：

% 生成数据
rng(1); % 设置随机数种子，以确保结果可重复
mu1 = [1 2];
sigma1 = [1 0; 0 1];
mu2 = [-3 -5];
sigma2 = [2 0; 0 2];
X = [mvnrnd(mu1, sigma1, 100); mvnrnd(mu2, sigma2, 100)];

% 初始化参数
K = 2; % GMM的组件数
N = size(X, 1); % 数据点数目
pi = ones(1, K) / K; % 混合系数
mu = [1 1; -1 -1]; % 均值
sigma = repmat(eye(2), 1, 1, K); % 协方差矩阵

% EM算法迭代
max_iter = 100; % 最大迭代次数
tol = 1e-6; % 收敛阈值
log_likelihood = zeros(max_iter, 1); % 保存每次迭代的对数似然值
for iter = 1:max_iter
    % E步：计算后验概率
    gamma = zeros(N, K);
    for k = 1:K
        gamma(:, k) = pi(k) * mvnpdf(X, mu(k, :), sigma(:, :, k));
    end
    gamma = gamma ./ sum(gamma, 2);
    
    % M步：更新参数
    Nk = sum(gamma, 1);
    pi = Nk / N;
    for k = 1:K
        mu(k, :) = sum(gamma(:, k) .* X) / Nk(k);
        X_centered = X - mu(k, :);
        sigma(:, :, k) = (X_centered' * (gamma(:, k) .* X_centered)) / Nk(k);
    end
    
    % 计算对数似然值
    log_likelihood(iter) = sum(log(sum(bsxfun(@times, pi, mvnpdf(X, mu, sigma))), 2));
    
    % 判断是否收敛
    if iter > 1 && abs(log_likelihood(iter) - log_likelihood(iter-1)) < tol
        break;
    end
end

% 绘制数据点和估计的高斯分布
figure;
scatter(X(:, 1), X(:, 2), 'filled');
hold on;
for k = 1:K
    plot_gaussian_ellipsoid(mu(k, :), sigma(:, :, k));
end
hold off;
xlabel('x');
ylabel('y');
title('GMM Parameter Estimation using EM Algorithm');

% 绘制对数似然值随迭代次数的变化
figure;
plot(1:iter, log_likelihood(1:iter));
xlabel('Iteration');
ylabel('Log Likelihood');
title('Log Likelihood vs. Iteration');

% 高斯分布椭圆绘制函数
function plot_gaussian_ellipsoid(mu, sigma)
    [V, D] = eig(sigma);
    t = linspace(0, 2*pi);
    a = (V * sqrt(D)) * [cos(t(:))'; sin(t(:))'];
    plot(mu(1) + a(1, :), mu(2) + a(2, :), 'r');
end

这段代码演示了如何使用EM算法对一个由两个高斯分布组成的数据集进行参数估计，并绘制出估计的高斯分布。代码首先生成一个由两个高斯分布组成的数据集，然后使用EM算法迭代估计混合系数、均值和协方差矩阵，直到收敛或达到最大迭代次数。最后，代码绘制出数据点和估计的高斯分布，并绘制出对数似然值随迭代次数的变化曲线。