⽤python求期望_Python——EM（期望极⼤算法）教学（附详
细代码与注解）
今天，我们详细的讲⼀下EM算法。
前提准备
Jupyter notebook 或 Pycharm
⽕狐浏览器或⾕歌浏览器
win7或win10电脑⼀台
⽹盘提取csv数据
需求分析
实现⾼斯混合模型的 EM 算法(GMM_EM)
⾼斯混合模型是多个⾼斯模型的线性叠加⽽成的，⾼斯混合模型的概率分布表⽰如下：
其中，k表⽰模型的个数，αkα_kαk是第 k 个模型的系数，表⽰出现该模型的概率，ϕ(x;μk,Σk) 是第 k 个⾼斯模型的概率分布。
E步：样本 xix_ixi来⾃于第 k 个模型的概率，我们把这个概率称为模型 k 对样本 xix_ixi的“责任”，也叫“响应度”，记作
γ(ik)γ_(ik)γ(i k)，计算公式如下：
M步：根据样本和当前 γ 矩阵重新估计参数，注意这⾥ x 为列向量，计算公式如下：
【⽬标】给定⼀堆没有标签的样本和模型个数 K，以此求得混合模型的参数，然后就可以⽤这个模型来对样本进⾏聚类。
python代码如下：
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal #本问题考虑的是⾼斯混合模型，所以导⼊多元⾼斯分布multivariate_normal
def prob_Y_k(Y,mu_k,cov_k): #Y为样本矩阵
norm = multivariate_normal(mean = mu_k , cov = cov_k) #⽣成多元正太分布，mu为第k个模型的均值，cov为第k个模型的协⽅差矩阵(协⽅差矩阵必须是实对称矩阵)
return norm.pdf(Y) #返回样本Y来⾃于第k个模型的概率
def Estep(Y,mu,cov,alpha): #Y为样本矩阵，alpha为权重
N = Y.shape[0] #样本数
K = alpha.shape[0] #模型数
assert N>1 , "There must be more than one sample!" #为避免单个样本导致返回的结果的类型不⼀致，因此要求样本数必须⼤于⼀
assert K>1 , "There must be more than one gaussian model!" #为避免单个模型结果的类型不⼀致，因此要求模型须⼤于⼀
gamma = np.s((N,K))) #初始化响应度矩阵，⾏对应样本数，列对应模型数
prob = np.zeros((N,K)) #初始化所有样本出现的概率矩阵，⾏对应样本数，列对应响应度
在线考试程序源码
for k in range(K):
prob[:,k] = prob_Y_k(Y,mu[k],cov[k]) #第k个模型的概率prob_Y_k
prob = np.mat(prob) #K个prob放⼊数组中
for k in range(K):
gamma[:,k] = alpha[k] * prob[:,k] #计算模型k对样本i的响应度
for i in range(N):
gamma[i,:] /= np.sum(gamma[i,:]) #第i个样本的占总样本的响应程度
return gamma #gamma为响应度矩阵
def Mstep(Y,gamma): #传⼊样本矩阵Y和Estep得到的gamma响应度矩阵
N, D = Y.shape #N为样本数，D为特征数
K = gamma.shape[1] #模型数
mu = np.zeros((K,D)) #初始化参数均值mu，每个模型的D维各有均值故mu的矩阵为K⾏D列
cov = [] #初始化参数协⽅差矩阵
alpha = np.zeros(K) # 初始化权重数组，每个模型都有权值
#接下来是更新每个模型的参数
sass产品for k in range(K):
Nk = np.sum(gamma[:,k]) #第k个模型所有样本的响应度之和
mu[k,:] = np.sum(np.multiply(Y, gamma[:,k]),axis=0)/Nk #更新参数均值mu，对每个特征求均值
cov_k = (Y - mu[k]).T * np.multiply((Y - mu[k]), gamma[:,k]) / Nk #更新cov
cov = np.append(cov_k)
alpha[k] = Nk / N
cov = np.array(cov)
return mu, cov, alpha
def normalize_data(Y): #将所有数据进⾏归⼀化处理，
for i in range(Y.shape[1]):
max_data = Y[:,i].max()网页代码中的单标签
min_data = Y[:,i].min()
Y[:,i] = (Y[:,i] - min_data)/(max_data - min_data) #此处⽤到min-max归⼀化
debug("Data Normalized")
return Y
def init_params(shape,K): #在执⾏该算法之前，需要先给出⼀个初始化的模型参数。我们让每个模型的μ为随机值,Σ 为单位矩阵，α为 1/K，即每个模型初始时都是等概率出现的。
mu = np.random.rand(K, D) #⽣成⼀个K⾏D列的[0,1)之间的数组
cov = np.array([np.eye(D)] * K) #⽣成K个D维的对⾓矩阵
alpha = np.array([1.0 / K] * K) #⽣成K个权重
debug("Parameters initialized.")
debug("mu:",mu, "cov:",cov ,"alpha:",alpha,sep = "\n" )
return mu, cov, alpha
def GMM_EM(Y, K, times): #⾼斯混合EM算法,Y为给定样本矩阵，K为模型个数，times为迭代次数，⽬的是求该模型的参数Y = normalize_data(Y) #调⽤前⾯定义的normalize_data函数，归⼀化样本矩阵Y
mu, cov, alpha = init_params(Y.shape, K) #调⽤init_params函数得到初始化的参数mu,cov,alpha
for i in range(times):
gamma = Estep(Y, mu, cov, alpha) #调⽤Estep得到响应度矩阵
mu, cov, alpha = Mstep(Y, gamma) #调⽤Mstep得到更新后的参数mu,cov,alpha
debug("{sep} Result {sep}".format(sep="-"*20))
debug("mu:", mu , "cov:",cov , "alpha:",alpha , sep="\n")
return mu,cov,alpha
import matplotlib.pyplot as plt
from gmm import *
DEBUG = True
舍曲林让人感觉快乐高兴Y = np.loadtxt("gmm.data") #载⼊数据
matY = np.matrix(Y ,copy = True)
python入门教程网盘
K = 2 #模型个数(相当于聚类的类别个数)
mu, cov, alpha = GMM_EM(matY , K , 100) #调⽤GMM_EM函数，计算GMM模型参数
N = Y.shape[0]
gamma = Estep(matY, mu, cov, alpha) #求当前模型参数下，各模型对样本的响应矩阵
category = gamma.argmax(axis = 1).flatten().tolist()[0] #对每个样本，求响应度最⼤的模型下标，作为其类别标识
class1 = np.array([Y[i] for i in range(N) if category[i] == 0]) #将每个样本放⼊对应样本的列表中
class2 = np.array([Y[i] for i in range(N) if category[i] == 1])
plt.plot(class1[:,0],class1[:,1], 'rs' ,label = "class1")
plt.plot(class2[:,0],class2[:,1], 'bo' ,label = "class2")
plt.legend(loc = "best")
plt.title("GMM Clustering By EM Algorithm")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
民间小调二人转送情郎cov1 = np.mat("0.3 0 ; 0 0.1") #2维协⽅差矩阵(必须是对⾓矩阵)
cov2 = np.mat("0.2 0 ; 0 0.3")
mu1 = np.array([0,1])
mu2 = np.array([2,1])
sample = np.zeros((100,2)) #初始化100个样本，样本特征为2
sample[:30, :] = np.random.multivariate_normal(mean=mu1, cov=cov1, size=30) #⽣成多元正态分布矩阵sample[30:, :] = np.random.multivariate_normal(mean=mu2, cov=cov2, size=70)
np.savetxt("sample.data",sample) # 将array保存到txt⽂件中
plt.plot(sample[:30, 0], sample[:30, 1], "bo") #30个样本⽤蓝⾊圆圈标记
plt.plot(sample[30:, 0], sample[30:, 1], "rs") #70个样本⽤红⾊⽅块标记
plt.title("sample_data")
plt.show()

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。