0%

统计学习方法-EM算法笔记

EM算法是一种迭代算法,由E步(求期望)和M步(求极值)两步组成,因此EM算法全称为期望极大算法(Expectation Maximum algorithm)。由于其是基础算法,因此在其他机器学习方法中常会用到,比如HMM EM算法的应用场合:如果我们要估计概率模型的参数,已知观测数据,那么一般会用极大似然估计法或者贝叶斯估计法,这在之前的逻辑斯蒂回归和朴素贝叶斯算法中提到。但当除了观测数据外还有隐藏数据,那么这时则要用EM算法

在生物信息分析算法中,用到EM算法的方法也不少,比如在转录组定量方法中,RSEM和Salmon都用到了EM算法,其需要将一些多重比对的序列通过EM算法评估其在不同比对位置上的可能性

EM算法主要用了启发式的迭代方法,比如我们并不知道隐变量是什么,那么就先通过初始模型参数和观测值估计隐变量(E步);在得到预估的隐变量后,将其和观测数据一起通过极大化似然函数来估计模型参数(M步),依次迭代上述过程直至模型参数变化小于指定阈值;以上可看出,我们先固定了模型参数来估计隐变量,然后再固定隐变量来估计模型参数,依次循环。。。

这思路在K-means算法中有直观体现:我们不知道某个数据点来自于哪个质心(隐变量),因此我们先初始化质心,然后对价值函数J极小化,计算数据点所属的类(模型参数);然后根据数据点的类重新计算质心,迭代至质心变化很小,小于指定阈值等

对于EM算法的理解,网上有张来自于EM文章(What is the expectation maximization algorithm?)的图片解释有助于理解,如下:

EM-Figure1-1

对于图片的解释可看这篇文章EM算法实例及python实现

EM算法还有个常用应用是在GMM(高斯混合模型,Gaussian misture model),《统计学习方法》9.3节对于高斯混合模型的EM算法做了很详细的介绍,主要是用EM算法来估计GMM的参数;GMM在生信信息中也比较常见,GATK的VQSR就用到了GMM,其通过过GMM模型构建一个分类器来对变异数据进行打分,从而评估每个位点的可行度

因此以《统计学习方法》习题9.3来了解下EM算法是如何实现GMM参数估计的(公式推导可看书。。。),以Python实现:

# 公式9.25,计算高斯分布密度
def gaussian_pro(x, mu, sigma):
    gsi = np.exp(-1 * (x - mu) * (x - mu) / (2 * sigma**2)) * \
    (1 / (math.sqrt(2 * math.pi) * sigma))
    return(gsi)

# 算法9.2的E步,计算分模型k对观测数据yi的响应度gamma
def E_step(x, alpha_s, mu_s, sigma_s, K):
    gamma = []
    for i in range(K): 
        gamma.append(alpha_s[i] * gaussian_pro(x, mu_s[i], sigma_s[i]))

    sum_gamma = sum(gamma)

    for i in range(K): 
        gamma[i] /= sum_gamma
    return(gamma)

# 算法9.2的M步,更新模型参数mu,sigma,alpha
def M_step(x, mu_s, gamma, K):
    mu_k = []
    sigma_k = []
    alpha_k = []

    for i in range(K): 
        # 更新 mu
        mu_k.append(np.dot(gamma[i], x) / np.sum(gamma[i]))
        # 更新 sigam
        sigma_k.append(math.sqrt(np.dot(gamma[i], (x - mu_s[i])**2) / np.sum(gamma[i])))
        # 更新 alpha
        alpha_k.append(np.sum(gamma[i]) / len(gamma[i]))
    return mu_k, sigma_k, alpha_k

# 迭代E步和M步,直至收敛
def EM(x, K, max_iter = 500):
    # 初始化alpha, mu, sigma
    alpha = [1 / K for i in range(K)]
    random.seed(12345)
    mu = [random.random() for i in range(K)]
    sigma = [random.random()*10 for i in range(K)]

    gamma = np.zeros((K, len(x)))
    iter = 0
    while iter < max_iter:
        gamma_old = copy.deepcopy(gamma)
        gamma = E_step(x, alpha, mu, sigma, K)
        mu, sigma, alpha = M_step(x, mu, gamma, K)
        # gamma变化小于1e-10则停止迭代
        if np.sum(abs(gamma[0] - gamma_old[0])) < 1e-10:
            break
        else:
            iter += 1
    return alpha, mu, sigma

按照习题例子输入观测数据(这里是简单的一维数据),通过EM函数估计两个分量的高斯混合模型的参数

import numpy as np
import math
import random
import copy

alpha, mu, sigma = EM(dataset, 2)
print("alpha1:%.3f, alpha2:%.3f" %(alpha[0], alpha[1]))
print("mu1:%.1f, mu2:%.1f" %(mu[0], mu[1]))
print("sigma1:%.1f, sigma2:%.1f" %(sigma[0], sigma[1]))
# alpha1:0.658, alpha2:0.342
# mu1:21.6, mu2:19.7
# sigma1:44.5, sigma2:8.3 

参考资料:
Python AI 教学 | EM算法(Expectation Maximization Algorithm)及应用
机器学习系列-EM算法
EM算法原理总结

本文出自于http://www.bioinfo-scrounger.com转载请注明出处