EM算法用于含有隐变量的概率参数模型的最大似然估计或极大后验概率估计。它在机器学习和统计学中有着广泛的应用，比如在高斯混合模型（GMM）、隐马尔可夫模型（HMM）以及各种聚类和分类问题中。
EM算法的基本思想是：首先根据已经给出的观测数据，估计出模型参数的值；然后再根据上一步估计出的参数值估计缺失数据的值，再根据估计出的缺失数据加上之前已经观测到的数据重新再对参数值进行估计，然后反复迭代，直至最后收敛，迭代结束
EM算法的迭代过程分为两个步骤：
- E步（Expectation Step）：在这一步中，算法会估计缺失的或隐藏的变量，即基于当前模型参数的估计来计算隐藏变量的期望值。这一步实际上是在计算一个下界（lower bound）或者说是对缺失数据的期望可能性（expected likelihood）的估计。
- M步（Maximization Step）：在这一步中，算法会优化模型参数以最大化在E步中计算出的期望可能性。这一步实际上是在最大化一个代理函数（surrogate function），这个函数依赖于隐藏变量的期望值。
EM算法的收敛性通常是有保证的，但收敛速度可能是一个问题，特别是在高维数据和复杂模型中。此外，EM算法能够估计复杂模型的参数，但这种复杂性可能会导致模型解释性降低。在实际应用中，我们需要仔细考虑这种权衡。
下面是一个简单的EM算法实现的例子，用于高斯混合模型（GMM）的参数估计
1. 首先，我们需要导入一些必要的库：
import numpy as np
from scipy.stats import multivariate_normalscipy.stats是SciPy库中的一个子库，它提供了统计学和概率论的函数，这里我们使用其中的multivariate_normal类来表示多变量高斯分布。
2. 然后，我们定义一个类来表示高斯混合模型：
class GaussianMixture:
def __init__(self, n_components=3, covariance_type='full'):
self.n_components = n_components
self.covariance_type = covariance_type
self.weights = None
self.means = None
self.covariances = None
def fit(self, X, n_iter=100):
n_samples, n_features = X.shape
self.weights = np.ones(self.n_components) / self.n_components
self.means = np.random.rand(self.n_components, n_features)
self.covariances = np.array([np.eye(n_features)] * self.n_components)
for i in range(n_iter):
# E-step: Compute responsibilities
responsibilities = self._compute_responsibilities(X)
# M-step: Update parameters
self._update_parameters(X, responsibilities)
def _compute_responsibilities(self, X):
responsibilities = np.zeros((X.shape[0], self.n_components))
for i in range(self.n_components):
rv = multivariate_normal(self.means[i], self.covariances[i])
responsibilities[:, i] = self.weights[i] * rv.pdf(X) / np.sum([self.weights[j] * multivariate_normal(self.means[j], self.covariances[j]).pdf(X) for j in range(self.n_components)], axis=0)
return responsibilities
def _update_parameters(self, X, responsibilities):
self.weights = np.mean(responsibilities, axis=0)
for i in range(self.n_components):
self.means[i] = np.sum(responsibilities[:, i].reshape(-1, 1) * X, axis=0) / np.sum(responsibilities[:, i])
self.covariances[i] = np.dot((X - self.means[i]).T * responsibilities[:, i], (X - self.means[i])) / np.sum(responsibilities[:, i])
def predict(self, X):
return np.argmax([self._compute_responsibilities(X) for i in range(self.n_components)], axis=1)这个简单的实现包括了：
- 初始化模型参数。
- E步：计算每个数据点属于每个高斯分量的“责任”（即概率）。
- M步：根据这些责任更新模型参数（权重、均值、协方差）。
我们来分析一下各段代码
2.1 定义高斯混合模型类
class GaussianMixture:
def __init__(self, n_components=3, covariance_type='full'):
self.n_components = n_components #混合模型中的高斯分量数量，默认为3。
self.covariance_type = covariance_type #指定协方差矩阵的类型，这里使用 'full' 表示每个分量都有自己的完整协方差矩阵。
#分别存储每个高斯分量的权重、均值和协方差矩阵，它们在模型训练过程中会被更新。
self.weights = None
self.means = None
self.covariances = None2.2 训练模型
def fit(self, X, n_iter=100):
n_samples, n_features = X.shape
#初始化模型的混合权重
self.weights = np.ones(self.n_components) / self.n_components
#随机初始化每个组件的均值
self.means = np.random.rand(self.n_components, n_features)
#初始化每个组件的协方差矩阵。np.eye(n_features) 生成一个单位矩阵，表示每个组件的初始协方差矩阵是单位矩阵。
self.covariances = np.array([np.eye(n_features)] * self.n_components)
for i in range(n_iter):
# E-step: Compute responsibilities
#计算责任度。根据当前的模型参数（均值、协方差和权重）来计算每个数据点属于每个组件的概率。
responsibilities = self._compute_responsibilities(X)
# M-step: Update parameters
#根据计算出的责任度来更新均值、协方差和权重，以最大化数据的似然函数。
self._update_parameters(X, responsibilities)2.3 E步：计算责任
def _compute_responsibilities(self, X):
responsibilities = np.zeros((X.shape[0], self.n_components))
for i in range(self.n_components):
rv = multivariate_normal(self.means[i], self.covariances[i])#高斯正态分布
responsibilities[:, i] = self.weights[i] * rv.pdf(X) / np.sum([self.weights[j] * multivariate_normal(self.means[j], self.covariances[j]).pdf(X) for j in range(self.n_components)], axis=0)
return responsibilities- 对于每个分量，使用
multivariate_normal.pdf方法计算数据点的密度，然后通过权重和归一化因子计算责任。
2.4 M步：更新参数
def _update_parameters(self, X, responsibilities):
self.weights = np.mean(responsibilities, axis=0)
for i in range(self.n_components):
self.means[i] = np.sum(responsibilities[:, i].reshape(-1, 1) * X, axis=0) / np.sum(responsibilities[:, i])
self.covariances[i] = np.dot((X - self.means[i]).T * responsibilities[:, i], (X - self.means[i])) / np.sum(responsibilities[:, i])- 更新权重：每个分量的权重是该分量责任的平均值。
- 更新均值：每个分量的均值是数据点的加权平均，权重是责任。
- 更新协方差：每个分量的协方差是加权的数据点偏差的外积，权重是责任。
2.5 预测数据点的类别
def predict(self, X):
return np.argmax([self._compute_responsibilities(X) for i in range(self.n_components)], axis=1)3 使用示例
# 生成一些模拟数据
np.random.seed(0)
data = np.concatenate([np.random.normal(0, 1, (100, 2)), np.random.normal(5, 2, (100, 2))])
data[:, 1] += data[:, 0] * 3
# 训练模型
gmm = GaussianMixture(n_components=2)
gmm.fit(data)
# 预测数据点的类别
labels = gmm.predict(data)
print(labels)- 生成模拟数据：两个高斯分布，分别在 (0, 1) 和 (5, 2) 附近。
- 训练模型：使用
GaussianMixture类训练模型。 - 预测类别：使用训练好的模型预测数据点的类别。
