1.算法仿真效果
matlab2022a仿真结果如下：
2.算法涉及理论知识概要
GMM，高斯混合模型，也可以简写为MOG。高斯模型就是用高斯概率密度函数（正态分布曲线）精确地量化事物，将一个事物分解为若干的基于高斯概率密度函数（正态分布曲线）形成的模型。GMMs已经在数值逼近、语音识别、图像分类、图像去噪、图像重构、故障诊断、视频分析、邮件过滤、密度估计、目标识别与跟踪等领域取得了良好的效果 。
高斯混合模型 (GMM) 是一种机器学习算法。它们用于根据概率分布将数据分类为不同的类别。高斯混合模型可用于许多不同的领域，包括金融、营销等等！这里要对高斯混合模型进行介绍以及真实世界的示例、它们的作用以及何时应该使用GMM。
高斯混合模型 (GMM) 是一个概率概念，用于对真实世界的数据集进行建模。GMM是高斯分布的泛化，可用于表示可聚类为多个高斯分布的任何数据集。高斯混合模型是一种概率模型，它假设所有数据点都是从具有未知参数的高斯分布的混合中生成的。
高斯混合模型可用于聚类，这是将一组数据点分组为聚类的任务。GMM 可用于在数据集中可能没有明确定义的集群中查找集群。此外，GMM 可用于估计新数据点属于每个集群的概率。高斯混合模型对异常值也相对稳健，这意味着即使有一些数据点不能完全适合任何集群，它们仍然可以产生准确的结果。这使得 GMM 成为一种灵活而强大的数据聚类工具。它可以被理解为一个概率模型，其中为每个组假设高斯分布，并且它们具有定义其参数的均值和协方差。
GMM 由两部分组成——均值向量 (μ) 和协方差矩阵 (Σ)。高斯分布被定义为呈钟形曲线的连续概率分布。高斯分布的另一个名称是正态分布。这是高斯混合模型的图片：它可以被理解为一个概率模型，其中为每个组假设高斯分布，并且它们具有定义其参数的均值和协方差。GMM 由两部分组成——均值向量 (μ) 和协方差矩阵 (Σ)。高斯分布被定义为呈钟形曲线的连续概率分布。高斯分布的另一个名称是正态分布。这是高斯混合模型的图片：
EM (Expectation Maximization)算法是由Dempster、Laind和Rubin在1977年提出的一种求参数的极大似然估计方法，可以广泛地应用于处理缺损数据、截尾数据等带有噪声的不完整数据。针对不完整数据集，EM算法主要应用于以下两种情况的参数估计：第一，由于观测过程中本身的错误或局限性导致的观测数据自身不完整；第二，数据没有缺失，但是无法得到似然函数的解析解，或似然函数过于复杂，难以直接优化分析，而引入额外的缺失参数能使得简化后的似然函数便于参数估计。
最大期望算法（Expectation-Maximization algorithm, EM），或Dempster-Laird-Rubin算法 ，是一类通过迭代进行极大似然估计（Maximum Likelihood Estimation, MLE）的优化算法，通常作为牛顿迭代法（Newton-Raphson method）的替代用于对包含隐变量（latent variable）或缺失数据（incomplete-data）的概率模型进行参数估计 。
EM算法的标准计算框架由E步（Expectation-step）和M步（Maximization step）交替组成，算法的收敛性可以确保迭代至少逼近局部极大值。EM算法是MM算法（Minorize-Maximization algorithm）的特例之一，有多个改进版本，包括使用了贝叶斯推断的EM算法、EM梯度算法、广义EM算法等。
由于迭代规则容易实现并可以灵活考虑隐变量，EM算法被广泛应用于处理数据的缺测值 ，以及很多机器学习（machine learning）算法，包括高斯混合模型（Gaussian Mixture Model, GMM）和隐马尔可夫模型（Hidden Markov Model, HMM）的参数估计。
3.MATLAB核心程序
z2 = gaussian2D(gridX, mu2, sigma2);
Z1 = reshape(z1, gridSize, gridSize);
Z2 = reshape(z2, gridSize, gridSize);
[C, h] = contour(u, u, Z1);
[C, h] = contour(u, u, Z2);
axis([-6 6 -6 6])
title('Original Data and PDFs');
m = size(X, 1);
k = 2;
n = 2;
indeces = randperm(m);
mu = X(indeces(1:k), :);
sigma = [];
for (j = 1 : k)
sigma{j} = cov(X);
end
phi = ones(1, k) * (1 / k);
%Run Expectation Maximization
W = zeros(m, k);
for (iter = 1:1000)
%Expectation
pdf = zeros(m, k);
for (j = 1 : k)
pdf(:, j) = gaussian2D(X, mu(j, :), sigma{j});
end
pdf_w = bsxfun(@times, pdf, phi);
W = bsxfun(@rdivide, pdf_w, sum(pdf_w, 2));
%Maximization
prevMu = mu;
for (j = 1 : k)
phi(j) = mean(W(:, j), 1);
mu(j, :) = weightedAverage(W(:, j), X);
sigma_k = zeros(n, n);
Xm = bsxfun(@minus, X, mu(j, :));
for (i = 1 : m)
sigma_k = sigma_k + (W(i, j) .* (Xm(i, :)' * Xm(i, :)));
end
sigma{j} = sigma_k ./ sum(W(:, j));
end
if (mu == prevMu)
break
end
end
figure(2);
hold off;
plot(X1(:, 1), X1(:, 2), 'bo');
hold on;
plot(X2(:, 1), X2(:, 2), 'ro');