EM算法

最大似然估计

Maximum Likelihood Estimation，最大似然估计，即利用已知的样本结果，反推最有可能(最大概率)导致这样结果的参数值的计算过程。

直白来讲，就是给定了一定的数据，假定知道数据是从某种分布中随机抽取出来的，但是不知道这个分布具体的参数值，即：模型已知，参数未知，而MLE就是用来估计模型的参数。

MLE的目标是找出一组参数(模型中的参数)，使得模型产出观察数据的概率最大。

\[arg~max_θP(X;θ)
\]

MLE求解过程

• 写出似然函数(即联合概率函数)

似然函数：在样本固定的情况下，样本出现的概率与参数θ之间的函数

• 对似然函数取对数，并整理；
• 求导；
• 解似然方程。

\[L(X;θ)→l(θ)=ln (L(X;θ))→\frac{∂l}{∂θ}
\]

举个例子：

假定盒子中有黑白两种球，数目未知，黑白球比例也未知，现只知道随机的10次有放回抽样情况，求盒子中抽取出白球的概率：

盒子编号

1

2

3

4

5

6

7

8

9

10

1

黑

黑

黑

黑

黑

黑

黑

黑

黑

黑

2

黑

黑

白

白

黑

黑

黑

白

黑

黑

3

黑

白

黑

黑

白

黑

白

黑

白

白

4

白

黑

白

白

黑

白

黑

白

白

白

5

白

白

白

白

白

白

白

白

白

白

现假定白球的比例为p，则黑球的比例为1-p，有如下：

\[L(X;p)=P(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10};p)=\prod_{n=1}^{10}{P(x_i;p)}
\]

\[l(P)=ln(L(X;p))=\sum_{i=1}^{10}{P(x_i;p)}
\]

\[令：\frac{∂l}{∂p}=0
\]

\[盒子1中抽取白球概率：L(X;p)=(1-p)^{10}→l(P)=ln(L(X;p))=10ln(1-p)→由于0<p<1,因此当p=0时，l取最大值
\]

\[盒子2中抽取白球概率：L(X;p)=p^3(1-p)^{7}→l(P)=3lnp+7ln(1-p)→\frac{∂l}{∂p}=\frac{3}{p}-\frac{7}{1-p}=0→p=0.3
\]

\[盒子3中抽取白球概率：L(X;p)=p^5(1-p)^{5}→l(P)=5lnp+5ln(1-p)→\frac{∂l}{∂p}=\frac{5}{p}-\frac{5}{1-p}=0→p=0.5
\]

\[盒子4中抽取白球概率：L(X;p)=p^7(1-p)^{3}→l(P)=7lnp+3ln(1-p)→\frac{∂l}{∂p}=\frac{7}{p}-\frac{3}{1-p}=0→p=0.7
\]

\[盒子5中抽取白球概率：L(X;p)=p^{10}→l(P)=ln(L(X;p))=10lnp→由于0<p<1,因此当p=1时，l取最大值
\]

盒子编号

白球

黑球

1

0

1

2

0.3

0.7

3

0.5

0.5

4

0.7

0.3

5

1

0

贝叶斯算法估计

贝叶斯算法估计是一种从先验概率和样本分布情况来计算后验概率的一种方式。

常见概念：

P(A)事件A的先验概率或者边缘概率；

P(A|B)已知B发生后A发生的条件概率，也称A的后验概率；

P(B|A)已知A发生后B发生的条件概率，也称B的后验概率；

P(B)事件B的先验概率或者边缘概率；

\[P(AB)=P(A)P(B|A)=P(B)P(A|B)⇒P(A|B)=\frac{P(B|A)P(A)}{P(B)}
\]

\[离散条件下的贝叶斯公式：P(A_i|B)=\frac{P(B|A_i)P(A_i)}{\sum_{j}{P(B|A_i)P(A_i)}}
\]

再举个例子：

现在有5个盒子，假定每个盒子中都有黑白两种球，其比例如下；现已知从这5个盒子中的任意一个盒子中有放回的抽取两个球，且均为白球，问这两个球是从哪个盒子中抽取出来的？

盒子编号

白球

黑球

1

0

1

2

0.3

0.7

3

0.5

0.5

4

0.7

0.3

5

1

0

使用MLE估计的话，如下：

\[L(X;p)=p^2⇒p=1
\]

使用贝叶斯算法估计(执果索因)，假定抽出白球的事件为B，从第i个盒子中抽取为事件\(A_i\)，有如下：

\[P(A_1|B)=\frac{P(A_1)P(B|A_1)}{P(B)}=\frac{0.2*0*0}{0.2*0^2+0.2*0.3^2+0.2*0.5^2+0.2*0.7^2+0.2*1^2}=\frac{0}{0.366}=0
\]

\[P(A_2|B)=\frac{P(A_2)P(B|A_2)}{P(B)}=0.049~~~P(A_3|B)=\frac{P(A_3)P(B|A_3)}{P(B)}=0.137
\]

\[P(A_4|B)=\frac{P(A_4)P(B|A_4)}{P(B)}=0.268~~~P(A_5|B)=\frac{P(A_5)P(B|A_5)}{P(B)}=0.564
\]

最大后验概率估计

Maximum a posteriori estimation，最大后验概率(MAP)和最大先验概率(MLE)一样，都是通过样本估计参数θ的值；

在MLE中，是使似然函数P(x|θ)最大的时候参数θ的值，MLE中假设先验概率是一个等值的；

而MAP中，则是求θ使P(x|θ)P(θ)的值最大，这也就是要求θ不仅仅是让似然函数函数最大，同时要求θ本身出现的先验概率也得比较大。

可以认为MAP是贝叶斯算法的一种应用，只是去除了贝叶斯算法中分母的部分，如下：

\[P(θ’|X)=\frac{P(θ')P(X|θ')}{P(X)}⇒argmax_{θ'}P(θ'|X)⇒argmax_{θ'}P(θ')P(X|θ')
\]

引入一个例子

背景：公司中有男同事=[A, B, C]，同时有很多漂亮的女职员=[小甲，小章，小乙]。你怀疑这些男同事和这些女职员有"问题"。为了科学的验证你的猜想，你进行了细致的观察。

有观察数据如下：

1. A，小甲，小乙一起出门了；
2. B，小甲，小章一起出门了；
3. B，小章，小乙一起出门了；
4. C，小乙一起出门了；

收集到数据后，进行了EM计算：

初始化：所有人的条件都一样，每个人和每个人都有关系。因此，每个男同事和每个女职员"有问题"的概率都是1/3；

E-step：

1. A跟小甲出去过1/2*1/3=1/6次，和小乙也出去了1/6次；
2. B跟小甲，小章也出去了1/6次；
3. B跟小章，小乙又出去了1/6次；
4. C跟小乙出去了1/3次；

M-step：更新你的八卦

\[A跟小甲，小乙有问题的概率都是\frac{\frac{1}{6}}{\frac{1}{6}+\frac{1}{6}}=\frac{1}{2}
\]

\[B跟小甲，小乙有问题的概率：\frac{\frac{1}{6}}{\frac{1}{6}*4}=\frac{1}{4}；跟小章有问题的概率是：\frac{\frac{1}{6}*2}{\frac{1}{6}*4}=\frac{1}{2}
\]

\[C跟小乙又问题的概率是：1
\]

E-step：根据最新的概率进行计算

1. A跟小甲出去过1/2*1/2=1/4次，和小乙也出去了1/4次；
2. B跟小甲出去了1/2*1/4=1/8次，跟小章出去了1/2 * 1/2 = 1/4次；
3. B跟小乙出去了1/2*1/4=1/8次，跟小章又出去了1/2 * 1/2 = 1/4次；
4. C跟小乙出去了1次；

M-step：重新更新你的八卦

\[A跟小甲，小乙有问题的概率都是\frac{\frac{1}{4}}{\frac{1}{4}+\frac{1}{4}}=\frac{1}{2}
\]

\[B跟小甲，小乙有问题的概率：\frac{\frac{1}{8}}{\frac{1}{8}*2+\frac{1}{4}*2}=\frac{1}{6}；跟小章有问题的概率是：\frac{2}{3}
\]

\[C跟小乙又问题的概率是：1
\]

好了，至此，你似乎已经得到了真相。

EM算法，Expectation Maximization Algorithm，最大期望算法，这是一种迭代类型的算法，是一种在概率模型中寻找参数最大似然估计或者最大后验估计的算法，但是在模型中存在着无法观测的隐藏变量。

EM算法流程：

• 初始化分布参数；
• 重复下列两个操作直至收敛：
  • E步骤：估计隐藏变量的概率分布期望函数；
  • M步骤：根据期望函数重新估计分布参数。

EM算法原理

给定m个训练样本，如下：

\[\{ x^{(1)},x^{(2)},…,x^{(m)} \}
\]

样本间独立，找出样本的模型参数θ，极大化模型分布的对数似然函数，如下：

\[θ=arg max_θ\sum_{i=1}^{m}{log(P(x^{(i)};θ))}
\]

假定样本数据中存在隐含数据Z，如下：

\[Z=\{ z^{(1)},z^{(2)},…,z^{(k)} \}
\]

此时极大化模型分布的对数似然函数修改如下：

\[θ=arg max_θ\sum_{i=1}^{m}{log(P(x^{(i)};θ))}
\]

\[=arg max_θ\sum_{i=1}^{m}{log[\sum_{z^{(j)}}{P(z^{(i)})P(x^{(i)}|z^{(j)};θ)]}}
\]

\[=argmax_θ\sum_{i=1}^{m}{log[\sum_{z^{(j)}}{P(x^{(i)},z^{(j)};θ)}]}
\]

令Z的分布为Q(z;θ)，那么有如下公式：

\[其中Q(z^{(j)};θ)满足如下条件：Q(z^{(j)};θ)≥0，\sum_{z^{(j)}}{Q(z^{(j)};θ)}=1
\]

\[l(θ)=\sum_{i=1}^{m}{log[\sum_{z^{(j)}}{P(x^{(i)},z^{(j)};θ)}]} \geq \sum_{i=1}^{m}{\sum_{z^{j}}{Q(z^{(j)};θ)log(\frac{P(x^{(i)},z^{(j)};θ)}{Q(z^{(j)};θ)})}}
\]

推导过程如下：

\[l(θ)=\sum_{i=1}^{m}{log[\sum_{z^{(j)}}{Q(z^{(j)};θ) \cdot \frac{P(x^{(i)},z^{(j)};θ)}{Q(z^{(j)};θ)}}]}
\]

\[由上Jensen不等式⇒原始=\sum_{i=1}^{m}{log[E_Q(\frac{P(x^{(i)},z;θ)}{Q(z;θ)})]}
\]

\[由于f函数为log函数⇒原始\geq \sum_{i=1}^{m}{E_Q[log(\frac{P(x^{(i)},z;θ)}{Q(z;θ)})}]
\]

\[将上式还原⇒上式=\sum_{i=1}^{m}{\sum_{z^{j}}{Q(z^{(j)};θ)log(\frac{P(x^{(i)},z^{(j)};θ)}{Q(z^{(j)};θ)})}}
\]

Jensen不等式介绍：

若函数f下凸，如下图：

则存在如下公式：

\[f(θx+(1-θ)y) \leq θf(x)+(1-θ)f(y)
\]

泛化后：

\[f(θ_1x_1+…+θ_kx_k) \leq θ_1f(x_1)+…+θ_kf(θ_k)
\]

\[f(E(x)) \leq E(f(x))
\]

\[其中：θ_1,…,θ_k \geq 0,θ_1+…+θ_k=1
\]

在上述式子成立的前提下，由Jensen不等式可推出：当下列式子的值为常数时，上式才能取等号，即l(θ)的下界，如下，

\[\frac{P(x,z;θ)}{Q(z;θ)}=c
\]

\[⇒Q(z;θ)=\frac{P(x,z;θ)}{c}
\]

\[由于\sum_{z^{(j)}}{Q(z^{(j)};θ)}=1⇒Q(z;θ)=\frac{P(x,z;θ)}{c\cdot\sum_{z^{(j)}}{Q(z^{(j)};θ)}}
\]

\[=\frac{P(x,z;θ)}{\sum_{z^{(j)}}{cP(x,z^{(j)};θ)}}=\frac{P(x,z;θ)}{P(x;θ)}=\frac{P(z|x;θ)P(x;θ)}{P(x;θ)}=P(z|x;θ)
\]

即，Z的分布其实是在θ确定的情况下，x给定下z的概率分布

而我们最初的目的是：找出样本的模型参数θ，极大化模型分布的对数似然函数，经过上述推导，有如下

\[θ^{new}=argmax_θl(θ)=argmax_θ\sum_{i=1}^{m}{\sum_{z^{j}}{Q(z^{(j)};θ^{old})log(\frac{P(x^{(i)},z^{(j)};θ)}{Q(z^{(j)};θ^{old})})}}
\]

\[=argmax_θ\sum_{i=1}^{m}{\sum_{z^{(j)}}{P(z^{(j)}|x^{(i)};θ^{old})log(\frac{P(x^{(i)},z^{(j)};θ)}{P(z^{(j)}|x;θ^{old})})}}
\]

\[=argmax_θ\sum_{i=1}^{m}{\sum_{z^{j}}{P(z^{(j)}|x^{(i)};θ^{old})log(P(x^{(i)},z^{(j)};θ))}}-C
\]

EM算法流程

\[有样本数据x=\{ x_1,x_2,…,x_m \},联合分布P(x,z;θ),条件分布P(z|x;θ)，最大迭代次数J
\]

• 随机初始化模型参数θ的初始值\(θ^0\)

• 开始EM算法的迭代处理：

  • E步：计算联合分布的条件概率期望

  \[Q^j=P(Z|X;θ^j)~~~l(θ)=\sum_{i=1}^{m}{\sum_{z^{j}}{P(z^{(j)}|x^{(i)};θ^{old})log(P(x^{(i)},z^{(j)};θ))}}
  \]

  • M步：极大化l函数，得到新的θ值

  \[θ^{j+1}=arg max_θl(θ)
  \]

  • 如果得出的新的θ已经收敛，则算法结束，输出最终的模型参数θ，否则继续迭代处理

EM算法直观案例：

假设现有两个装有不定数量的黑白球盒子，随机从盒子中抽取出一个白球的概率为\(p_1,p_2\)；为了估计这两个概率，每次选择一个盒子，有放回的连续随机抽取5个球，记录如下：

盒子编号

1

2

3

4

5

统计

1

白

白

黑

白

黑

3白-2黑

2

黑

黑

白

白

黑

2白-3黑

1

白

黑

黑

黑

黑

1白-4黑

2

白

黑

白

黑

白

3白-2黑

1

黑

白

黑

白

黑

2白-3黑

使用MLE最大似然估计：

\[l(p_1)=log(p_1^6(1-p_1)^9)=6logp_1+9log(1-p_1)
\]

\[\frac{∂l(p_1)}{∂p_1}=0→p_1=0.4
\]

\[同理，可得p_2=0.5
\]

若此时，不知道具体盒子的编号，但是同样为了求解\(p_1,p_2\)的值，此时就相当于多了一个隐藏变量z，z表示的是每次抽取的时候选择的盒子编号，比如z1就表示第一次抽取的时候选择的盒子1还是盒子2，如下

盒子编号

1

2

3

4

5

统计

z1

白

白

黑

白

黑

3白-2黑

z2

黑

黑

白

白

黑

2白-3黑

z3

白

黑

黑

黑

黑

1白-4黑

z4

白

黑

白

黑

白

3白-2黑

z5

黑

白

黑

白

黑

2白-3黑

• 随机初始一个概率值：盒子1中取出白球概率：p1=0.1，盒子2中取出白球概率：p2=0.9；然后使用MLE计算每轮操作中两个盒子中抽取的最大概率，计算出z值，之后重新使用最大似然估计法估计概率值

\[L(z_1=1|x;p_1)=p_1^3 \times p_2^2=0.1^3 \times 0.9^2=0.00081
\]

\[L(z_1=2|x;p_2))=p_1^3 \times p_2^2=0.9^3 \times 0.1^2=0.00729
\]

轮数

盒子1概率

盒子2概率

归一化：盒1

归一化：盒2

1

0.00081

0.00729

0.1

0.9

2

0.00729

0.00081

0.9

0.1

3

0.06561

0.00009

0.999

0.001

4

0.00081

0.00729

0.1

0.9

5

0.00729

0.00081

0.9

0.1

• 重新计算p的概率值

\[l(p_1)=log[p_1^{0.1×3+0.9×2+0.999×1+0.1×3+0.9×2}(1-p_1)^{0.1×2+0.9×3+0.999×4+0.1×2+0.9×3}]
\]

\[log[p_1^{5.199}(1-p_1)^{9.796}]=5.199logp_1+9.796log(1-p_1)
\]

\[\frac{∂l(p_1)}{∂p_1}=0→p_1=0.347
\]

\[同理，计算得p_2=0.58
\]

• 根据p的概率值，再次计算在p的条件下，从每个盒子中抽出的概率，如下：

轮数

盒子1概率

盒子2概率

归一化：盒1

归一化：盒2

1

0.0178

0.0344

0.34

0.66

2

0.0335

0.0249

0.57

0.43

3

0.0630

0.0180

0.78

0.22

4

0.0178

0.0344

0.34

0.66

5

0.0335

0.0249

0.57

0.43

• 根据新的z值，采用MLE进行计算新的p值，如下

\[p_1=0.392~~~~~p_2=0.492
\]

继续迭代，一直迭代到收敛，此时的p值即为所求。

EM算法收敛证明

EM算法的本质为：寻找参数最大似然估计。因此在每次迭代的过程中，只需要迭代后的参数\(θ^{j+1}\)计算的似然函数大于迭代前参数\(θ^{j}\)计算的似然函数即可，如下：

\[\sum_{i=1}^{m}log(p(x^{(i)};θ^{j+1})) \geq \sum_{i=1}^{m}log(p(x^{(i)};θ^{j}))
\]

具体的证明流程，略

引入例子：

随机选择1000名用户，测量用户的身高；若样本中存在男性和女性，身高分别服从高斯分布\(N(μ_1,σ_1)和N(μ_2,σ_2)\)的分布，试估计参数：\(μ_1,σ_1,μ_2,σ_2\)；

• 若明确知道样本的情况(即男女性数据是分开的)，那么我们使用极大似然估计来估计这个参数值；
• 如果样本是混合而成的，不能明确的区分开，那么就没法直接使用极大似然估计来进行参数估计，此时就引出了GMM

GMM(Gaussian Mixture Model,高斯混合模型)是指该算法由多个高斯模型线性叠加混合而成。每个高斯模型称之为component(成分)。GMM算法描述的是数据本身存在的一种分布。

GMM算法常用于聚类应用中，component的个数就可以认为是类别的数量。

假定GMM有k个Gaussian分布线性叠加而成，那么概率密度函数如下：

\[p(x)=\sum_{k=1}^{K}{p(k)p(x|k)}=\sum_{k=1}^{K}{π_kp(x;μ_k,Σ_k)}~~~;π_k:选择第k个类别的概率,μ_k,Σ_k:均值和方差矩阵
\]

对数似然函数如下：

\[l(π,μ,Σ)=\sum_{i=1}^{N}{log[\sum_{k=1}^{K}{π_kp(x^i;μ_k,Σ_k)}]}
\]

GMM求解过程

E-step，在给定x的条件下，数据属于第j个类别的概率：

\[w_j^{(i)}=Q_i(z^{(i)}=j)=p(z^{(i)}=j|x^{(i)};π,μ,Σ)
\]

M-step：极大化对数似然函数l(π,μ,Σ)，更新参数π,μ,Σ：

具体的推导步骤过于繁琐，故省略

\[μ_j=\frac{\sum_{i=1}^{m}{w_j^{(i)}x^{(i)}}}{\sum_{i=1}^{m}{w_j^{(i)}}}
\]

\[Σ_j=\frac{\sum_{i=1}^{m}{w_j^{(i)}(x^{(i)}-μ_l)[(x^{(i)}-μ_j)]^T}}{\sum_{i=1}^{m}{w_j^{(i)}}}
\]

\[π_j=\frac{1}{m}\sum_{i=1}^{m}{w_j^{(i)}}
\]

在π,μ,Σ更新完成后，又可进行E-step，不停的进行迭代，直至收敛；

在收敛后，用收敛后的参数，对于待测样本x都可以得出在不同的类别j下的概率值，选一个最大的概率值即为预测值。

• EM1_案例一：EM分类初识及GMM算法实现
• EM2_案例二：GMM算法分类及参数选择案例
• EM3_案例三：GMM的不同参数
• EM4_案例四：EM无监督算法分类鸢尾花数据

GitHub