基于R的高级统计（EM算法）

作者: Shaoqian_Ma | 来源:发表于2020-04-26 11:38 被阅读0次

The EM Algorithm

EM算法全称叫做Expectation-Maximization，也就是说它分成基本的两个步骤就是The “E-Step” and the “M-step”

EM算法的基本思想如下：

一部分观测数据用Y来表示，同时有一部分缺失数据用Z来表示。这样观测数据Y和不可观测的缺失数据Z共同构成了完整的数据集X=(Y,Z) 或者，Z可以理解成“隐变量”

假设完整的数据X有一个概率密度函数g(y,z∣θ), θ为一个参数向量。因为含缺失数据，我们不能估计g
可观测的数据具有的概率密度为：把g对于z积分即可
$f(y∣θ)=∫g(y,z∣θ)dz$
这样一来，仅考虑可观测数据的似然函数就是：
$ℓ(θ∣y)=logf(y∣θ)$

EM algorithm工作流程：

E-step: 给定θ一个初值θ0，定义：
$Q(θ∣θ_0)=E[log{g(y,z∣θ)}∣y,θ_0]$
M-step: 对于函数Q(θ∣θ0)，求θ=argmax(Q），更新θ
重复步骤1和2直至结果收敛

在上述的E步骤中, 隐变量的分布可以写成：（如果把θ都去掉其实就是最简单的条件概率：g是隐变量分布和可观测分布的联合概率）
$h(z∣y,θ)=\frac {g(y,z∣θ)}{f(y∣θ)}$

总的来说，只要给定了初值θ，就可以计算出上面的h(z∣y,θ)（这时获得的就是z的后验分布），进而获得Q(θ∣θ0)关于θ的函数，再对其求θ=argmax进行迭代

详细可参考：https://zhuanlan.zhihu.com/p/40991784

通俗理解就是E-step（猜）、M-step（反思）：

你知道一些东西（观察的到的数据），你不知道一些东西（观察不到的），你很好奇，想知道点那些不了解的东西。怎么办呢，你就根据一些假设（parameter）先猜（E-step），把那些不知道的东西都猜出来，假装你全都知道了; 然后有了这些猜出来的数据，你反思一下，更新一下你的假设（parameter）, 让你观察到的数据更加可能(Maximize likelihood; M-stemp); 然后再猜，在反思，最后，你就得到了一个可以解释整个数据的假设了。

引自：https://www.zhihu.com/question/27976634

我看到的关于EM算法讲的最详细的帖子：请务必理解下面的部分，后面的代码才好理解（因为图片问题不太好直接展示出来，直接看链接里的内容吧）

引自作者：彭一洋
链接：https://www.zhihu.com/question/27976634/answer/163164402

Canonical Examples

Two-Part Normal Mixture Model

假设我们有一个混合高斯分布：由两个不同的高斯分布混合而成：
$f(y∣θ)=λφ(y∣μ_1,σ_2^1)+(1−λ)φ(y∣μ_2,σ_2^2)$
对应的对数似然函数为：
$logf(y_1,…,y_n∣θ)=log∑_{i=1}^nλφ(y_i∣μ_1,σ_1)+(1−λ)φ(y_i∣μ_2,σ_2)$
对于上述函数的优化，其实可以用牛顿法实现，不过EM算法提供了一种更好更稳定的途径。这里的缺失数据，也就是隐变量，就是混合模型中每个样本分别来自哪个高斯分布

The “missing data” in this case are the labels identifying which observation came from which population. Therefore, we assert that there are missing data z1,…,znsuch that
$z_i∼Bernoulli(λ)$

zi=1时样本yi来自第一个高斯分布，zi=0时样本yi来自第二个高斯分布

也就是说我们把这个混合高斯分布的生成分成两个步骤：首先我们随机抽样，看这个样本来自哪个分布zi，然后在给定z的前提下，再从对应的高斯分布里取样yi。于是完整数据的联合分布表示为：
$g(y,z∣θ)=φ(y∣μ_1,σ^2_1)^zφ(y∣μ_2,σ^2_2)^{1−z}λ^z(1−λ)^{1−z}$
R语言实现：

首先构造模拟的混合高斯分布数据：

mu1 <- 1
s1 <- 2
mu2 <- 4
s2 <- 1
lambda0 <- 0.4
n <- 100
set.seed(2017-09-12)
z <- rbinom(n, 1, lambda0)     ## "Missing" data  假设隐变量z服从伯努利分布
x <- rnorm(n, mu1 * z + mu2 * (1-z), s1 * z + (1-z) * s2)#x的均值和方差
hist(x)
rug(x)

mixed_norm_simulate.png

我们知道混合数据的对数似然为：
$logf(y_1,…,y_n∣λ)=log∑_{i=1}^nλφ(y_i∣μ_1,σ_1)+(1−λ)φ(y_i∣μ_2,σ_2)$
所以我们写一个用来生成混合概率密度的函数：

f <- function(x, lambda) {
        lambda * dnorm(x, mu1, s1) + (1-lambda) * dnorm(x, mu2, s2)
}

然后对上面那个函数求对数似然并加和：

loglike <- function(lambda) {
        sum(log(f(x, lambda)))
}
loglike <- Vectorize(loglike, "lambda")  ## Vectorize for plotting
par(mar = c(5,4, 1, 1))
curve(loglike, 0.01, 0.95, n = 200, ylab = "Log-likelihood", 
      xlab = expression(lambda))#可视化

logGMD.png

注意我们的模拟数据里的真实的lambda0=0.4，这里我们手动计算一下loglike的最大值点：

op <- optimize(loglike, c(0.1, 0.9), maximum = TRUE)
op$maximum
[1] 0.3097435

这里最大似然估计的值似乎有点bias，不过现在我们还不用去考虑这些

接下来我们构造EM算法里用Jensen（琴生）不等式推导的下界函数，设定初值λ0=0.8，就是我们先验地给定混合高斯分布的混合比来构造下界函数

lam0 <- 0.8
minor <- function(lambda) {
        p1 <- sum(log(f(x, lam0)))
        pi <- lam0 * dnorm(x, mu1, s1) / (lam0 * dnorm(x, mu1, s1) 
                                          + (1 - lam0) * dnorm(x, mu2, s2))#pi是隐变量z的后验分布，是可以算出来的
        p2 <- sum(pi * dnorm(x, mu1, s1, log = TRUE) 
                  + (1-pi) * dnorm(x, mu2, s2, log = TRUE)
                  + pi * log(lambda)
                  + (1-pi) * log(1-lambda))
        p3 <- sum(pi * dnorm(x, mu1, s1, log = TRUE) 
                  + (1-pi) * dnorm(x, mu2, s2, log = TRUE)
                  + pi * log(lam0)
                  + (1-pi) * log(1-lam0))#
        p1 + p2 - p3
}
minor <- Vectorize(minor, "lambda")
#然后可视化minorising函数
par(mar = c(5,4, 1, 1))
curve(loglike, 0.01, 0.95, ylab = "Log-likelihood", 
      xlab = expression(lambda))
curve(minor, 0.01, 0.95, add = TRUE, col = "red")
legend("topright", c("obs. log-likelihood", "minorizing function"), 
       col = 1:2, lty = 1, bty = "n")

minorization.jpg

上面关于minor的函数可能一眼看上去比较懵，帖子里的θ对应代码中的λ，pi就是P(Z|X;θ)，代表先验性给定θ0以后隐变量z的后验分布，p2就是P(Z|X;θ)logP(X,Z;θ)，P3就是P(Z|X;θ)logP(Z|X;θ)

$pi=P(Z|X;θ);p2=P(Z|X;θ)logP(X,Z;θ);p3=P(Z|X;θ)logP(Z|X;θ)$

equation.jpg

上式变成：
$θ^*=p2-p3$
而p1是一个定值，教材中的推导有一点点不同，多了一个p1=logf(y∣θ0)，不过这是一个定值，对于求θ*影响不大

既然我们有了下界函数，接下来对其求最值的话，就可以获得更新的λ：

par(mar = c(5,4, 2, 1))
curve(loglike, 0.01, 0.95, ylab = "Log-likelihood", 
      xlab = expression(lambda), xlim = c(-0.5, 1), #这是最开始的loglike
      ylim = c())
abline(v = lam0, lty = 2)
mtext(expression(lambda[0]), at = lam0, side = 3)
curve(minor, 0.01, 0.95, add = TRUE, col = "red", lwd = 2)#这是minor函数，还没求最值点的
op <- optimize(minor, c(0.1, 0.9), maximum = TRUE)#对红色的minor在(0.1, 0.9)范围内求最值点
abline(v = op$maximum, lty = 2)#添加第一个最值点对应的垂线，也就是λ1
lam0 <- op$maximum#然后重新把λ1赋值给初值，实现参数更新
curve(minor, 0.01, 0.95, add = TRUE, col = "blue", lwd = 2)#更新以后再画蓝色的下界函数，重复操作
abline(v = lam0, lty = 2)
mtext(expression(lambda[1]), at = lam0, side = 3)
op <- optimize(minor, c(0.1, 0.9), maximum = TRUE)
abline(v = op$maximum, lty = 2)
mtext(expression(lambda[2]), at = op$maximum, side = 3)
legend("topleft", 
       c("obs. log-likelihood", "1st minorizing function", "2nd minorizing function"), 
       col = c(1, 2, 4), lty = 1, bty = "n")

minorization2.jpg

这样我们就实现了对原函数进行一步步地优化，越来越逼近我们想要的那个点

Constrained Minimization With and Adaptive Barrier

实际情况经常是比高斯混合分布的问题要复杂的，假设我们要对f(θ)求最值，但是θ有范围，比如参数向量θ经常要受到条件限制gi(θ)≥0：
$g_i(θ)=u′_iθ−c_i$
其中ui是一个向量，长度和θ相同，ci是一个常数。i=1,…,ℓ 在这种情况下，优化函数变成：
$R(θ∣θn)=f(θ)−λ∑_{i=1}^ℓg_i(θ_n)logg_i(θ)−u′_iθ；λ>0$
可以理解成加了个拉格朗日乘子

参考教材：Advanced Statistical Computing

基于R的高级统计（EM算法）

The EM Algorithm

Canonical Examples

Two-Part Normal Mixture Model

Constrained Minimization With and Adaptive Barrier

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

生物信息学与算法