每日算法 — 使用java实现EM算法：期望最大化与高斯混合模型参数估计

引言

期望最大化（Expectation-Maximization, EM）算法是统计机器学习中最优雅的参数估计方法之一，由Dempster、Laird和Rubin于1977年系统提出。EM算法通过迭代执行期望步（E步）和最大化步（M步），有效解决含有隐变量的概率模型参数估计问题。本文将用Java完整实现EM算法，并以高斯混合模型（GMM）为例，深入讲解隐变量推断与参数优化过程。

算法原理

问题背景

当概率模型包含隐变量（未观测变量）时，直接使用最大似然估计会遇到困难。EM算法通过引入隐变量的后验分布，将复杂的优化问题分解为两个交替执行的简单步骤。

数学推导

观测数据：$X = {x_1, x_2, …, x_n}$

隐变量：$Z = {z_1, z_2, …, z_n}$

目标：最大化对数似然函数

$$\log P(X|\theta) = \log \sum_Z P(X, Z|\theta)$$

由于log-sum结构难以直接优化，EM算法使用Jensen不等式构造下界：

$$\log P(X|\theta) \geq \sum_Z Q(Z) \log \frac{P(X, Z|\theta)}{Q(Z)} = \mathcal{L}(Q, \theta)$$

E步：固定参数 $\theta^{(t)}$，计算隐变量的后验分布

$$Q^{(t+1)}(Z) = P(Z|X, \theta^{(t)})$$

M步：固定 $Q^{(t+1)}$，最大化下界

$$\theta^{(t+1)} = \arg\max_\theta \sum_Z Q^{(t+1)}(Z) \log P(X, Z|\theta)$$

高斯混合模型（GMM）

GMM模型

$$P(x) = \sum_{k=1}^{K} \pi_k \mathcal{N}(x|\mu_k, \Sigma_k)$$

其中：
– $K$ 为高斯成分数
– $\pi_k$ 为混合系数（$\sum_k \pi_k = 1$）
– $\mu_k$ 为均值向量
– $\Sigma_k$ 为协方差矩阵
– $\mathcal{N}$ 为高斯分布

GMM的EM算法

E步：计算每个样本属于每个高斯成分的后验概率（责任值）

$$\gamma_{ik} = \frac{\pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(x_i|\mu_j, \Sigma_j)}$$

M步：使用责任值重新估计参数

$$N_k = \sum_{i=1}^{n} \gamma_{ik}$$

$$\pi_k = \frac{N_k}{n}$$

$$\mu_k = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik} x_i$$

$$\Sigma_k = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik} (x_i – \mu_k)(x_i – \mu_k)^T$$

Java实现

项目结构

src/
├── GaussianMixtureModel.java // GMM核心类
├── GaussianComponent.java    // 单个高斯成分
├── EMOptimizer.java          // EM优化器
├── DataPoint.java            // 数据点类
├── EMDemo.java               // 演示入口
└── DatasetGenerator.java     // 数据集生成器

单个高斯成分

/**
 * 高斯分布成分
 * 包含均值向量、协方差矩阵和混合系数
 */
public class GaussianComponent {
    private double[] mean;          // 均值向量
    private double[][] covariance;  // 协方差矩阵
    private double weight;          // 混合系数 π_k
    private int dimension;

    public GaussianComponent(int dimension) {
        this.dimension = dimension;
        this.mean = new double[dimension];
        this.covariance = new double[dimension][dimension];
        this.weight = 1.0;
    }

    /**
     * 计算多维高斯分布的概率密度
     */
    public double probabilityDensity(double[] x) {
        double[] diff = new double[dimension];
        for (int i = 0; i < dimension; i++) {
            diff[i] = x[i] - mean[i];
        }

        // 计算协方差矩阵的逆（简化：假设对角协方差矩阵）
        double det = 1.0;
        double[] invCovDiag = new double[dimension];
        for (int i = 0; i < dimension; i++) {
            det *= covariance[i][i];
            invCovDiag[i] = 1.0 / covariance[i][i];
        }

        // 计算马氏距离: (x-μ)^T Σ^{-1} (x-μ)
        double mahalanobis = 0;
        for (int i = 0; i < dimension; i++) {
            mahalanobis += diff[i] * diff[i] * invCovDiag[i];
        }

        // 高斯概率密度
        double norm = Math.pow(2 * Math.PI, dimension / 2.0) * Math.sqrt(det);
        return Math.exp(-0.5 * mahalanobis) / norm;
    }

    /**
     * 使用Cholesky分解处理一般协方差矩阵的概率密度计算
     * （简化版本，实际应使用数值稳定的矩阵运算库）
     */
    public double probabilityDensityStable(double[] x) {
        // 添加正则化防止奇异矩阵
        double[][] regCov = new double[dimension][dimension];
        for (int i = 0; i < dimension; i++) {
            System.arraycopy(covariance[i], 0, regCov[i], 0, dimension);
            regCov[i][i] += 1e-6; // 正则化
        }

        double[] diff = new double[dimension];
        for (int i = 0; i < dimension; i++) {
            diff[i] = x[i] - mean[i];
        }

        // 解线性方程组 Σ^{-1} * diff
        double[] invDiff = solveLinearSystem(regCov, diff);

        double mahalanobis = 0;
        for (int i = 0; i < dimension; i++) {
            mahalanobis += diff[i] * invDiff[i];
        }

        double det = computeDeterminant(regCov);
        double norm = Math.pow(2 * Math.PI, dimension / 2.0) * Math.sqrt(det);
        return Math.exp(-0.5 * mahalanobis) / norm;
    }

    /**
     * 使用高斯消元法解线性方程组
     */
    private double[] solveLinearSystem(double[][] A, double[] b) {
        int n = A.length;
        double[][] M = new double[n][n + 1];
        for (int i = 0; i < n; i++) {
            System.arraycopy(A[i], 0, M[i], 0, n);
            M[i][n] = b[i];
        }

        // 高斯消元
        for (int i = 0; i < n; i++) {
            // 选主元
            int maxRow = i;
            for (int k = i + 1; k < n; k++) {
                if (Math.abs(M[k][i]) > Math.abs(M[maxRow][i])) {
                    maxRow = k;
                }
            }
            double[] temp = M[i];
            M[i] = M[maxRow];
            M[maxRow] = temp;

            // 消元
            for (int k = i + 1; k < n; k++) {
                double factor = M[k][i] / M[i][i];
                for (int j = i; j <= n; j++) {
                    M[k][j] -= factor * M[i][j];
                }
            }
        }

        // 回代
        double[] x = new double[n];
        for (int i = n - 1; i >= 0; i--) {
            x[i] = M[i][n];
            for (int j = i + 1; j < n; j++) {
                x[i] -= M[i][j] * x[j];
            }
            x[i] /= M[i][i];
        }
        return x;
    }

    /**
     * 计算行列式（简化版，假设正定矩阵）
     */
    private double computeDeterminant(double[][] A) {
        int n = A.length;
        double[][] M = new double[n][n];
        for (int i = 0; i < n; i++) {
            System.arraycopy(A[i], 0, M[i], 0, n);
        }

        double det = 1;
        for (int i = 0; i < n; i++) {
            det *= M[i][i];
            for (int j = i + 1; j < n; j++) {
                double factor = M[j][i] / M[i][i];
                for (int k = i; k < n; k++) {
                    M[j][k] -= factor * M[i][k];
                }
            }
        }
        return det;
    }

    public void setMean(double[] mean) {
        this.mean = mean.clone();
    }

    public void setCovariance(double[][] covariance) {
        this.covariance = new double[dimension][dimension];
        for (int i = 0; i < dimension; i++) {
            System.arraycopy(covariance[i], 0, this.covariance[i], 0, dimension);
        }
    }

    public void setWeight(double weight) {
        this.weight = weight;
    }

    public double[] getMean() {
        return mean.clone();
    }

    public double[][] getCovariance() {
        double[][] copy = new double[dimension][dimension];
        for (int i = 0; i < dimension; i++) {
            System.arraycopy(covariance[i], 0, copy[i], 0, dimension);
        }
        return copy;
    }

    public double getWeight() {
        return weight;
    }

    @Override
    public String toString() {
        StringBuilder sb = new StringBuilder();
        sb.append("  混合系数: ").append(String.format("%.4f%n", weight));
        sb.append("  均值: [");
        for (int i = 0; i < dimension; i++) {
            sb.append(String.format("%.4f", mean[i]));
            if (i < dimension - 1) sb.append(", ");
        }
        sb.append("]\n");
        sb.append("  协方差对角: [");
        for (int i = 0; i < dimension; i++) {
            sb.append(String.format("%.4f", covariance[i][i]));
            if (i < dimension - 1) sb.append(", ");
        }
        sb.append("]");
        return sb.toString();
    }
}

GMM核心类

import java.util.ArrayList;
import java.util.List;
import java.util.Random;

/**
 * 高斯混合模型（GMM）
 * 使用EM算法进行参数估计
 */
public class GaussianMixtureModel {
    private int nComponents;           // 高斯成分数K
    private int dimension;             // 数据维度
    private GaussianComponent[] components; // 高斯成分数组
    private double[][] responsibilities; // 责任值 γ_ik
    private double logLikelihood;      // 对数似然
    private Random random;

    public GaussianMixtureModel(int nComponents, int dimension, long seed) {
        this.nComponents = nComponents;
        this.dimension = dimension;
        this.components = new GaussianComponent[nComponents];
        this.random = new Random(seed);

        for (int k = 0; k < nComponents; k++) {
            components[k] = new GaussianComponent(dimension);
        }
    }

    /**
     * 使用EM算法拟合数据
     */
    public void fit(double[][] data, int maxIterations, double tolerance) {
        int n = data.length;
        responsibilities = new double[n][nComponents];

        System.out.println("=== GMM-EM 训练开始 ===");
        System.out.println("样本数: " + n);
        System.out.println("数据维度: " + dimension);
        System.out.println("高斯成分数: " + nComponents);

        // 初始化参数
        initializeParameters(data);
        System.out.println("\n初始参数：");
        printParameters();

        double prevLogLikelihood = Double.NEGATIVE_INFINITY;

        for (int iter = 0; iter < maxIterations; iter++) {
            // E步：计算责任值
            expectationStep(data);

            // M步：更新参数
            maximizationStep(data);

            // 计算对数似然
            logLikelihood = computeLogLikelihood(data);
            double improvement = logLikelihood - prevLogLikelihood;

            if ((iter + 1) % 10 == 0 || iter < 5) {
                System.out.printf("迭代 %d: 对数似然 = %.4f, 提升 = %.6f%n",
                    iter + 1, logLikelihood, improvement);
            }

            // 检查收敛
            if (improvement < tolerance) {
                System.out.println("\n收敛！共迭代 " + (iter + 1) + " 次");
                break;
            }

            prevLogLikelihood = logLikelihood;
        }

        System.out.println("\n最终参数：");
        printParameters();
        System.out.printf("最终对数似然: %.4f%n", logLikelihood);
    }

    /**
     * 使用K-Means风格初始化
     */
    private void initializeParameters(double[][] data) {
        int n = data.length;

        // 随机选择K个中心点作为初始均值
        boolean[] selected = new boolean[n];
        for (int k = 0; k < nComponents; k++) {
            int idx;
            do {
                idx = random.nextInt(n);
            } while (selected[idx]);
            selected[idx] = true;

            components[k].setMean(data[idx].clone());

            // 初始化协方差为单位矩阵的缩放
            double[][] cov = new double[dimension][dimension];
            for (int d = 0; d < dimension; d++) {
                cov[d][d] = 1.0;
            }
            components[k].setCovariance(cov);
        }

        // 均匀初始化混合系数
        for (int k = 0; k < nComponents; k++) {
            components[k].setWeight(1.0 / nComponents);
        }
    }

    /**
     * E步：计算每个样本属于每个成分的后验概率
     */
    private void expectationStep(double[][] data) {
        int n = data.length;

        for (int i = 0; i < n; i++) {
            double[] weightedProbs = new double[nComponents];
            double sum = 0;

            for (int k = 0; k < nComponents; k++) {
                weightedProbs[k] = components[k].getWeight() *
                                   components[k].probabilityDensity(data[i]);
                sum += weightedProbs[k];
            }

            // 归一化得到责任值
            for (int k = 0; k < nComponents; k++) {
                responsibilities[i][k] = sum > 1e-10 ? weightedProbs[k] / sum : 0;
            }
        }
    }

    /**
     * M步：使用责任值更新模型参数
     */
    private void maximizationStep(double[][] data) {
        int n = data.length;

        for (int k = 0; k < nComponents; k++) {
            double Nk = 0;
            double[] newMean = new double[dimension];
            double[][] newCov = new double[dimension][dimension];

            // 累加责任值和加权样本
            for (int i = 0; i < n; i++) {
                double gamma = responsibilities[i][k];
                Nk += gamma;

                for (int d = 0; d < dimension; d++) {
                    newMean[d] += gamma * data[i][d];
                }
            }

            // 更新均值
            for (int d = 0; d < dimension; d++) {
                newMean[d] /= Nk;
            }
            components[k].setMean(newMean);

            // 更新协方差
            for (int i = 0; i < n; i++) {
                double gamma = responsibilities[i][k];
                for (int d1 = 0; d1 < dimension; d1++) {
                    for (int d2 = d1; d2 < dimension; d2++) {
                        double diff1 = data[i][d1] - newMean[d1];
                        double diff2 = data[i][d2] - newMean[d2];
                        newCov[d1][d2] += gamma * diff1 * diff2;
                    }
                }
            }

            for (int d1 = 0; d1 < dimension; d1++) {
                for (int d2 = d1; d2 < dimension; d2++) {
                    newCov[d1][d2] /= Nk;
                    newCov[d2][d1] = newCov[d1][d2];
                }
                // 对角线添加正则化
                newCov[d1][d1] += 1e-6;
            }
            components[k].setCovariance(newCov);

            // 更新混合系数
            components[k].setWeight(Nk / n);
        }
    }

    /**
     * 计算对数似然
     */
    private double computeLogLikelihood(double[][] data) {
        double logLik = 0;
        for (double[] point : data) {
            double prob = 0;
            for (int k = 0; k < nComponents; k++) {
                prob += components[k].getWeight() *
                        components[k].probabilityDensity(point);
            }
            logLik += Math.log(prob + 1e-10);
        }
        return logLik;
    }

    /**
     * 预测样本最可能属于的成分
     */
    public int predictCluster(double[] point) {
        double[] probs = new double[nComponents];
        for (int k = 0; k < nComponents; k++) {
            probs[k] = components[k].getWeight() *
                       components[k].probabilityDensity(point);
        }

        int bestCluster = 0;
        double maxProb = probs[0];
        for (int k = 1; k < nComponents; k++) {
            if (probs[k] > maxProb) {
                maxProb = probs[k];
                bestCluster = k;
            }
        }
        return bestCluster;
    }

    private void printParameters() {
        for (int k = 0; k < nComponents; k++) {
            System.out.println("成分 " + (k + 1) + ":");
            System.out.println(components[k]);
        }
    }

    public GaussianComponent[] getComponents() {
        return components;
    }
}

数据集生成器

import java.util.Random;

/**
 * 生成用于GMM演示的数据集
 */
public class DatasetGenerator {
    private Random random;

    public DatasetGenerator(long seed) {
        this.random = new Random(seed);
    }

    /**
     * 从多个高斯分布生成混合数据
     */
    public double[][] generateGMMData(int n, double[][] means, double[][][] covariances) {
        int k = means.length;
        int dim = means[0].length;
        double[][] data = new double[n][dim];

        int samplesPerCluster = n / k;
        for (int cluster = 0; cluster < k; cluster++) {
            int start = cluster * samplesPerCluster;
            int end = (cluster == k - 1) ? n : start + samplesPerCluster;

            for (int i = start; i < end; i++) {
                for (int d = 0; d < dim; d++) {
                    // 简化：使用独立高斯采样
                    double std = Math.sqrt(covariances[cluster][d][d]);
                    data[i][d] = means[cluster][d] + random.nextGaussian() * std;
                }
            }
        }

        return data;
    }

    /**
     * 生成三个分离的二维高斯簇
     */
    public double[][] generateThreeClusters(int n) {
        double[][] means = {
            {-3, -3},
            {3, 3},
            {0, 4}
        };
        double[][][] covs = {
            {{1, 0}, {0, 1}},
            {{1, 0}, {0, 1}},
            {{0.8, 0}, {0, 0.8}}
        };
        return generateGMMData(n, means, covs);
    }
}

演示入口

/**
 * EM算法与GMM演示
 */
public class EMDemo {
    public static void main(String[] args) {
        demoClustering();
    }

    /**
     * 演示：GMM聚类
     */
    private static void demoClustering() {
        System.out.println("=== 演示：高斯混合模型聚类 ===\n");

        DatasetGenerator generator = new DatasetGenerator(42);
        double[][] data = generator.generateThreeClusters(300);

        System.out.println("生成数据：3个高斯簇，共300个样本");
        System.out.println("数据维度: 2\n");

        // 使用GMM进行聚类
        GaussianMixtureModel gmm = new GaussianMixtureModel(3, 2, 42);
        gmm.fit(data, 100, 1e-6);

        // 统计每个簇的样本数
        int[] clusterCounts = new int[3];
        for (double[] point : data) {
            clusterCounts[gmm.predictCluster(point)]++;
        }

        System.out.println("\n聚类结果分布：");
        for (int i = 0; i < 3; i++) {
            System.out.println("  簇 " + (i + 1) + ": " + clusterCounts[i] + " 个样本");
        }
    }
}

复杂度分析

算法特性

优点

• 处理隐变量：天然适合含有隐变量的概率模型
• 保证收敛：每次迭代不降低对数似然
• 灵活性高：可应用于多种概率模型

缺点

• 可能陷入局部最优：对初始化敏感
• 计算成本高：需要多次迭代
• 需要预设成分数：GMM需要预先指定K值

算法应用

缺失值填补

EM算法可用于处理数据中的缺失值，将缺失值视为隐变量进行估计。

隐马尔可夫模型训练

Baum-Welch算法是EM算法在HMM中的特例，用于估计状态转移概率和发射概率。

主题模型

潜在狄利克雷分配（LDA）使用变分EM算法进行参数估计。

总结

EM算法通过E步和M步的交替执行，巧妙地解决了含隐变量的参数估计难题。本文以高斯混合模型为例，完整展示了责任值计算、参数更新的全过程。EM算法不仅是统计学习的基石，更是理解变分推断、EM-style深度学习等前沿方法的重要基础。