EM算法以及python例子实现

参考文章
模拟两个正态分布的均值估计。θ是我们要估计的均值
在这里插入图片描述
模拟两个正态分布的均值估计,由于我们使用的是高斯分布,即p服从高斯分布
在这里插入图片描述
在这里插入图片描述
由上面这张图:E步:固定θ,优化Q;M步:固定Q,优化θ;交替将极值推向最大

# # 模拟两个正态分布的均值估计
#
from numpy import *
import numpy as np
import random
import copy

SIGMA = 6
EPS = 0.0001


# 生成方差相同,均值不同的样本
def generate_data():
    Miu1 = 20	# 先假设分布1均值为20
    Miu2 = 40	# 分布2为40
    N = 1000	# 1k个样本,随机服从一个分布
    X = mat(zeros((N, 1)))  # [1000,1]
    for i in range(N):
        temp = random.uniform(0, 1)
        if (temp > 0.5):
            X[i] = temp * SIGMA + Miu1
        else:
            X[i] = temp * SIGMA + Miu2
    return X            # [1000,1]


# EM算法
def my_EM(X):
    k = 2
    N = len(X)
    Miu = np.random.rand(k, 1)
    Posterior = mat(zeros((N, 2))) 
    # 先求后验概率Qi(z^(i))
    for iter in range(1000):    # 最大迭代次数
     	# E-step
        for i in range(N):
            dominator = 0
            # 对应第一张图紫色框内的公式
            for j in range(k):  # k类分布 z^(1),z^(2),...,z^(k)
                dominator = dominator + np.exp(-1.0 / (2.0 * SIGMA ** 2) * (X[i] - Miu[j]) ** 2)
            for j in range(k):
                numerator = np.exp(-1.0 / (2.0 * SIGMA ** 2) * (X[i] - Miu[j]) ** 2)
                Posterior[i, j] = numerator / dominator
        oldMiu = copy.deepcopy(Miu)
        # M-step
        for j in range(k):
            numerator = 0
            dominator = 0
            for i in range(N):
                numerator = numerator + Posterior[i, j] * X[i]
                dominator = dominator + Posterior[i, j]
            Miu[j] = numerator / dominator
        print((abs(Miu - oldMiu)).sum())

        if (abs(Miu - oldMiu)).sum() < EPS:
            print('-----------')
            print('两个正态分布的均值估计是',Miu)
            print('迭代次数', iter)
            break


if __name__ == '__main__':
    X = generate_data()
    my_EM(X)

在这里插入图片描述

  • 3
    点赞
  • 21
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
GMM(Gaussian Mixture Model)是一种基于高斯分布的概率模型,常用于聚类或密度估计。EM(Expectation-Maximization)算法是一种迭代算法,通常用于GMM的参数估计。下面是使用Python实现GMM和EM算法的示例代码: ``` import numpy as np from sklearn.mixture import GaussianMixture # 生成随机数据 np.random.seed(0) X = np.concatenate([np.random.randn(100, 2) + [2, 2], np.random.randn(100, 2) + [-2, -2], np.random.randn(100, 2) + [2, -2]]) # 初始化GMM模型 gmm = GaussianMixture(n_components=3, covariance_type='full') # 训练模型 gmm.fit(X) # 打印聚类结果 print(gmm.predict(X)) # 打印GMM模型参数 print('Means:') print(gmm.means_) print('Covariances:') print(gmm.covariances_) print('Weights:') print(gmm.weights_) ``` 这段代码使用了`sklearn.mixture.GaussianMixture`类,它可以方便地进行GMM模型的训练和参数估计。其中,`n_components`参数指定了聚类个数,`covariance_type`参数指定了协方差矩阵类型。在上面的例子中,我们使用了`'full'`类型,即完整协方差矩阵。 下面是使用Python实现EM算法的示例代码: ``` import numpy as np # 初始化参数 np.random.seed(0) K = 3 N = 300 mu = np.array([[-2, 2], [2, 2], [0, -2]]) sigma = np.array([[[1, 0], [0, 1]], [[1, 0.5], [0.5, 1]], [[0.5, 0], [0, 0.5]]]) alpha = np.ones(K) / K x = np.zeros((N, 2)) for i in range(K): x[i * 100:(i + 1) * 100, :] = np.random.multivariate_normal(mu[i, :], sigma[i, :, :], 100) # EM算法迭代 for t in range(10): # E步:计算后验概率 gamma = np.zeros((N, K)) for k in range(K): gamma[:, k] = alpha[k] * np.exp(-0.5 * np.sum((x - mu[k, :]) ** 2 / sigma[k, :, :], axis=1)) / np.sqrt(np.linalg.det(sigma[k, :, :])) gamma /= np.sum(gamma, axis=1, keepdims=True) # M步:更新模型参数 for k in range(K): Nk = np.sum(gamma[:, k]) mu[k, :] = np.sum(gamma[:, k].reshape(-1, 1) * x, axis=0) / Nk sigma[k, :, :] = np.sum(gamma[:, k].reshape(-1, 1, 1) * np.matmul((x - mu[k, :]).reshape(-1, 2, 1), (x - mu[k, :]).reshape(-1, 1, 2)), axis=0) / Nk alpha[k] = Nk / N # 打印模型参数 print('Iteration', t + 1) print('Means:') print(mu) print('Covariances:') print(sigma) print('Weights:') print(alpha) ``` 这段代码使用了EM算法来估计GMM模型的参数。其中,`mu`、`sigma`和`alpha`分别表示高斯分布的均值、协方差矩阵和权重,`gamma`表示后验概率。在每一轮迭代中,首先计算后验概率,然后根据后验概率更新模型参数。迭代结束后,打印出模型参数。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值