项目简介

高斯混合模型（GMM, Gaussian Mixture Model） 是一种概率模型，用于表示一个具有多个高斯分布成分的整体数据分布。它是一种非常常用的聚类和密度估计方法，在许多应用场景中，如图像处理、语音识别、异常检测等领域都有广泛的应用。

GMM的基本思想是通过多个高斯分布的加权和来描述整个数据集的分布。每个高斯分布有自己的均值（mean）、方差（variance）和权重（weight），这些参数需要通过最大化似然函数来估计。

项目目标

本项目的目标是实现一个简化版本的高斯混合模型（GMM）。我们将手动实现其核心算法：期望最大化（EM, Expectation-Maximization）算法，并使用C++代码对给定数据进行聚类。

GMM算法概述

GMM算法主要包括两个步骤：E步（期望步骤）和M步（最大化步骤）。

1. E步：计算每个数据点属于每个高斯分布的概率（也叫“责任”）。
2. M步：通过最大化期望对数似然函数来更新高斯分布的参数（均值、方差、权重）。

这些步骤交替进行，直到收敛为止。

算法步骤

1. 初始化：随机初始化高斯分布的参数（均值、方差、权重）。
2. E步：计算每个数据点属于每个高斯分布的概率（责任）。
3. M步：根据E步计算出的责任，更新高斯分布的参数。
4. 重复E步和M步，直到模型收敛。

C++实现代码

下面是一个简化版本的GMM实现，使用**期望最大化（EM）**算法。

#include <iostream>
#include <vector>
#include <cmath>
#include <random>
#include <limits>

using namespace std;

// 定义一个二维数据点结构
struct DataPoint {
    double x, y;
    DataPoint(double _x, double _y) : x(_x), y(_y) {}
};

// 高斯分布类，用于描述每个高斯分布的参数
class Gaussian {
public:
    double mean_x, mean_y;    // 均值
    double variance_x, variance_y; // 方差
    double weight;             // 权重

    Gaussian(double mx, double my, double vx, double vy, double w)
        : mean_x(mx), mean_y(my), variance_x(vx), variance_y(vy), weight(w) {}

    // 计算给定点属于当前高斯分布的概率密度函数值
    double pdf(const DataPoint& p) const {
        double norm_x = exp(-(p.x - mean_x) * (p.x - mean_x) / (2 * variance_x)) / sqrt(2 * M_PI * variance_x);
        double norm_y = exp(-(p.y - mean_y) * (p.y - mean_y) / (2 * variance_y)) / sqrt(2 * M_PI * variance_y);
        return weight * norm_x * norm_y;
    }
};

// 高斯混合模型（GMM）类
class GMM {
private:
    vector<Gaussian> gaussians;    // 存储各个高斯分布
    int k;                         // 高斯分布的数量
    vector<vector<double>> responsibilities; // 存储责任矩阵
    vector<DataPoint> data;        // 输入数据

public:
    GMM(int k) : k(k) {}

    // 随机初始化高斯分布的参数
    void initialize(const vector<DataPoint>& dataset) {
        data = dataset;
        gaussians.clear();

        // 随机初始化均值、方差和权重
        random_device rd;
        mt19937 gen(rd());
        uniform_real_distribution<> dis(0.0, 1.0);

        // 均值初始化
        for (int i = 0; i < k; ++i) {
            double mean_x = dis(gen) * 10;
            double mean_y = dis(gen) * 10;
            double variance_x = dis(gen) + 0.1;  // 保证方差不为零
            double variance_y = dis(gen) + 0.1;
            double weight = 1.0 / k;
            gaussians.push_back(Gaussian(mean_x, mean_y, variance_x, variance_y, weight));
        }
    }

    // E步：计算每个数据点对每个高斯分布的责任
    void expectation() {
        responsibilities.clear();
        responsibilities.resize(data.size(), vector<double>(k, 0.0));

        // 计算每个数据点的责任（属于每个高斯分布的概率）
        for (size_t i = 0; i < data.size(); ++i) {
            double total_prob = 0.0;
            for (int j = 0; j < k; ++j) {
                total_prob += gaussians[j].pdf(data[i]);
            }
            for (int j = 0; j < k; ++j) {
                responsibilities[i][j] = gaussians[j].pdf(data[i]) / total_prob;
            }
        }
    }

    // M步：根据E步的责任，更新高斯分布的参数
    void maximization() {
        for (int j = 0; j < k; ++j) {
            double sum_responsibility = 0.0;
            double sum_x = 0.0, sum_y = 0.0;
            double sum_xx = 0.0, sum_yy = 0.0;

            // 计算当前高斯分布的参数
            for (size_t i = 0; i < data.size(); ++i) {
                double resp = responsibilities[i][j];
                sum_responsibility += resp;
                sum_x += resp * data[i].x;
                sum_y += resp * data[i].y;
                sum_xx += resp * data[i].x * data[i].x;
                sum_yy += resp * data[i].y * data[i].y;
            }

            // 更新均值、方差和权重
            gaussians[j].mean_x = sum_x / sum_responsibility;
            gaussians[j].mean_y = sum_y / sum_responsibility;
            gaussians[j].variance_x = sum_xx / sum_responsibility - gaussians[j].mean_x * gaussians[j].mean_x;
            gaussians[j].variance_y = sum_yy / sum_responsibility - gaussians[j].mean_y * gaussians[j].mean_y;
            gaussians[j].weight = sum_responsibility / data.size();
        }
    }

    // 训练GMM
    void train(const vector<DataPoint>& dataset, int max_iterations = 100, double tol = 1e-6) {
        initialize(dataset);

        for (int iter = 0; iter < max_iterations; ++iter) {
            expectation();
            maximization();

            // 计算对数似然估计（检查收敛）
            double log_likelihood = 0.0;
            for (size_t i = 0; i < data.size(); ++i) {
                double total_prob = 0.0;
                for (int j = 0; j < k; ++j) {
                    total_prob += gaussians[j].pdf(data[i]);
                }
                log_likelihood += log(total_prob);
            }

            // 如果变化小于tol，说明收敛，停止迭代
            if (iter > 0 && abs(log_likelihood - prev_log_likelihood) < tol) {
                break;
            }
            prev_log_likelihood = log_likelihood;
        }
    }

    // 打印模型参数
    void printModel() {
        for (int i = 0; i < k; ++i) {
            cout << "Gaussian " << i + 1 << " parameters:" << endl;
            cout << "Mean: (" << gaussians[i].mean_x << ", " << gaussians[i].mean_y << ")" << endl;
            cout << "Variance: (" << gaussians[i].variance_x << ", " << gaussians[i].variance_y << ")" << endl;
            cout << "Weight: " << gaussians[i].weight << endl;
        }
    }

private:
    double prev_log_likelihood = -numeric_limits<double>::infinity();
};

int main() {
    // 创建一个包含二维数据的样本集
    vector<DataPoint> dataset = {
        DataPoint(1.0, 2.0), DataPoint(1.2, 2.1), DataPoint(3.0, 4.0),
        DataPoint(3.2, 4.1), DataPoint(10.0, 10.0), DataPoint(10.2, 10.1)
    };

    // 创建一个GMM对象，并训练模型
    GMM gmm(2);  // 假设我们要拟合两个高斯分布
    gmm.train(dataset);

    // 打印训练后的模型参数
    gmm.printModel();

    return 0;
}

代码解读

1. DataPoint 结构：用于表示二维数据点。
2. Gaussian 类：每个高斯分布都有均值（mean_x 和 mean_y），方差（variance_x 和 variance_y），以及该分布的权重（weight）。pdf 方法计算给定数据点在该高斯分布下的概率密度。
3. GMM 类：
   • initialize：初始化高斯分布的参数（均值、方差、权重）。
   • expectation：计算每个数据点对每个高斯分布的责任。
   • maximization：根据责任矩阵更新每个高斯分布的参数。
   • train：训练模型，通过反复执行E步和M步来逼近最优解。
   • printModel：打印每个高斯分布的参数（均值、方差、权重）。
4. 主函数：创建一个数据集，训练GMM模型并输出结果。

项目总结

1. 高斯混合模型（GMM） 是一种强大的概率模型，用于聚类和密度估计。它通过多个高斯分布的加权和描述数据的总体分布。
2. EM算法 是解决GMM参数估计问题的常用方法，通过交替执行E步（计算责任）和M步（更新参数）来优化模型。
3. 应用场景：GMM可以用于聚类、异常检测、密度估计等任务，广泛应用于图像处理、语音识别、金融建模等领域。

通过本项目的实现，我们深入了解了高斯混合模型的核心思想，并掌握了如何使用C++实现GMM算法。