机器学习之极大期望 EM,估算多维高斯混合模型参数（C语言描述）

统计学习7

第七章 极大期望估计 EM 算法

---

---

前言

最近拜读了机器学习领域，经典书籍，李航老师的《统计学习方法》，深受裨益。现实现算法，记录分享。本章，还结合就另外一本比较好的统计学习的书 《机器学习经典算法–基于opcv》里面关于极大期望 EM 算法。

---

一、什么是 EM 算法

EM 算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或者极大后验概率估计。算法由两步迭代组成：E 步，求期望（expectation）；M 步，求极大（maximization）。所以这一算法称作极大期望算法，简称 EM 算法。

概率模型有时候既有观测变量，又有隐变量。如果这个概率模型的变量都是观测变量的，那么给定数据，可以直接使用极大似然估计法去估计模型的参数。什么叫有隐含变量？常见的例子是，有一堆数据，但是不知道这堆数据的分类。不知道这堆数据的就不知道它们的参数。例如前面所说一个正太贝叶斯分类器： y ^ = a r g m a x k = { 1 , 2.. , k } p ( x 1 , x 2 . . . , x n ∣ C k ) = 1 ( 2 π ) n ∣ Σ k ∣ e x p ( − 1 2 ( x − μ k ) T Σ k − 1 ( x − μ k ) ) \begin{aligned} \widehat{y} & = \mathop{argmax}\limits_{k=\{1,2..,k\}}p(x_1,x_2...,x_n|C_k) \\ & = \frac{1}{\sqrt{(2\pi)^n|\Sigma_k|} } exp(-\frac{1}{2}(\boldsymbol{x-\mu_k})^T\Sigma_k^{-1}(\boldsymbol{x-\mu_k})) \end{aligned} y ​​=k={1,2..,k}argmax​p(x1​,x2​...,xn​∣Ck​)=(2π)n∣Σk​∣ ​1​exp(−21​(x−μk​)TΣk−1​(x−μk​))​要使上述的分类器能工作，那必须知道 $C_k$ 也就是每个数据分类数据。EM 算法就是在没有分类数据的情况下，估算这里面的参数：$\mu ,\Sigma$。

二、EM 算法分析t

推理过程免了，直接上算法过程，现有一堆数据，$X$ 可观察数据，$Z$ 为隐藏数据，$\theta$ 为概率分布的参数， $X$ 与 $Z$ 合在一起成为完整数据。现在设 $P(X,Z|\theta )$, 为联合分布的似然函数，$P(Z|X,\theta )$ 为 Z 在 关于 X 与 $\theta$ 条件下的概率分布。然后：

0、给出 ${\theta }^0$ 的初始化。

1、给出完全数据的对数似然函数 $logP(X,Z|\theta )$ 关于在给定观察数据 X 和 当前 $\theta$ 下对未知数据 Z 的条件概率分布 $P(Z|X,\theta )$ 的期望定义，即 E 步：$$Q({\theta }^i)=E_z(logP(X,Z|{\theta }^i)|X,{\theta }^i)$$
2、 求使得 ${\theta }^i$ 最大化的 参数 ${\theta }^{i+1}$，完成迭代 ${\theta }^i\to {\theta }^{i+1}$，即 M 步：$${\theta }^{i+1}=\underset{\theta }{\mathrm{argmax}}Q({\theta }^i)$$3、重复 1、2 步骤，直到 $\theta$ 收敛。

EM 算法在应用上通常配合高斯混合模型。EM 算法可以估计出，高斯混合模型的模型比例 $\alpha$、各个模型的均值 $\mu$ 以及方差。若是多维高斯混合模型，则可以估计出模型的均值向量 $\mathbf{\mu }$ 和 协方差矩阵 $\Sigma$。下面是 K 个多维高斯混合模型：$$\begin{array}{rl}p(X|\theta )& =\sum _{k=1}^K{\alpha }_{k}N(X|{\mu }_k,{\Sigma }_k)\\ & =\sum _{k=1}^K{\alpha }_k\frac{1}{(2\pi {)}^{d/2}|{\Sigma }_k{|}^{1/2}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }_k^{-1}(x-\mu )}\end{array}$$结合上述 EM 算法得出：

1、E 步：$$E_{z_k}({\alpha }^{(t-1)},{\mu }^{(t-1)},{\Sigma }^{(t-1)})=\frac{{\alpha }_k^{(t-1)}N(x_i|{\mu }_k^{(t-1)},{\Sigma }_k^{(t-1)})}{{\sum }_{k=1}^K{\alpha }_k^{(t-1)}N(x_i|{\mu }_k^{(t-1)},{\Sigma }_k^{(t-1)})}$$2、M步：$$\begin{array}{rl}& {\alpha }_k^{(t)}=\frac{1}{n}\sum _{i=1}^n{E}_{z_k}^{(i)}\\ & \\ & {\mu }_k^{(t)}=\frac{{\sum }_{i=1}^n[x_i{E}_{z_k}^{(i)}]}{{\sum }_{i=1}^n{E}_{z_k}^{(i)}}\\ & \\ & {\Sigma }_k^{(t)}=\frac{{\sum }_{i=1}^n[E_{z_k}^{(i)}(x_i-{\mu }^t)(x_i-{\mu }^t{)}^T]}{{\sum }_{i=1}^n}\end{array}$$这里说明一下：$E_{z_k}$ 是一个 n * 1 维向量，n等于数据集合 X 的行数，

三、上代码

/*
 * @Author: zuweie jojoe.wei@gmail.com
 * @Date: 2023-08-15 14:48:47
 * @LastEditors: zuweie jojoe.wei@gmail.com
 * @LastEditTime: 2023-09-26 08:42:34
 * @FilePath: /boring-code/src/statistical_learning/em.c
 * @Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
 */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <float.h>
#include "em.h"

static int __init_apK_data(matrix2_t* apK)
{
    //srand((unsigned) time (NULL));

    for (int i=0; i<apK->rows * apK->cols; ++i) {
        apK->pool[i] = (double)rand() / RAND_MAX;
    }
    return 0;
}

static int __init_mu_data(matrix2_t* mu) 
{
    //Mat2_fill(mu, 0.f);
    srand((unsigned) time (NULL));

    for (int i=0; i<mu->rows * mu->cols; ++i) {
        mu->pool[i] = 0.5;//(double) rand() / RAND_MAX;
    }
    return 0;
}

/**
 * @brief 
 * 
 * @param sigma_data 
 * @param size 
 * @return int 
 */
static int __init_sigma_data(matrix2_t* sigma) 
{
    //memset(sigma_data, 0, size*size*sizeof(vfloat_t));

    Mat2_fill(sigma, 0.f);

    MAT2_POOL_PTR(sigma, sigma_ptr);

    srand((unsigned) time (NULL));

    for (int i=0; i<sigma->rows; ++i) {
        sigma_ptr[i][i] = 0.99f;//(double) rand() / RAND_MAX;
    }
    return 0;
}
/**
 * @brief 书本中 9-82 式, 计算第 K 模型下第 i 条数据的后验证概率。
 * 
 * @param _X 
 * @param K 
 * @param mu 
 * @param sigema 
 * @param alphaK 
 * @return int 
 */
static double __calculate_apK(matrix2_t* _Xi, matrix2_t* sigma, double alpha)
{
    double det_sigma = 1.f;
    double pK        = 0.f;

    matrix2_t* w_mat         = Mat2_create(1,1);
    matrix2_t* u_mat         = Mat2_create(1,1);
    matrix2_t* _Xi_T         = Mat2_create(1,1);
    matrix2_t* inverse_sigma = Mat2_create_cpy(sigma);
    
    Mat2_inverse(inverse_sigma);

    Mat2_eig(w_mat, u_mat, inverse_sigma);

    Mat2_dot(_Xi, u_mat);

    Mat2_cpy(_Xi_T, _Xi);

    Mat2_T(_Xi_T);

    Mat2_hadamard_product(_Xi, w_mat);
    
    Mat2_dot(_Xi, _Xi_T);

    for (int k=0; k<w_mat->cols; ++k) {
        det_sigma *= 1.f / w_mat->pool[k];
    }

    pK += log(alpha);
    pK -= 0.5 * log(det_sigma);
    pK -= 0.5 * _Xi->pool[0];

    Mat2_destroy(w_mat);
    Mat2_destroy(u_mat);
    Mat2_destroy(_Xi_T);
    Mat2_destroy(inverse_sigma);
    
    double ret = exp(pK);
    //printf("ret: %lf \n", ret);
    return ret;
}

/**
 * @brief 使用原始的算法试试。
 * 
 * @param _Xi 
 * @param sigma 
 * @param alpha 
 */
static double __calculate_apK_2(matrix2_t* _Xi, matrix2_t* mu,  matrix2_t* sigma, double alpha) 
{
    // 此处的 _Xi mu, sigma 对象不能被改变，因为在外面还要用到 
    // matrix2_t* _Xi_cpy = Mat2_create_cpy(_Xi); 
    matrix2_t* _Xi_T   = Mat2_create(1,1);
    // matrix2_t* mu_cpy  = Mat2_create_cpy(mu);
    // matrix2_t* sigma_cpy = Mat2_create_cpy(sigma);

    int d = _Xi->cols;

    vfloat_t det_sigma;
    Mat2_det(sigma, &det_sigma);
    Mat2_inverse(sigma);

    Mat2_T(mu);
    Mat2_sub(_Xi, mu);
    
    Mat2_cpy(_Xi_T, _Xi);
    Mat2_T(_Xi_T);

    Mat2_dot(_Xi, sigma);


    Mat2_dot(_Xi, _Xi_T);

    double pk = 0.f;
    pk += log(alpha);
    pk -= 0.5f * d * log(2*3.1415926);
    pk -= 0.5f * log(det_sigma);
    pk -= 0.5f * _Xi->pool[0];

    double ret = exp(pk);
    //printf("<alpha: %lf, det_sigm: %lf, xi[0]: %lf, pk: %lf, ret: %lf> \n", alpha, det_sigma, _Xi->pool[0], pk, ret);
    //printf("ret: %lf ", ret);

    Mat2_destroy(_Xi_T);
    // Mat2_destroy(mu_cpy);
    // Mat2_destroy(sigma_cpy);
    // Mat2_destroy(_Xi_cpy);


    return ret;
}


//  计算 期望 步。
/**
 * @brief 计算各个模型下，p(y|X) 的先验概率。
 * 
 * @param _X 
 * @param K 
 * @param alphas 
 * @param sigmas 
 * @param apK 
 * @return int 
 */
static int __e_step(matrix2_t* _X, int K, double* alphas, matrix2_t** mus, matrix2_t** sigmas, matrix2_t** apKs)
{
    matrix2_t* _Xi   = Mat2_create(1,1);
    matrix2_t* sigma = Mat2_create(1,1);
    matrix2_t* mu    = Mat2_create(1,1);
    matrix2_t* apK;
    double alpha;

    for (int i=0; i<_X->rows; ++i) {

       
        double p = 0.f;

        for (int j=0; j<K; ++j) {

            Mat2_slice_row_to(_Xi, _X, i);

            alpha  = alphas[j];
            apK    = apKs[j];

            Mat2_cpy(mu, mus[j]);
            Mat2_cpy(sigma, sigmas[j]);

            //apK->pool[i] = __calculate_apK(_Xi, sigma, alpha);
            apK->pool[i] = __calculate_apK_2(_Xi, mu, sigma, alpha);
            p += apK->pool[i];
        }

        // 归一化。
        for (int l=0; l<K; ++l) {
            apK = apKs[l];
            apK->pool[i] /= p;
        }
    }

    Mat2_destroy(_Xi);
    Mat2_destroy(mu);
    Mat2_destroy(sigma);
    return 0;
}

// 计算 M 步 
/**
 * @brief 计算 M 步，即更新 alpha、mu、sigma
 * 
 * @param _X 
 * @param K 
 * @param alphas 
 * @param mus 
 * @param sigmas 
 * @param apK 
 * @return int 
 */
static int __m_step(matrix2_t* _X, int K, double* alphas, matrix2_t** mus, matrix2_t** sigmas, matrix2_t** apKs)
{
    matrix2_t* apK_cpy = Mat2_create(1,1);
    matrix2_t* _Xi     = Mat2_create(1,1);
    matrix2_t* _Xi_T   = Mat2_create(1,1);
    double sum_apk     = 0.f;
    for (int i=0; i<K; ++i) {

        matrix2_t* apK   = apKs[i];
        matrix2_t* mu    = mus[i];
        matrix2_t* sigma = sigmas[i];
        
        Mat2_cpy(apK_cpy, apKs[i]);
        //int n = apK_cpy->rows;
        Mat2_sum(apK_cpy, 0);
        sum_apk = apK_cpy->pool[0];
        
        // TODO: 1 各个模型的 alpah 
        alphas[i] =  sum_apk / apK->rows;
        // end 1




        // TODO： 3 各个模型中 sigma 计算 sigma
        Mat2_fill(sigma, 0.f);

        for (int r=0; r<_X->rows; ++r) {

            Mat2_slice_row_to(_Xi, _X, r);

            Mat2_T(_Xi);

            Mat2_sub(_Xi, mu);

            Mat2_cpy(_Xi_T, _Xi);
            
            Mat2_T(_Xi_T);

            Mat2_dot(_Xi, _Xi_T);

            Mat2_scalar_multiply(_Xi, apK->pool[r]);
            
            Mat2_add(sigma, _Xi);
        }

        Mat2_scalar_multiply(sigma, 1 / sum_apk);
        // end 3

        
        // TODO： 2 各个模型中的 mu 
        Mat2_cpy(mu, _X);

        MAT2_POOL_PTR(mu, mu_ptr);

        for (int m=0; m<mu->rows; ++m) {
            for (int l=0; l<mu->cols; ++l) {
                mu_ptr[m][l] *= apK->pool[m];
            }
        }
        
        Mat2_sum(mu, 0);
        
        Mat2_scalar_multiply(mu,  1 / sum_apk);

        Mat2_T(mu);

        // end 2
    
    }

    Mat2_destroy(apK_cpy); 
    Mat2_destroy(_Xi);    
    Mat2_destroy(_Xi_T); 

    return 0;
}

// EM 的训练。

/**
 * @brief 通过 EM 算法， 计算 K 个 alpha，K 个 mu，K 个 sigema
 * 
 * @param _X D 维向量输入数据。
 * @param K 高斯模型的个数，人工输入。
 * @param eps EM 训练的终止条件。
 * @param alphas 输出的 alpha，有 K 个。
 * @param mus 输出的 mu 向量，有 K 个，每个 D 维
 * @param sigemas 
 * @return int 
 */
int EM_train(matrix2_t* _X, int K, int Max_iter, double epsilon, double** alphas, matrix2_t** mus, matrix2_t** sigmas, void (*progress)(const char*, double, unsigned long, unsigned long))
{
    // TODO : 非常简单，就是不断的进行 E 步 与 M 的轮流计算，获取高斯分布的那几个关键的参数。nu sigma，与 alpha。
    matrix2_t** __apKs     = (matrix2_t**)malloc (K * sizeof(matrix2_t*));
    matrix2_t** __mus      = (matrix2_t**)malloc (K * sizeof(matrix2_t*));
    matrix2_t** __sigmas   = (matrix2_t**) malloc (K * sizeof(matrix2_t*));
   

    double*    __alphas    = (double*) malloc (K * sizeof(double));

    double last_apks_norm[K];

    for (int i=0; i<K; ++i) {
        __mus[i]    = Mat2_create(_X->cols, 1);
        //__init_mu_data(__mus[i]);

        __sigmas[i] = Mat2_create(_X->cols, _X->cols);
        //__init_sigma_data(__sigmas[i]);

        //MAT2_INSPECT(__sigmas[i]);

        __apKs[i]   = Mat2_create(_X->rows, 1);
        __init_apK_data(__apKs[i]);

        __alphas[i] = 1.f / (double) K;     

        last_apks_norm[i] = FLT_MAX; 
    }

    int iter = 0;
    double sum_diff = FLT_MAX;
    double curr_norm = 0.f;
    double last_sum = 0.f;

    // 用于记录每次迭代后的 apk 的模长度，当长度不在变化，那么表明 apk 不再有变化于是收敛了。


    while ( ++iter <= Max_iter && sum_diff  > epsilon) {

        if (progress)
            progress("M step ... ", 0.f, iter, Max_iter);

        __m_step(_X, K, __alphas, __mus, __sigmas, __apKs);

        if (progress)
            progress("E step ... ", 0.f, iter, Max_iter);

        __e_step(_X, K, __alphas, __mus, __sigmas, __apKs);

        // __apks 被更新后，即可计算它的模长度，看看跟上一轮有没有变化，有变化即可继续优化。
        sum_diff = 0.f;

        for (int i=0; i<K; ++i) {
            
            curr_norm = Mat2_norm(__apKs[i]);
            sum_diff += fabs( curr_norm - last_apks_norm[i] );
            last_apks_norm[i] = curr_norm;

        }
        if (progress)
            progress(" epsilon ... ", sum_diff, iter, Max_iter);
    }
    // 输出结果。
    *alphas = __alphas;
    *mus    = __mus;
    *sigmas = __sigmas;
    return 0;
}

int EM_recycle(int K, double** alphas, matrix2_t** mus, matrix2_t** sigmas)
{
    for (int i=0; i<K; ++i) {
        matrix2_t* m = mus[i];
        Mat2_destroy(m);

        matrix2_t* s = sigmas[i];
        Mat2_destroy(s);
    }

    free(alphas);
    free(mus);
    free(sigmas);
    return 0;
}

四、看结果

先使用 Python 生成两组数据:

第一组的 $\mu 和\Sigma$ 是，[5,3,7,8]、[25, 64, 49, 81]，共3000条数据。
第二组的 $\mu 和\Sigma$ 是，[2,1,8,5]、[121,64,25,81]，共9000条数据。
使用精度是：1e-5。

经过376秒的估算得到以下结果：

K0 是第二组数据，占比是 0.7574。估算的结果：$$\begin{array}{rl}& \mu =[2.09,0.85,7.89,5.109]\\ & \\ & \Sigma =\left[\begin{array}{cccc}120.17,& -1.06,& -1.48,& 0.79\\ -1.06,& 65.3,& 0.78,& 0.09\\ -1.48,& 0.78,& 25.31,& 0.95\\ 0.79,& 0.09,& 0.95,& 84.34\end{array}\right]\end{array}$$

K1 是第一组数据，占比是 0.24251。估算的结果：$$\begin{array}{rl}& \mu =[4.92,3.18,6.979,7.5728]\\ & \\ & \Sigma =\left[\begin{array}{cccc}23.51,& 0.33,& 0.51,& 1.43\\ 0.33,& 60.48,& -1.55,& 0.02\\ 0.51,& -1.55,& 47.42,& -1.64\\ 1.43,& 0.02,& -1.64,& 80.44\end{array}\right]\end{array}$$

源码：http://github.com/zuweie/boringcode/blob/main/src/statistical_learning/em.c

---

总结

李航老师中的《统计学习》给出的估算高斯混合模型的算法是一维的，也就是其中的高斯分布是一维的。本文根据使用的是多维高斯分布混合模型，更具有实用性。另外书中提到，在算法开始时对参数的初始化，会对最后的估算结果有影响，在本文的算法中，对 $\mu 与\Sigma$ 的初始化只是作出简单的随机值的赋值。得到的结果只与预测的结果有 90% 相似。至于怎么样对参数初始值进行初始化，会更有利于估算结果，书中并没有提及。

完