高斯混合模型(GMM)和EM算法 —— python实现

一、EM算法

EM算法是一种迭代算法,用于含有隐含变量的概率模型参数的极大似然估计。设Y为观测随机变量的数据,Z为隐藏的随机变量数据,Y和Z一起称为完全数据。

观测数据的似然函数为:

\[P(Y | \theta) = \sum_{Z} P(Y, Z | \theta) = \sum_{Z} P(Z | \theta)\ P(Y | Z, \theta) \]

模型参数θ的极大似然估计为:

\[\hat{\theta} = \underset{\theta}{\arg\max} \log P(Y | \theta) \]

这个问题只有通过迭代求解,下面给出EM算法的迭代求解过程:

step1、选择合适的参数初值θ(0),开始迭代

step2、E步,求期望。θ(i)为第i次迭代θ的估计值,在第i+1步,计算下面的Q函数:

\[\begin{aligned} Q(\theta, \theta^{(i)}) &= E_{Z}[\log P(Y, Z | \theta) | Y, \theta^{(i)}] \\ &= \sum_{Z} \log P(Y, Z | \theta) P(Z | Y, \theta^{(i)}) \end{aligned} \]

Q函数为logP(Y,Z|θ)关于在给定观测数据Y和当前参数θ(i)下对隐藏变量Z的条件概率分布P(Z|Y,θ(i))的期望。

step3、M步,求极大化。求使Q函数极大化的θ,确定第i+1次迭代的参数估计:

step4、重复第2、3步,直到收敛。

EM算法对初值的选取比较敏感,且不能保证找到全局最优解。

二、在高斯混合模型(GMM)中的应用

一维高斯混合模型:

\[P(x)=\sum_{k = 1}^{K} w_{k} \cdot \frac{1}{\sqrt{2\pi\sigma_{k}}} \exp\left(-\frac{(x - \mu_{k})^{2}}{2\sigma_{k}^{2}}\right) \]

多维高斯混合模型:

\[P(x)=\sum_{k = 1}^{K} w_{k} \cdot \frac{1}{(2\pi)^{d/2}|\Sigma|^{1/2}} \exp\left(-\frac{1}{2}(x - \mu_{k})^{T}\Sigma^{-1}(x - \mu_{k})\right) \]

其中,\(w_{k}\)\(k = 1, 2, \ldots, K\))为混合项系数,\(\sum_{k = 1}^{K} w_{k}=1\)\(\Sigma\)为协方差矩阵。\(\theta=(w_{k}, \mu_{k}, \sigma_{k})\)

设有N个可观测数据yi,它们是这样产生的:先根据概率wk选择第k个高斯分布模型,生成观测数据yi。yi是已知的,但yi属于第j个模型是未知的,是隐藏变量。用Zij表示隐藏变量,含义是第i个数据属于第j个模型的概率。先写出完全数据的似然函数,然后确定Q函数,要最大化期望,对wk、uk、σk求偏导并使其为0。可得高斯混合模型参数估计的EM算法(以高维数据为例):

step1、参数赋初始值,开始迭代

step2、E步,计算混合项系数Zij的期望 \(E[Z_{ij}]\)

\[E\left[Z_{i j}\right]=\frac{\left|\Sigma_{j}\right|^{-1 / 2} e^{-\frac{1}{2}\left(x^{(i)}-\mu_{j}\right)^{T} \Sigma_{j}^{-1}\left(x^{(i)}-\mu_{j}\right)}}{\sum_{k=1}^{K}\left|\Sigma_{k}\right|^{-1 / 2} e^{-\frac{1}{2}\left(x^{(i)}-\mu_{k}\right)^{T} \Sigma_{k}^{-1}\left(x^{(i)}-\mu_{k}\right)}} \]

step3、M步,计算新一轮迭代的参数模型:

\[\mu_j=\frac{\sum_{i=1}^NE[Z_{ij}]x^{(i)}}{\sum_{i=1}^NE[Z_{ij}]} \]

\[\sum_{j} = \frac{\sum_{i = 1}^{N} E[Z_{ij}] \cdot \{(x^{(i)} - \mu_{j})^T (x^{(i)} - \mu_{j})\}}{\sum_{i = 1}^{N} E[Z_{ij}]} \]

\[w_j = \frac{\sum_{i = 1}^{N} E[Z_{ij}]}{N} \]

step4、重复第2、3步,直到收敛。

前文搬运自链接:https://blog.csdn.net/u014157632/article/details/65442165

三、python程序示例

首先初始化分布。

MU1 = [1 2];
SIGMA1 = [1 0; 0 0.5];
MU2 = [-1 -1];
SIGMA2 = [1 0; 0 1];

并用这两个分布各生成1000个散点。

点击查看代码
import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture

# 初始化分布参数
MU1 = np.array([1, 2])
SIGMA1 = np.array([[1, 0], [0, 0.5]])
MU2 = np.array([-1, -1])
SIGMA2 = np.array([[1, 0], [0, 1]])

# 生成数据点
data1 = np.random.multivariate_normal(MU1, SIGMA1, 1000)
data2 = np.random.multivariate_normal(MU2, SIGMA2, 1000)

X = np.concatenate([data1, data2])  # [2000, 2]

#对X进行随机打乱,此步复现时不可忽略
np.random.shuffle(X)

接下来绘制二维散点图。

点击查看代码
plt.scatter(data1[:, 0], data1[:, 1], label='Class 1')
plt.scatter(data2[:, 0], data2[:, 1], label='Class 2')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Scatter Plot of Data Points')
plt.legend()
plt.show()

image

EM算法实现高斯混合模型的拟合

点击查看代码

import math


def fit_gaussian_mixture_model(X, n_components=2, max_iter=100):
    # 合并数据
    gmm = GaussianMixture(n_components=n_components, max_iter=max_iter)
    gmm.fit(X)

    return gmm.means_, gmm.covariances_, gmm.weights_


###拟合函数!
def my_fit_GMM(X, n_components=2, max_iter=1000):

    # 给数据排序
    n_samples, n_features = X.shape
    X_SortIndex = np.lexsort((X[:, 0], X[:, 1]))
    newArr = []
    for id in X_SortIndex:
        newArr.append(X[id])
    X = np.array(newArr)

    # 初始均值向量与协方差矩阵
    mu = X[:n_components, :]                  # [n_components, n_features]
    # cov = np.array(([[1.,0.],[0.,1.]],[[1.,0.],[0.,1.]]))
    cov = np.ones(shape=(n_components, n_features, n_features))        # [n_components, n_features, n_features]
    interval = n_samples // n_components

    for i in range(n_components):
        cov[i, ...] = np.cov(np.transpose(X[i*interval: (i+1)*interval, ...]))
        mu[i, ...] = np.mean(X[i*interval: (i+1)*interval, ...], axis=0)
    weights = np.ones(shape=(n_components, )) / n_components                      # [n_components, ]


    # 直接编写即可,不一定要使用结构体。
    expectation = np.zeros(shape=(n_samples, n_components))
    old_mu = np.copy(mu)

    X = np.matrix(X)
    mu = np.matrix(mu)
    expectation = np.matrix(expectation)

    for iteration in range(max_iter):
        old_mu = np.copy(mu)
        # 更新期望E
        for i in range(n_samples):
            for j in range(n_components):
                expectation[i, j] = math.exp(-(X[i,:]-mu[j,:])*np.linalg.inv(cov[j,...])*np.transpose(X[i,:]-mu[j,:])/2)\
                                    / np.sqrt(np.linalg.det(cov[j,...]))
                expectation[i, j] *= weights[j]

        t = np.sum(expectation, axis=1).reshape(-1, 1)
        expectation = expectation / t

        # 更新模型(均值向量)
        sum_E = np.zeros(shape=(n_components, ))
        for j in range(n_components):
            sum_E[j] = np.sum(expectation[:, j])
            mu[j, ...] = np.sum(np.multiply(expectation[:, j], X), axis=0) / sum_E[j]

        # 更新模型(协方差矩阵),混合项系数
        for j in range(n_components):
            cov_s = np.matrix(np.zeros((n_features, n_features)))
            for i in range(n_samples):
                cov_i = (X[i, ...] - mu[j, ...]).T * (X[i, ...] - mu[j, ...])
                cov_s = cov_s + np.multiply(expectation[i, j], cov_i)
            cov[j, ...] = np.array(cov_s) / sum_E[j]
            weights[j] = sum_E[j] / n_samples


        # 停止迭代的条件
        if np.abs(np.linalg.det(old_mu - mu)) < 1e-6:
            break

    means = np.array(mu)
    return means, cov, weights

# 请改为用my_fit_GMM拟合散点分布。并输出估计的均值,方差和混合项系数。

# means, cov, weights = fit_gaussian_mixture_model(X)
#
means, cov, weights = my_fit_GMM(X)
print("Data1.Means:")
print(means[0])
print("Data2.Means:")
print(means[1])

print("Data1.Covariances:")
print(cov[0])
print("Data2.Covariances:")
print(cov[1])

print("Weights:")
print(weights)

Data1.Means:
[-0.98994072 -0.99780233]
Data2.Means:
[1.00491397 1.98848228]
Data1.Covariances:
[[ 1.00762612 -0.04130926]
[-0.04130926 0.93607269]]
Data2.Covariances:
[[1.07037303 0.04512588]
[0.04512588 0.47164631]]
Weights:
[0.49222433 0.50777567]

点击查看代码
#绘制结果的概率密度函数,用于验证算法效果

from scipy.stats import multivariate_normal

def draw_results(m,s,w):

    # 定义两个二维高斯分布的参数
    m1 = m[0]
    s1 = s[0]

    m2 = m[1]
    s2 = s[1]

    # 生成网格点
    x, y = np.mgrid[-5:5:.01, -5:5:.01]
    pos = np.dstack((x, y))

    # 计算每个点的概率密度值
    rv1 = multivariate_normal(m1, s1)
    rv2 = multivariate_normal(m2, s2)
    z1 = rv1.pdf(pos)
    z2 = rv2.pdf(pos)

    # 混合两个高斯分布的概率密度函数
    z = w[0] * z1 + w[1] * z2  # 设置混合比例

    # 绘制3D图
    fig = plt.figure(figsize=(5, 10))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(x, y, z, cmap='viridis')

    # 设置坐标轴标签和标题
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Probability Density')
    ax.set_title('3D Plot of Mixture Gaussian Probability Density')

    plt.show()

draw_results(means, cov, weights)
#看到拟合的可视化结果了

image

posted @   小··明  阅读(627)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~
点击右上角即可分享
微信分享提示