期望极大算法（Expectation Maximization Algorithm,EM）

定义

输入:观测变量数据Y,隐变量数据Z,联合分布P(Y,Z|$\theta$),条件分布PP(Z,Y|$\theta$);
输出:模型参数$\theta$
(1)选择参数的初值${\theta }^{(0)},开始迭代;$
(2)E步:记${\theta }^{(i)}为第i次迭代参数$\theta$的估计值,在第$i+1$次迭代的E步,计算

$$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$

      $P(Z|Y,{\theta }^{(i)})$:给定观测数据$Y$和当前的参数估计${\theta }^{(i)}$下隐变量数据$Z$的条件概率分布；
(3)M步:求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$,确定第$i+1$次迭代的参数的估计值${\theta }^{(i+1)}$
$${\theta }^{(i+1)}=arg\ast \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
(4)重复第(2)步和第(3)步，直到收敛。

输入空间

T = { ( x 1 , x 2 , … , x N } T=\left\{(x_1,x_2,\dots,x_N\right\} T={(x1​,x2​,…,xN​}

import numpy as np
import random
import math
import timedef loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):'''初始化数据集这里通过服从高斯分布的随机函数来生成数据集:param mu0: 高斯0的均值:param sigma0: 高斯0的方差:param mu1: 高斯1的均值:param sigma1: 高斯1的方差:param alpha0: 高斯0的系数:param alpha1: 高斯1的系数:return: 混合了两个高斯分布的数据'''#定义数据集长度为1000length = 1000#初始化第一个高斯分布，生成数据，数据长度为length * alpha系数，以此来#满足alpha的作用data0 = np.random.normal(mu0, sigma0, int(length * alpha0))#第二个高斯分布的数据data1 = np.random.normal(mu1, sigma1, int(length * alpha1))#初始化总数据集#两个高斯分布的数据混合后会放在该数据集中返回dataSet = []#将第一个数据集的内容添加进去dataSet.extend(data0)#添加第二个数据集的数据dataSet.extend(data1)#对总的数据集进行打乱（其实不打乱也没事，只不过打乱一下直观上让人感觉已经混合了# 读者可以将下面这句话屏蔽以后看看效果是否有差别）random.shuffle(dataSet)#返回伪造好的数据集return dataSet

# mu0是均值μ
# sigmod是方差σ
#在设置上两个alpha的和必须为1，其他没有什么具体要求，符合高斯定义就可以
alpha0 = 0.3; mu0 = -2; sigmod0 = 0.5
alpha1 = 0.7; mu1 = 0.5; sigmod1 = 1#初始化数据集
dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)

np.shape(dataSetList)

print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f'%(alpha0, mu0, sigmod0, alpha1, mu1, sigmod1))

统计学习方法

模型

$$arg\ast \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

策略

$$L(\theta )=log(\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))$$

算法

高斯混合模型
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k),{\alpha }_k:系数,{\alpha }_k\ge 0,\sum _{k=1}^K{\alpha }_k=1;\varphi (y|{\theta }_k):高斯分布密度,{\theta }_k=({\mu }_k,{\sigma }_k^2)$$
$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}exp(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2})$$

def calcGauss(dataSetArr, mu, sigmod):'''根据高斯密度函数计算值:param dataSetArr: 可观测数据集:param mu: 均值:param sigmod: 方差:return: 整个可观测数据集的高斯分布密度（向量形式）'''result = (1 / (math.sqrt(2 * math.pi) * sigmod)) * \np.exp(-1 * (dataSetArr - mu) * (dataSetArr - mu) / (2 * sigmod**2))#返回结果return result

$$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
$$P(Z|Y,{\theta }^{(i)}):给定观测数据Y和当前的参数估计{\theta }^{(i)}下隐变量数据Z的条件概率分布；$$

def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):'''依据当前模型参数，计算分模型k对观数据y的响应度:param dataSetArr: 可观测数据y:param alpha0: 高斯模型0的系数:param mu0: 高斯模型0的均值:param sigmod0: 高斯模型0的方差:param alpha1: 高斯模型1的系数:param mu1: 高斯模型1的均值:param sigmod1: 高斯模型1的方差:return: 两个模型各自的响应度'''#计算y0的响应度#先计算模型0的响应度的分子gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)#模型1响应度的分子gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)#两者相加为E步中的分布sum = gamma0 + gamma1#各自相除，得到两个模型的响应度gamma0 = gamma0 / sumgamma1 = gamma1 / sum#返回两个模型响应度return gamma0, gamma1

$${\theta }^{(i+1)}=arg\ast \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$

def M_step(muo, mu1, gamma0, gamma1, dataSetArr):mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0))sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1))alpha0_new = np.sum(gamma0) / len(gamma0)alpha1_new = np.sum(gamma1) / len(gamma1)#将更新的值返回return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new

def EM_Train(dataSetList, iter = 500):'''根据EM算法进行参数估计:param dataSetList:数据集（可观测数据）:param iter: 迭代次数:return: 估计的参数'''#将可观测数据y转换为数组形式，主要是为了方便后续运算dataSetArr = np.array(dataSetList)#步骤1：对参数取初值，开始迭代alpha0 = 0.5; mu0 = 0; sigmod0 = 1alpha1 = 0.5; mu1 = 1; sigmod1 = 1#开始迭代step = 0while (step < iter):#每次进入一次迭代后迭代次数加1step += 1#步骤2：E步：依据当前模型参数，计算分模型k对观测数据y的响应度gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)#步骤3：M步mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = \M_step(mu0, mu1, gamma0, gamma1, dataSetArr)#迭代结束后将更新后的各参数返回return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1

alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)
print('Parameters predict:')
print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (alpha0, mu0, sigmod0, alpha1, mu1, sigmod1))

假设空间（Hypothesis Space）

$$\left\{arg\ast \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\right\}$$

输出

$$\theta$$