机器学习之协方差矩阵、黑塞矩阵、标准差椭圆和EM算法

一、统计学的基本概念
二、混合高斯函数
- 1、标准差椭圆的生成算法
- 2、高斯混合模型
- 高斯混合模型（GMM）
- GMM参数估计过程
  - GMM的贝叶斯理解
- EM算法估计GMM参数
三、EM算法

一、统计学的基本概念

学过概率统计的孩子都知道，统计里最基本的概念就是样本的均值，方差，或者再加个标准差。首先我们给你一个含有n个样本的集合X={X1,…,Xn}，依次给出这些概念的公式描述，这些高中学过数学的孩子都应该知道吧，一带而过。

均值：

X¯=∑ni=1XinX¯=∑i=1nXin

\bar{X}=\frac{\sum_{i=1}^n X_{i}}{n}
标准差：

s=∑ni=1(Xi−X¯)2n−1−−−−−−−−−−−−−√s=∑i=1n(Xi−X¯)2n−1

s=\sqrt{\frac{\sum_{i=1}^n (X_{i}-\bar{X})^2}{n-1}}
方差：

s2=∑ni=1(Xi−X¯)2n−1s2=∑i=1n(Xi−X¯)2n−1

s^2=\frac{\sum_{i=1}^n (X_{i}-\bar{X})^2}{n-1}
很显然，均值描述的是样本集合的中间点，它告诉我们的信息是很有限的，而标准差给我们描述的则是样本集合的各个样本点到均值的距离之平均。以这两个集合为例，[0，8，12，20]和[8，9，11，12]，两个集合的均值都是10，但显然两个集合差别是很大的，计算两者的标准差，前者是8.3，后者是1.8，显然后者较为集中，故其标准差小一些，标准差描述的就是这种“散布度”。之所以除以n-1而不是除以n，是因为这样能使我们以较小的样本集更好的逼近总体的标准差，即统计上所谓的“无偏估计”。而方差则仅仅是标准差的平方。

为什么需要协方差？
上面几个统计量看似已经描述的差不多了，但我们应该注意到，标准差和方差一般是用来描述一维数据的，但现实生活我们常常遇到含有多维数据的数据集，最简单的大家上学时免不了要统计多个学科的考试成绩。面对这样的数据集，我们当然可以按照每一维独立的计算其方差，但是通常我们还想了解更多，比如，一个男孩子的猥琐程度跟他受女孩子欢迎程度是否存在一些联系啊，嘿嘿~协方差就是这样一种用来度量两个随机变量关系的统计量，我们可以仿照方差的定义：

var(X)=∑ni=1(Xi−X¯)(Xi−X¯)n−1var(X)=∑i=1n(Xi−X¯)(Xi−X¯)n−1

var(X)=\frac{\sum_{i=1}^n (X_{i}-\bar{X})(X_{i}-\bar{X})}{n-1}
来度量各个维度偏离其均值的程度，标准差可以这么来定义：

cov(X,Y)=∑ni=1(Xi−X¯)(Yi−Y¯)n−1cov(X,Y)=∑i=1n(Xi−X¯)(Yi−Y¯)n−1

cov(X,Y)=\frac{\sum_{i=1}^n (X_{i}-\bar{X})(Y_{i}-\bar{Y})}{n-1}
协方差的结果有什么意义呢？如果结果为正值，则说明两者是正相关的(从协方差可以引出“相关系数”的定义)，也就是说一个人越猥琐就越受女孩子欢迎，嘿嘿，那必须的~结果为负值就说明负相关的，越猥琐女孩子越讨厌，可能吗？如果为0，也是就是统计上说的“相互独立”。

从协方差的定义上我们也可以看出一些显而易见的性质，如：

1.cov(X,X)=var(X)

2.cov(X,Y)=cov(Y,X)
协方差多了就是协方差矩阵
上一节提到的猥琐和受欢迎的问题是典型二维问题，而协方差也只能处理二维问题，那维数多了自然就需要计算多个协方差，比如n维的数据集就需要计算n!\(n−2)!∗2个协方差，那自然而然的我们会想到使用矩阵来组织这些数据。给出协方差矩阵的定义：

Cn×n=(ci,j,ci,j=cov(Dimi,Dimj))Cn×n=(ci,j,ci,j=cov(Dimi,Dimj))

C_{n\times n}=(c_{i,j},c_{i,j}=cov(Dim_{i},Dim_{j}))
这个定义还是很容易理解的，我们可以举一个简单的三维的例子，假设数据集有{x,y,z}三个维度，则协方差矩阵为

C.=⎡⎣⎢cov(x,x)cov(y,x)cov(z,x)cov(x,y)cov(y,y)cov(z,y)cov(x,z)cov(y,z)cov(z,z)⎤⎦⎥C.=[cov(x,x)cov(x,y)cov(x,z)cov(y,x)cov(y,y)cov(y,z)cov(z,x)cov(z,y)cov(z,z)]

C.=\left[\begin{matrix}cov(x,x)&cov(x,y)&cov(x,z) \\cov(y,x)&cov(y,y)&cov(y,z) \\cov(z,x)&cov(z,y)&cov(z,z) \end{matrix}\right]
可见，协方差矩阵是一个对称的矩阵，而且对角线是各个维度上的方差。

总结
理解协方差矩阵的关键就在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间，拿到一个样本矩阵，我们最先要明确的就是一行是一个样本还是一个维度，心中明确这个整个计算过程就会顺流而下，这么一来就不会迷茫了~

以上内容来自：http://pinkyjie.com/2010/08/31/covariance/

协方差矩阵
多维随机变量的协方差矩阵
对多维随机变量X=[X1,X2,X3,…,Xn]T，我们往往需要计算各维度两两之间的协方差，这样各协方差组成了一个n×n的矩阵，称为协方差矩阵。协方差矩阵是个对称矩阵，对角线上的元素是各维度上随机变量的方差。我们定义协方差矩阵为Σ，这个符号与求和∑相同，需要根据上下文区分。矩阵内的元素Σij为

Σij=cov(Xi,Xj)=E[(Xi−E[Xi])(Xj−E[Xj])]Σij=cov(Xi,Xj)=E[(Xi−E[Xi])(Xj−E[Xj])]

\Sigma_{ij}={cov}(X_i,X_j)={E}\big[(X_i-{E}[X_i])(X_j-{E}[X_j])\big]
这样这个矩阵为

Σ=E[(X−E[X])(X−E[X])T]Σ=E[(X−E[X])(X−E[X])T]

\Sigma={E}\big[(\textbf X-{E}[\textbf X]\big)(\textbf X-{E}[\textbf X])^T]

=cov(X1,X1)cov(X2,X1)⋮cov(Xn,X1)cov(X1,X2)cov(X2,X2)⋮cov(Xn,X2)⋯⋯⋱⋯cov(X1,Xn)cov(X2,Xn)⋮cov(Xn,Xn)=cov(X1,X1)cov(X1,X2)⋯cov(X1,Xn)cov(X2,X1)cov(X2,X2)⋯cov(X2,Xn)⋮⋮⋱⋮cov(Xn,X1)cov(Xn,X2)⋯cov(Xn,Xn)

=\begin{matrix} {cov}(X_1, X_1) (X_1, X_2) &\cdots (X_1, X_n) \\{cov}(X_2, X_1) (X_2, X_2) &\cdots (X_2, X_n) \\\vdots &\vdots &\ddots &\vdots \\{cov}(X_n, X_1) (X_n, X_2) &\cdots (X_n, X_n) \end{matrix}

=E[(X1−E[X1])(X1−E[X1])]E[(X2−E[X2])(X1−E[X1])]⋮E[(Xn−E[Xn])(X1−E[X1])]E[(X1−E[X1])(X2−E[X2])]E[(X2−E[X2])(X2−E[X2])]⋮E[(Xn−E[Xn])(X2−E[X2])]⋯⋯⋱⋯E[(X1−E[X1])(Xn−E[Xn])]E[(X2−E[X2])(Xn−E[Xn])]⋮E[(Xn−E[Xn])(Xn−E[Xn])]=E[(X1−E[X1])(X1−E[X1])]E[(X1−E[X1])(X2−E[X2])]⋯E[(X1−E[X1])(Xn−E[Xn])]E[(X2−E[X2])(X1−E[X1])]E[(X2−E[X2])(X2−E[X2])]⋯E[(X2−E[X2])(Xn−E[Xn])]⋮⋮⋱⋮E[(Xn−E[Xn])(X1−E[X1])]E[(Xn−E[Xn])(X2−E[X2])]⋯E[(Xn−E[Xn])(Xn−E[Xn])]

=\begin{matrix} {E}\big[(X_1-{E}[X_1])(X_1-{E}[X_1])\big] \big[(X_1-{E}[X_1])(X_2-{E}[X_2])\big] &\cdots \big[(X_1-{E}[X_1])(X_n-{E}[X_n])\big] \\{E}\big[(X_2-{E}[X_2])(X_1-{E}[X_1])\big] \big[(X_2-{E}[X_2])(X_2-{E}[X_2])\big] &\cdots \big[(X_2-{E}[X_2])(X_n-{E}[X_n])\big] \\\vdots &\vdots &\ddots &\vdots \\{E}\big[(X_n-{E}[X_n])(X_1-{E}[X_1])\big] \big[(X_n-{E}[X_n])(X_2-{E}[X_2])\big] &\cdots \big[(X_n-{E}[X_n])(X_n-{E}[X_n])\big] & \end{matrix}
样本的协方差矩阵

与上面的协方差矩阵相同，只是矩阵内各元素以样本的协方差替换。样本集合为{x⋅j=[x1j,x2j,…,xnj]T|1⩽j⩽m}，m为样本数量，所有样本可以表示成一个n×m的矩阵。我们以Σ^表示样本的协方差矩阵，与Σ区分。

Σ^=q11q21⋮qn1q12q21⋮qn2⋯⋯⋱⋯q1nq2n⋮qnnΣ^=q11q12⋯q1nq21q21⋯q2n⋮⋮⋱⋮qn1qn2⋯qnn

\hat \Sigma=\begin{matrix} q_{11} & q_{12} & \cdots & q_{1n} \\q_{21} & q_{21} & \cdots & q_{2n} \\\vdots & \vdots & \ddots & \vdots \\q_{n1} & q_{n2} & \cdots & q_{nn} \end{matrix}

=1m−1∑mj=1(x1j−x¯1)(x1j−x¯1)∑mj=1(x2j−x¯2)(x1j−x¯1)⋮∑mj=1(xnj−x¯n)(x1j−x¯1)∑mj=1(x1j−x¯1)(x2j−x¯2)∑mj=1(x2j−x¯2)(x2j−x¯2)⋮∑mj=1(xnj−x¯n)(x2j−x¯2)⋯⋯⋱⋯∑mj=1(x1j−x¯1)(xnj−x¯n)∑mj=1(x2j−x¯2)(xnj−x¯n)⋮∑mj=1(xnj−x¯n)(xnj−x¯n)=1m−1∑j=1m(x1j−x¯1)(x1j−x¯1)∑j=1m(x1j−x¯1)(x2j−x¯2)⋯∑j=1m(x1j−x¯1)(xnj−x¯n)∑j=1m(x2j−x¯2)(x1j−x¯1)∑j=1m(x2j−x¯2)(x2j−x¯2)⋯∑j=1m(x2j−x¯2)(xnj−x¯n)⋮⋮⋱⋮∑j=1m(xnj−x¯n)(x1j−x¯1)∑j=1m(xnj−x¯n)(x2j−x¯2)⋯∑j=1m(xnj−x¯n)(xnj−x¯n)

=\frac {1}{m-1} \begin{matrix} {\sum_{j=1}^m{(x_{1j}-\bar x_1)(x_{1j}-\bar x_1)}} ^m{(x_{1j}-\bar x_1)(x_{2j}-\bar x_2)}} & \cdots ^m{(x_{1j}-\bar x_1)(x_{nj}-\bar x_n)}} \\ {\sum_{j=1}^m{(x_{2j}-\bar x_2)(x_{1j}-\bar x_1)}} ^m{(x_{2j}-\bar x_2)(x_{2j}-\bar x_2)}} & \cdots ^m{(x_{2j}-\bar x_2)(x_{nj}-\bar x_n)}} \\ \vdots & \vdots & \ddots & \vdots \\ {\sum_{j=1}^m{(x_{nj}-\bar x_n)(x_{1j}-\bar x_1)}} ^m{(x_{nj}-\bar x_n)(x_{2j}-\bar x_2)}} & \cdots ^m{(x_{nj}-\bar x_n)(x_{nj}-\bar x_n)}} \end{matrix}

=1m−1∑j=1m(x⋅j−x¯)(x⋅j−x¯)T=1m−1∑j=1m(x⋅j−x¯)(x⋅j−x¯)T

=\frac {1}{m-1} \sum_{j=1}^m (\textbf x_{\cdot j} - \bar {\textbf x}) (\textbf x_{\cdot j} - \bar {\textbf x})^T

公式中m为样本数量，x¯为样本的均值，是一个列向量，x⋅j为第j个样本，也是一个列向量。

在写程序计算样本的协方差矩阵时，我们通常用后一种向量形式计算。一个原因是代码更紧凑清晰，另一个原因是计算机对矩阵及向量运算有大量的优化，效率高于在代码中计算每个元素。

需要注意的是，协方差矩阵是计算样本不同维度之间的协方差，而不是对不同样本计算，所以协方差矩阵的大小与维度相同。

很多时候我们只关注不同维度间的线性关系，且要求这种线性关系可以互相比较。所以，在计算协方差矩阵之前，通常会对样本进行归一化，包括两部分：

y⋅j=x⋅j−x¯。即对样本进行平移，使其重心在原点；
zi⋅=yi⋅/σi。其中σi是维度i的标准差。这样消除了数值大小的影响。
这样，协方差矩阵Σ^可以写成

Σ^=1m−1∑j=1mz⋅jzT⋅jΣ^=1m−1∑j=1mz⋅jz⋅jT

\hat \Sigma=\frac {1}{m-1} \sum_{j=1}^{m}\textbf z_{\cdot j} \textbf z_{\cdot j}^T

该矩阵内的元素具有可比性。

以上内容来自 https://www.cnblogs.com/terencezhou/p/6235974.html

多元函数的黑塞矩阵

定理
设n多元实函数在点的邻域内有二阶连续偏导，若有：

二、混合高斯函数

1、标准差椭圆的生成算法

1、确定圆心。
2、确定旋转角度。
3、确定XY轴的长度。

首先是确定圆心，方向分布工具的圆心，直接利用的是算数平均中心来计算椭圆的圆心，公式如下：

SDEx=∑ni=1(xi−X¯)2n−−−−−−−−−−−−−√SDEx=∑i=1n(xi−X¯)2n

SDE_x=\sqrt{\frac{\sum_{i=1}^n(x_i-\bar{X})^2}{n}}

SDEx=∑ni=1(yi−Y¯)2n−−−−−−−−−−−−√SDEx=∑i=1n(yi−Y¯)2n

SDE_x=\sqrt{\frac{\sum_{i=1}^n(y_i-\bar{Y})^2}{n}}
其中，Xi和Yi是每个要素的空间位置坐标，X和Y是算数平均中心

SDEx和SDEy就是最后计算出来的椭圆的圆心。
然后确定椭圆的方向，以X轴为准，正北方（12点方向）为0度，顺时针旋转，计算公式如下：

最后确定XY轴的长度，公式如下：

把所有的数据都带入到公式中，就很容易的把所有的参数都计算出来，接下去只需要再地图上画出结果就行。
那么这个椭圆揭示了一些什么意义呢？
使用ArcGIS的话，方向分布工具除了生成这样一个椭圆以外，还会给出如下结果：

其中，Shape_Leng和Shape_Area是生成的椭圆的周长和面积，单位与你数据的单位相同，这里我的数据是经纬度的，所以生成的结果只能作为相对参考结果。
CenterX和CenterY表示的是椭圆的中心点。
XstdDist和YStdDist表示的X轴的长度和Y轴的长度。
Rotation表示的是椭圆的方向角度

2、高斯混合模型

高斯混合模型（Gaussian Mixed Model）指的是多个高斯分布函数的线性组合，理论上GMM可以拟合出任意类型的分布，通常用于解决同一集合下的数据包含多个不同的分布的情况（或者是同一类分布但参数不一样，或者是不同类型的分布，比如正态分布和伯努利分布）。

如图1，图中的点在我们看来明显分成两个聚类。这两个聚类中的点分别通过两个不同的正态分布随机生成而来。但是如果没有GMM，那么只能用一个的二维高斯分布来描述图1中的数据。图1中的椭圆即为二倍标准差的正态分布椭圆。这显然不太合理，毕竟肉眼一看就觉得应该把它们分成两类。
图1

这时候就可以使用GMM了！如图2，数据在平面上的空间分布和图1一样，这时使用两个二维高斯分布来描述图2中的数据，分别记为N(μ1,Σ1)和N(μ2,Σ2). 图中的两个椭圆分别是这两个高斯分布的二倍标准差椭圆。可以看到使用两个二维高斯分布来描述图中的数据显然更合理。实际上图中的两个聚类的中的点是通过两个不同的正态分布随机生成而来。如果将两个二维高斯分布N(μ1,Σ1)和N(μ2,Σ2)合成一个二维的分布，那么就可以用合成后的分布来描述图2中的所有点。最直观的方法就是对这两个二维高斯分布做线性组合，用线性组合后的分布来描述整个集合中的数据。这就是高斯混合模型（GMM）。

图2

高斯混合模型（GMM）

设有随机变量X，则混合高斯模型可以用下式表示：

p(x)=∑k=1KπkN(x|μk,Σk)p(x)=∑k=1KπkN(x|μk,Σk)

p(x) = \sum_{k=1}^K\pi_k \mathcal{N}({x}|{\mu}_k, {\Sigma}_k)
其中N(x|μk,Σk)称为混合模型中的第k个分量（component）。如前面图2中的例子，有两个聚类，可以用两个二维高斯分布来表示，那么分量数K=2. πk是混合系数（mixture coefficient），且满足：

∑k=1Kπk=1∑k=1Kπk=1

\sum_{k=1}^K\pi_k = 1
0≤πk≤1
可以看到πk相当于每个分量N(x|μk,Σk)的权重。

GMM的应用
GMM常用于聚类。如果要从 GMM 的分布中随机地取一个点的话，实际上可以分为两步：首先随机地在这 K 个 Component 之中选一个，每个 Component 被选中的概率实际上就是它的系数πk ，选中 Component 之后，再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为已知的问题。

将GMM用于聚类时，假设数据服从混合高斯分布（Mixture Gaussian Distribution），那么只要根据数据推出 GMM 的概率分布来就可以了；然后 GMM 的 K 个 Component 实际上对应K个 cluster 。根据数据来推算概率密度通常被称作 density estimation 。特别地，当我已知（或假定）概率密度函数的形式，而要估计其中的参数的过程被称作『参数估计』。

例如图2的例子，很明显有两个聚类，可以定义K=2. 那么对应的GMM形式如下：

p(x)=π1N(x|μ1,Σ1)+π2N(x|μ2,Σ2)p(x)=π1N(x|μ1,Σ1)+π2N(x|μ2,Σ2)

p({x}) =\pi_1 \mathcal{N}({x}|{\mu}_1, {\Sigma}_1) + \pi_2 \mathcal{N}({x}|{\mu}_2, {\Sigma}_2)
上式中未知的参数有六个：(π1,μ1,Σ1;π2,μ2,Σ2). 之前提到GMM聚类时分为两步，第一步是随机地在这K个分量中选一个，每个分量被选中的概率即为混合系数πk. 可以设定π1=π2=0.5，表示每个分量被选中的概率是0.5，即从中抽出一个点，这个点属于第一类的概率和第二类的概率各占一半。但实际应用中事先指定πk的值是很笨的做法，当问题一般化后，会出现一个问题：当从图2中的集合随机选取一个点，怎么知道这个点是来自N(x|μ1,Σ1)还是N(x|μ2,Σ2)呢？换言之怎么根据数据自动确定π1和π2的值？这就是GMM参数估计的问题。要解决这个问题，可以使用EM算法。通过EM算法，我们可以迭代计算出GMM中的参数：(πk,xk,Σk).

GMM参数估计过程

GMM的贝叶斯理解

在介绍GMM参数估计之前，我们先改写GMM的形式，改写之后的GMM模型可以方便地使用EM估计参数。GMM的原始形式如下：

p(x)=∑k=1KπkN(x|μk,Σk)(1)(1)p(x)=∑k=1KπkN(x|μk,Σk)

p({x}) = \sum_{k=1}^K\pi_k \mathcal{N}({x}|{\mu}_k, {\Sigma}_k) \tag{1}
前面提到πk可以看成是第k类被选中的概率。我们引入一个新的K维随机变量z. zk(1≤k≤K)只能取0或1两个值；zk=1表示第k类被选中的概率，即：p(zk=1)=πk；如果zk=0表示第k类没有被选中的概率。更数学化一点，zk要满足以下两个条件：

zk∈{0,1}zk∈{0,1}

z_k \in \{0,1\}

∑Kzk=1∑Kzk=1

\sum_K z_k = 1
例如图2中的例子，有两类，则z的维数是2. 如果从第一类中取出一个点，则z=(1,0)；，如果从第二类中取出一个点，则z=(0,1).

zk=1的概率就是πk，假设zk之间是独立同分布的（iid），我们可以写出z的联合概率分布形式：

p(z)=p(z1)p(z2)...p(zK)=∏k=1Kπzkk(2)(2)p(z)=p(z1)p(z2)...p(zK)=∏k=1Kπkzk

p({z}) = p(z_1) p(z_2)...p(z_K) = \prod_{k=1}^K \pi_k^{z_k} \tag{2}
因为zk只能取0或1，且z中只能有一个zk为1而其它zj(j≠k)全为0，所以上式是成立的。

图2中的数据可以分为两类，显然，每一類中的数据都是服从正态分布的。这个叙述可以用条件概率来表示：

p(x|zk=1)=N(x|μk,Σk)p(x|zk=1)=N(x|μk,Σk)

p({x}|z_k = 1) = \mathcal{N}({x}|{\mu}_k, {\Sigma}_k)

即第k类中的数据服从正态分布。进而上式有可以写成如下形式：

p(x|z)=∏k=1KN(x|μk,Σk)zk(3)(3)p(x|z)=∏k=1KN(x|μk,Σk)zk

p({x}| {z}) = \prod_{k=1}^K \mathcal{N}({x}|{\mu}_k, {\Sigma}_k)^{z_k} \tag{3}
上面分别给出了p(z)和p(x|z)的形式，根据条件概率公式，可以求出p(x)的形式：

p(x)=∑zp(z)p(x|z)=∑i=1KπkN(x|μk,Σk) (zk=0的项为1，省略)(4)(4)p(x)=∑zp(z)p(x|z)=∑i=1KπkN(x|μk,Σk)(zk=0的项为1，省略)

\begin{split} p(\boldsymbol{x}) &= \sum_{\boldsymbol{z}} p(\boldsymbol{z}) p(\boldsymbol{x}| \boldsymbol{z}) \\ &= \sum_{i=1}^K \pi_k \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \mbox{ ($z_k=0$的项为1，省略)} \end{split} \tag{4}
可以看到GMM模型的(1)式与(4)式有一样的形式，且(4)式中引入了一个新的变量z，通常称为隐含变量（latent variable）。对于图2中的数据，『隐含』的意义是：我们知道数据可以分成两类，但是随机抽取一个数据点，我们不知道这个数据点属于第一类还是第二类，它的归属我们观察不到，因此引入一个隐含变量z来描述这个现象。

注意到在贝叶斯的思想下，p(z)是先验概率， p(x|z)是似然概率，很自然我们会想到求出后验概率p(z|x)：

\begin{split} \gamma(z_k) &= p(z_k = 1| \boldsymbol{x}) \\ &= \frac{p(z_k = 1) p(\boldsymbol{x}|z_k = 1)}{p(\boldsymbol{x}, z_k = 1)} \\ &= \frac{p(z_k = 1) p(\boldsymbol{x}|z_k = 1)}{\sum_{j=1}^K p(z_j = 1) p(\boldsymbol{x}|z_j = 1)} \mbox{ (全概率公式)} \\ &= \frac{\pi_k \mathcal{N}(\boldsymbol{x|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k})}{\sum_{j=1}^K \pi_j \mathcal{N}(\boldsymbol{x|\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j})} \mbox{ (结合(3)(4))} \end{split} \tag{5}
上式中我们定义符号γ(zk)来表示来表示第k个分量的后验概率。在贝叶斯的观点下，πk可视为zk=1的先验概率。

上述内容改写了GMM的形式，并引入了隐含变量z和已知x后的的后验概率γ(zk)，这样做是为了方便使用EM算法来估计GMM的参数。

EM算法估计GMM参数

EM算法（Expectation-Maximization algorithm）分两步，第一步先求出要估计参数的粗略值，第二步使用第一步的值最大化似然函数。因此要先求出GMM的似然函数。

假设x={x1,x2,…,xN}，对于图2，x是图中所有点（每个点有在二维平面上有两个坐标，是二维向量，因此x1,x2等都用粗体表示）。GMM的概率模型如(1)式所示。GMM模型中有三个参数需要估计，分别是π，μ和Σ. 将(1)式稍微改写一下：

p(x|π,μ,Σ)=∑k=1KπkN(x|μk,Σk)(6)(6)p(x|π,μ,Σ)=∑k=1KπkN(x|μk,Σk)

p({x}|{\pi}, {\mu}, {\Sigma}) = \sum_{k=1}^K\pi_k \mathcal{N}({x}|{\mu}_k, {\Sigma}_k) \tag{6}
为了估计这三个参数，需要分别求解出这三个参数的最大似然函数。先求解μk的最大似然函数。对(6)式取对数后再对μk求导并令导数为0即得到最大似然函数。
(7)

0=−∑n=1NπkN(xn|μk,Σk)∑jπjN(xn|μj,Σj)Σk(xn−μk)(7)(7)0=−∑n=1NπkN(xn|μk,Σk)∑jπjN(xn|μj,Σj)Σk(xn−μk)

0 = -\sum_{n=1}^N \frac{\pi_k \mathcal{N}({x}_n | {\mu}_k, {\Sigma}_k)}{\sum_j \pi_j \mathcal{N}({x}_n|{\mu}_j, {\Sigma}_j)} {\Sigma}_k ({x}_n - {\mu}_k) \tag{7}
注意到上式中分数的一项的形式正好是(5)式后验概率的形式。两边同乘Σ−1k，重新整理可以得到：

μk=1Nk∑n=1Nγ(znk)xn(8)(8)μk=1Nk∑n=1Nγ(znk)xn

{\mu}_k = \frac{1}{N_k}\sum_{n=1}^N \gamma(z_{nk}) {x}_n \tag{8}
其中：

Nk=∑n=1Nγ(znk)(9)(9)Nk=∑n=1Nγ(znk)

N_k = \sum_{n=1}^N \gamma(z_{nk}) \tag{9}
(8)式和(9)式中，N表示点的数量。γ(znk)表示点n（xn）属于聚类k的后验概率。则Nk可以表示属于第k个聚类的点的数量。那么μk表示所有点的加权平均，每个点的权值是∑Nn=1γ(znk)，跟第k个聚类有关。

同理求Σk的最大似然函数，可以得到：

Σk=1Nk∑n=1Nγ(znk)(xn−μk)(xn−μk)T(10)(10)Σk=1Nk∑n=1Nγ(znk)(xn−μk)(xn−μk)T

{\Sigma}_k = \frac{1}{N_k} \sum_{n=1}^N \gamma(z_{nk}) (x_n - {\mu}_k)(x_n - {\mu}_k)^T \tag{10}
最后剩下πk的最大似然函数。注意到πk有限制条件∑Kk=1πk=1，因此我们需要加入拉格朗日算子：

lnp(x|π,μ,Σ)+λ(∑k=1Kπk−1)ln⁡p(x|π,μ,Σ)+λ(∑k=1Kπk−1)

\ln p({x}|{\pi}, {\mu}, {\Sigma}) + \lambda(\sum_{k=1}^K \pi_k -1)
求上式关于πk的最大似然函数，得到：

0=∑n=1NN(xn|μk,Σk)∑jπjN(xn|μj,Σj)+λ(11)(11)0=∑n=1NN(xn|μk,Σk)∑jπjN(xn|μj,Σj)+λ

0 = \sum_{n=1}^N \frac{\mathcal{N}({x}_n | {\mu}_k, {\Sigma}_k)}{\sum_j \pi_j \mathcal{N}({x}_n|{\mu}_j, {\Sigma}_j)} + \lambda \tag{11}
上式两边同乘πk，可以得到λ=−N，进而可以得到πk更简洁的表达式：

πk=NkN(12)(12)πk=NkN

\pi_k = \frac{N_k}{N} \tag{12}
EM算法估计GMM参数即最大化(8)，(10)和(12)。需要用到(5)，(8)，(10)和(12)四个公式。我们先指定π，μ和Σ的初始值，带入(5)中计算出γ(znk)，然后再将γ(znk)带入(8)，(10)和(12)，求得πk，μk和Σk；接着用求得的πk，μk和Σk再带入(5)得到新的γ(znk)，再将更新后的γ(znk)带入(8)，(10)和(12)，如此往复，直到算法收敛。

三、EM算法

定义分量数目K，对每个分量k设置πk，μk和Σk的初始值，然后计算(6)式的对数似然函数。
E step
根据当前的πk、μk、Σk计算后验概率γ(znk)

γ(znk)=πkN(xn|μn,Σn)∑Kj=1πjN(xn|μj,Σj)γ(znk)=πkN(xn|μn,Σn)∑j=1KπjN(xn|μj,Σj)

\gamma(z_{nk}) = \frac{\pi_k\mathcal{N}({x}_n| {\mu}_n, {\Sigma}_n)}{\sum_{j=1}^K \pi_j \mathcal{N}({x}_n| {\mu}_j, {\Sigma}_j)}
M step
根据E step中计算的γ(znk)再计算新的πk、μk、Σk

μnewkΣnewkπnewk=1Nk∑n=1Nγ(znk)xn=1Nk∑n=1Nγ(znk)(xn−μnewk)(xn−μnewk)T=NkNμknew=1Nk∑n=1Nγ(znk)xnΣknew=1Nk∑n=1Nγ(znk)(xn−μknew)(xn−μknew)Tπknew=NkN

\begin{aligned} {\mu}_k^{new} &= \frac{1}{N_k} \sum_{n=1}^N \gamma(z_{nk}) {x}_n \\ {\Sigma}_k^{new} &= \frac{1}{N_k} \sum_{n=1}^N \gamma(z_{nk}) ({x}_n - {\mu}_k^{new}) ({x}_n - {\mu}_k^{new})^T \\ \pi_k^{new} &= \frac{N_k}{N} \end{aligned}

其中：

Nk=∑n=1Nγ(znk)Nk=∑n=1Nγ(znk)

N_k = \sum_{n=1}^N \gamma(z_{nk})
计算(6)式的对数似然函数

lnp(x|π,μ,Σ)=∑n=1Nln{∑k=1KπkN(xk|μk,Σk)}ln⁡p(x|π,μ,Σ)=∑n=1Nln⁡{∑k=1KπkN(xk|μk,Σk)}

\ln p({x}|{\pi}, {\mu}, {\Sigma}) = \sum_{n=1}^N \ln \left\{\sum_{k=1}^K \pi_k \mathcal{N}({x}_k| {\mu}_k, {\Sigma}_k)\right\}
检查参数是否收敛或对数似然函数是否收敛，若不收敛，则返回第2步。
以上内容来自： https://blog.csdn.net/jinping_shi/article/details/59613054