em算法python代码实例

以下是一个简单的 EM 算法的 Python 代码实例，用于估计高斯混合模型的参数：

import numpy as np

# 初始化高斯混合模型参数
def init_params(data, k):
    n = data.shape[0]
    d = data.shape[1]
    # 初始化权重
    weights = np.ones(k) / k
    # 初始化均值
    means = np.random.randn(k, d)
    # 初始化协方差矩阵
    covs = np.zeros((k, d, d))
    for i in range(k):
        covs[i] = np.eye(d)
    return weights, means, covs

# E 步
def e_step(data, weights, means, covs):
    k = weights.shape[0]
    n = data.shape[0]
    # 初始化后验概率矩阵
    gamma = np.zeros((n, k))
    # 计算后验概率
    for i in range(n):
        for j in range(k):
            gamma[i, j] = weights[j] * normal_pdf(data[i], means[j], covs[j])
        gamma[i] /= np.sum(gamma[i])
    return gamma

# M 步
def m_step(data, gamma):
    k = gamma.shape[1]
    n = data.shape[0]
    d = data.shape[1]
    # 更新权重
    weights = np.sum(gamma, axis=0) / n
    # 更新均值
    means = np.zeros((k, d))
    for j in range(k):
        for i in range(n):
            means[j] += gamma[i, j] * data[i]
        means[j] /= np.sum(gamma[:, j])
    # 更新协方差矩阵
    covs = np.zeros((k, d, d))
    for j in range(k):
        for i in range(n):
            diff = data[i] - means[j]
            covs[j] += gamma[i, j] * np.outer(diff, diff)
        covs[j] /= np.sum(gamma[:, j])
    return weights, means, covs

# 计算多元高斯分布密度函数值
def normal_pdf(x, mean, cov):
    d = x.shape[0]
    coeff = 1.0 / np.sqrt((2*np.pi)**d * np.linalg.det(cov))
    diff = x - mean
    exponent = -0.5 * np.dot(np.dot(diff.T, np.linalg.inv(cov)), diff)
    return coeff * np.exp(exponent)

# EM 算法
def em_algorithm(data, k, max_iter):
    weights, means, covs = init_params(data, k)
    for i in range(max_iter):
        gamma = e_step(data, weights, means, covs)
        weights, means, covs = m_step(data, gamma)
    return weights, means, covs

其中，`data` 是输入数据，`k` 是高斯混合模型的个数，`max_iter` 是最大迭代次数。在 `init_params` 函数中，我们通过随机初始化来初始化高斯混合模型的参数。在 `e_step` 函数中，我们计算后验概率矩阵。在 `m_step` 函数中，我们使用后验概率矩阵来更新高斯混合模型的参数。在 `normal_pdf` 函数中，我们计算多元高斯分布密度函数值。最后，在 `em_algorithm` 函数中，我们使用 EM 算法来估计高斯混合模型的参数。