24/11/14 算法笔记 EM算法期望最大化算法

devtools/2024/11/17 0:48:19/

EM算法用于含有隐变量的概率参数模型的最大似然估计或极大后验概率估计。它在机器学习和统计学中有着广泛的应用,比如在高斯混合模型(GMM)、隐马尔可夫模型(HMM)以及各种聚类和分类问题中。

EM算法的基本思想是:首先根据已经给出的观测数据,估计出模型参数的值;然后再根据上一步估计出的参数值估计缺失数据的值,再根据估计出的缺失数据加上之前已经观测到的数据重新再对参数值进行估计,然后反复迭代,直至最后收敛,迭代结束

EM算法的迭代过程分为两个步骤:

  1. E步(Expectation Step):在这一步中,算法会估计缺失的或隐藏的变量,即基于当前模型参数的估计来计算隐藏变量的期望值。这一步实际上是在计算一个下界(lower bound)或者说是对缺失数据的期望可能性(expected likelihood)的估计。
  2. M步(Maximization Step):在这一步中,算法会优化模型参数以最大化在E步中计算出的期望可能性。这一步实际上是在最大化一个代理函数(surrogate function),这个函数依赖于隐藏变量的期望值。

EM算法的收敛性通常是有保证的,但收敛速度可能是一个问题,特别是在高维数据和复杂模型中。此外,EM算法能够估计复杂模型的参数,但这种复杂性可能会导致模型解释性降低。在实际应用中,我们需要仔细考虑这种权衡。

下面是一个简单的EM算法实现的例子,用于高斯混合模型(GMM)的参数估计

1. 首先,我们需要导入一些必要的库:

import numpy as np
from scipy.stats import multivariate_normal
  • scipy.stats 是SciPy库中的一个子库,它提供了统计学和概率论的函数,这里我们使用其中的 multivariate_normal 类来表示多变量高斯分布。

2. 然后,我们定义一个类来表示高斯混合模型:

class GaussianMixture:def __init__(self, n_components=3, covariance_type='full'):self.n_components = n_componentsself.covariance_type = covariance_typeself.weights = Noneself.means = Noneself.covariances = Nonedef fit(self, X, n_iter=100):n_samples, n_features = X.shapeself.weights = np.ones(self.n_components) / self.n_componentsself.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 responsibilitiesresponsibilities = self._compute_responsibilities(X)# M-step: Update parametersself._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 responsibilitiesdef _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 = Noneself.means = Noneself.covariances = None

2.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 类训练模型。
  • 预测类别:使用训练好的模型预测数据点的类别。

http://www.ppmy.cn/devtools/134577.html

相关文章

如何利用静态住宅IP提升TikTok营销效果:应对平台限制与账号安全的新利器

随着TikTok的全球化和日益严苛的运营政策,越来越多的品牌和商家开始面临着平台地域限制、账户管理及内容推广的挑战。特别是在多个账户管理和跨境营销中,如何避免账号封禁、提高内容曝光,成为了亟待解决的问题。在这种背景下,代理…

STM32 ADC --- DMA乒乓缓存

STM32 ADC — DMA乒乓缓存 文章目录 STM32 ADC --- DMA乒乓缓存软件切换实现乒乓利用DMA双缓冲实现乒乓 通过cubeMX配置生成HAL工程这里使用的是上篇文章(STM32 ADC — DMA采样)中生成的工程配置 软件切换实现乒乓 cubeMX默认生成的工程中是打开DMA中断…

【python系列】python数据类型之数字类型

1.定义 数字类型是编程中最常用的数据类型。什么是数字类型,下面是数字类型官方文档的解释:https://docs.python.org/zh-cn/3.10/library/stdtypes.html?highlightstr%20join#numeric-types-int-float-complex 以上可以知道: 数字类型包…

回调函数的概念、意义和应用场景

概念 回调函数,就是使用者自己定义一个函数,并实现函数的内容,然后把这个函数作为参数传入其它函数中,由其它函数在运行时来调用。 换句话说,函数是你实现的,但由别人的函数在运行时通过参数传递的方式调用…

flutter下拉刷新上拉加载的简单实现方式三

使用 CustomScrollView 结合 SliverList 实现了一个支持下拉刷新和上拉加载更多功能的滚动列表,对下面代码进行解析学习。 import dart:math;import package:flutter/material.dart;import custom_pull/gsy_refresh_sliver.dart; import package:flutter/cupertino…

hhdb数据库介绍(9-14)

SQL语法支持 DML语句 在关系集群数据库中,DML语句的逻辑将变的更为复杂。计算节点将DML语句分为两大类:单库DML语句与跨库DML语句。 单库DML语句,指SQL语句只需在一个节点上运行,即可计算出正确结果。假设分片表customer分片字…

杨中科 .Net Core 笔记 DI 依赖注入2

ServiceCollection services new ServiceCollection();//定义一个承放服务的集合 services.AddScoped<iGetRole, GetRole>();using (ServiceProvider serviceProvider services.BuildServiceProvider()) {var list serviceProvider.GetServices(typeof(iGetRole));//获…

力扣题目解析--合并两个链表

题目 将两个升序链表合并为一个新的 升序 链表并返回。新链表是通过拼接给定的两个链表的所有节点组成的。 示例 1&#xff1a; 输入&#xff1a;l1 [1,2,4], l2 [1,3,4] 输出&#xff1a;[1,1,2,3,4,4]示例 2&#xff1a; 输入&#xff1a;l1 [], l2 [] 输出&#xff…