import numpy as np
class Gauss_EM():
"""
算法9.2
"""
def __init__(self, yVector, parameters, maxIter,epsilon):
"""
:param yVector: 观测数据
:param parameters: 初始化的高斯混合模型参数,二维数组,每一组参数是一个数组[alpha,mu,sigma],共有K组
:param maxIter: 最大迭代次数
:param epsilon: 参数序列的精度阈值下界,即更新参数时所有参数的改变量都小于epsilon时停止
"""
self.yVector = np.array(yVector)
self.paramaters = np.array(parameters)
self.K = len(self.paramaters)
self.N = len(self.yVector)
self.train(maxIter,epsilon)
def Gauss(self, mu, sigma):
"""
:param mu: 高斯模型参数
:param sigma: 高斯模型参数
:return: 某个高斯分模型,公式中的y在这里换成向量,长度是N
"""
import math
result = (1 / (math.sqrt(2 * math.pi) * sigma)) * \
np.exp(-1 * (self.yVector - mu) * (self.yVector - mu) / (2 * sigma ** 2))
# 返回结果
return result
# return (1 / np.sqrt(2 * np.pi) * sigma) * np.exp(-(self.yVector - mu) ** 2 / (2 * sigma ** 2))
def E_step(self):
"""
:return:响应因子矩阵,N行K列
"""
gammaMatrix = []
for i in range(self.K):
gammaVector = self.paramaters[i, 0] * self.Gauss(self.paramaters[i, 1], self.paramaters[i, 2])
gammaMatrix.append(gammaVector)
gammaMatrix = np.asarray(gammaMatrix).T
# print(gammaMatrix/gammaMatrix.sum(axis=1).reshape(-1,1))
return gammaMatrix / gammaMatrix.sum(axis=1).reshape(-1, 1)
def M_step(self, gammaMatrix):
"""
:param gammaMatrix: E步得到的相应值
:return: 新计算得到的参数,格式与parameters一样
"""
nGammaSum = gammaMatrix.sum(axis=0)
alpha_new = (nGammaSum / self.N)
mu_new = (self.yVector @ gammaMatrix) / nGammaSum
# 暂时想不到能利用矩阵运算一步算出σ的
sigma_new = []
for k in range(self.K):
sigma_new.append((self.yVector - self.paramaters[k, 1]) ** 2 @ gammaMatrix[:, k])
sigma_new = np.sqrt(np.asarray(sigma_new) / nGammaSum) # 注意要开平方
# print(mu_new.shape,sigma_new.shape,alpha_new.shape)
paramerers_new = np.asarray(list(zip(alpha_new, mu_new, sigma_new)))
return paramerers_new
def train(self, maxIter,epsilon):
iter = 0
while (iter < maxIter):
iter += 1
gammaMatrix = self.E_step()
# print("gammaMatrix={}".format(gammaMatrix))
paramerers_new = self.M_step(gammaMatrix)
changeMatrix=paramerers_new - self.paramaters
print("迭代次数:{},更新改变:{}".format(iter,changeMatrix))
if (changeMatrix<epsilon).all():
break
self.paramaters = paramerers_new
paramerers = [[0.31, -3.122, 0.1234], [0.69, 0.124, 2.23]]
N = 1000
y0 = np.random.normal(paramerers[0][1], paramerers[0][2], int(N * paramerers[0][0]))
y1 = np.random.normal(paramerers[1][1], paramerers[1][2], int(N * paramerers[1][0]))
yVector = []
yVector.extend(y0)
yVector.extend(y1)
# print(yVector)
EM = Gauss_EM(yVector, [[0.5, 0, 1], [0.5, 1, 1]], 40,1e-4)
# print(EM.Gauss(0.5,0.5))
print("实际参数:\n{}".format(paramerers))
print("拟合参数:\n{}".format(EM.paramaters))
06-03
1605
![](https://csdnimg.cn/release/blogv2/dist/pc/img/readCountWhite.png)