一、EM算法
EM算法是一种迭代算法，用于含有隐含变量的概率模型参数的极大似然估计。设Y为观测随机变量的数据，Z为隐藏的随机变量数据，Y和Z一起称为完全数据。
观测数据的似然函数为：
模型参数θ的极大似然估计为：
这个问题只有通过迭代求解，下面给出EM算法的迭代求解过程：
step1、选择合适的参数初值θ(0)，开始迭代
step2、E步，求期望。θ(i)为第i次迭代θ的估计值，在第i+1步，计算下面的Q函数：
Q函数为logP(Y,Z|θ)关于在给定观测数据Y和当前参数θ(i)下对隐藏变量Z的条件概率分布P(Z|Y,θ(i))的期望。
step3、M步，求极大化。求使Q函数极大化的θ，确定第i+1次迭代的参数估计：
step4、重复第2、3步，直到收敛。
EM算法对初值的选取比较敏感，且不能保证找到全局最优解。
二、在高斯混合模型（GMM）中的应用
一维高斯混合模型：
多维高斯混合模型：
wk（k=1，2，……，K）为混合项系数，和为1。∑为协方差矩阵。θ=（wk，uk，σk）。
设有N个可观测数据yi，它们是这样产生的：先根据概率wk选择第k个高斯分布模型，生成观测数据yi。yi是已知的，但yi属于第j个模型是未知的，是隐藏变量。用Zij表示隐藏变量，含义是第i个数据属于第j个模型的概率。先写出完全数据的似然函数，然后确定Q函数，要最大化期望，对wk、uk、σk求偏导并使其为0。可得高斯混合模型参数估计的EM算法（以高维数据为例）：
step1、参数赋初始值，开始迭代
step2、E步，计算混合项系数Zij的期望 E[Zij]：
step3、M步，计算新一轮迭代的参数模型：
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()
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)
#看到拟合的可视化结果了
