KeepNotes blog

Stay hungry, Stay Foolish.

0%

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

Posted on Edited on In
Symbols count in article: 2.5k Reading time ≈ 2 mins.

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转载请注明出处

Related Posts