期望最大化(Expectation Maximization)硬币抛掷示例

10

最近我一直在自学期望最大化算法(Expectation Maximization),并通过一些简单的例子加深了对它的理解:

http://cs.dartmouth.edu/~cs104/CS104_11.04.22.pdf 有三枚硬币,编号为0、1、2,分别有P0、P1和P2的概率掷出正面。掷硬币0,如果结果是正面,则抛3次硬币1;否则,抛3次硬币2。硬币1和2掷出的数据如下:HHH、TTT、HHH、TTT、HHH。隐藏数据是硬币0的结果。估计P0、P1和P2。

http://ai.stanford.edu/~chuongdo/papers/em_tutorial.pdf 有两枚硬币A和B,分别有PA和PB的概率掷出正面。每轮随机选择一个硬币抛10次,然后记录结果。观测数据是这两枚硬币的投掷结果。但我们不知道每轮选择了哪个硬币。估计PA和PB。

虽然我能够进行计算,但我不知道它们是如何与原始EM理论相关联的。具体来说,在两个例子的M步骤中,我不明白它们如何最大化任何东西。它似乎只是重新计算参数,而新的参数比旧的更好。此外,这两个E步骤甚至看起来都不相似,更不用说原始理论的E步骤了。

那么这些例子到底是如何工作的呢?


2
"可能在计算机科学领域,http://cs.stackexchange.com/ 是一个更好的提问平台。" - hatchet - done with SOverflow
4个回答

8
第二个PDF文件无法下载,但我也访问了维基百科页面http://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm,其中有更多信息。http://melodi.ee.washington.edu/people/bilmes/mypapers/em.pdf(声称是温和的介绍)也值得一看。
EM算法的整个重点是找到最大化观察数据似然的参数。这是第一个PDF第8页唯一的要点,大写希腊字母Theta下标ML的方程式。
EM算法在存在隐藏数据的情况下非常有用,如果您知道隐藏数据,问题就会变得容易。在三枚硬币的例子中,这是抛掷硬币0的结果。如果您知道其结果,您可以(当然)估计硬币0正面朝上的概率。您还将知道在下一阶段中抛掷了硬币1或硬币2三次的情况,这将允许您对硬币1和硬币2正面朝上的概率进行估计。这些估计将被证明是最大化观察到的数据的似然性,其中不仅包括您获得的结果,还包括您没有获得的隐藏数据-来自硬币0的结果。对于一个正面朝上A次,反面朝上B次的硬币,您会发现正面朝上的概率的最大似然值为A /(A + B)-您可能值得详细计算此值,因为它是M步骤的基石。
在EM算法中,你说尽管你不知道隐藏数据,但可以使用概率估计来为它写出概率分布。对于每个可能的隐藏数据值,你都可以找到能够优化包括隐藏数据的数据的对数似然函数的参数值,而这几乎总是意味着要计算某种加权平均值(如果不是这样,EM步骤可能就太难以实现了)。
EM算法要求你找到最大化所有可能的隐藏数据值所给出的与观测数据相关的概率的加权对数似然之和的参数,其中权重由EM步骤开始时使用参数的情况下给定的相关隐藏数据的概率确定。这就是几乎所有人,包括维基百科算法,所称的Q函数。在维基百科文章中给出了EM算法的证明,该证明表明,如果改变参数以增加Q函数(这只是一种手段),则也将改变它们以增加观察数据的似然性(这是你关心的)。在实践中,你往往会发现可以使用隐藏数据的概率(在EM步骤开始时给出的估计)以某种方式加权观察结果来最大化Q函数。
在你的例子中,它意味着将每个硬币产生的正面和反面数量相加。在PDF中,他们计算出P(Y=H|X=) = 0.6967。这意味着您使用重量0.6967的情况Y=H,这意味着您通过0.6967增加Y=H的计数,并通过3*0.6967增加硬币1中X=H的计数,您通过0.3033增加Y=T的计数,并通过3*0.3033增加硬币2中X=H的计数。如果您对为什么A/(A+B)是标准情况下硬币概率的最大似然有详细的证明,那么您应该准备将其转化为为什么这种加权更新方案最大化Q函数的证明。
最后,观察到的数据的对数似然(您正在最大化的内容)可以为您提供非常有用的检查。它应该随着每个EM步骤而增加,至少直到您接近收敛时会出现舍入误差,在这种情况下,您可能会有一个非常小的减少,表示收敛。如果它急剧下降,则说明您的程序或数学存在错误。

7
侥幸的是,我最近也一直在苦恼这个材料。以下是我对它的理解:
考虑一个相关但不同的算法,称为classify-maximize算法,我们可以将其用作混合模型问题的解决技术。混合模型问题是指我们有一系列可能由N个不同过程中的任意一个产生的数据,我们知道这些过程的一般形式(例如高斯),但我们不知道过程的参数(例如均值和/或方差),甚至可能不知道各个过程的相对可能性。 (通常,我们至少知道过程的数量。如果没有这个,我们就进入了所谓的“非参数”领域。)从某种意义上说,生成每个数据的过程是该问题的“缺失”或“隐藏”数据。
现在,这个相关的classify-maximize算法所做的是从一些任意的过程参数猜测开始。每个数据点根据每个参数过程进行评估,并生成一组概率-数据点由第一个过程、第二个过程等等生成的概率,一直到最后一个第N个过程。然后,每个数据点根据最可能的过程进行分类。
此时,我们已经将数据分成了N个不同的类别。因此,对于每个数据类别,我们可以使用相对简单的微积分技术来优化该聚类的参数,使用最大似然技术进行优化。(如果我们在分类之前尝试对整个数据集这样做,通常会遇到解析上不可行的问题。)
然后我们更新我们的参数猜测值,重新分类,更新参数,重新分类等等,直到收敛。
期望最大化算法所做的是类似的,但更加通用:现在,我们使用软分类,其中每个数据点属于每个过程的概率都有一些可能性。(显然,每个点的概率需要总和为1,因此存在一些归一化处理。)我认为我们也可以认为每个过程/猜测对每个数据点具有一定的“解释力”。
因此,现在我们不是针对绝对属于每个类别的点(忽略绝对不属于某些点)优化猜测,而是在这些软分类或解释力的背景下重新优化猜测。如果您以正确的方式编写表达式,则您正在最大化其形式为期望值的函数。
话虽如此,还有一些注意事项:

1) 这听起来很容易,但对我来说并非如此。文献中充斥着各种特殊技巧和技术 - 使用似然表达式而不是概率表达式,转换为对数似然,使用指示变量,将它们放在基向量形式和将它们放在指数中等等。

这些技巧对于你已经掌握了一般思路可能更有帮助,但它们也可能使核心思想变得复杂难懂。

2) 无论问题的约束条件是什么,都可能很难纳入到框架中。尤其是,如果你知道每个过程的概率,则你可能做得很好。如果不知道,则需要估计这些概率,并且各个过程的概率之和必须为一;它们必须存在于一个概率单纯形上。如何保持这些约束不变并不总是明显的。

3) 这是一种非常通用的技术,我不知道如何编写通用代码。应用远远超出简单的聚类,并扩展到许多实际缺失数据的情况下,或者缺失数据的假设可能对你有所帮助。对于许多应用来说,这里有一种恶魔般的独创性。

4) 这种技术被证明能收敛,但并不一定是全局最优;需谨慎。

我发现以下链接对于上述解释很有帮助:统计学习幻灯片

而以下的文章详细介绍了一些痛苦的数学细节:迈克尔·科林斯的论文


3
我用Python写了下面的代码,解释了Do和Batzoglou在第二个示例论文中给出的例子。
我建议您先阅读链接,以清楚地解释为什么在下面的代码中获得'weightA''weightB'免责声明: 该代码确实有效,但我确定它没有被编码优化。我通常不是Python程序员,两周前开始使用它。
import numpy as np
import math

#### E-M Coin Toss Example as given in the EM tutorial paper by Do and Batzoglou* #### 

def get_mn_log_likelihood(obs,probs):
    """ Return the (log)likelihood of obs, given the probs"""
    # Multinomial Distribution Log PMF
    # ln (pdf)      =             multinomial coeff            *   product of probabilities
    # ln[f(x|n, p)] = [ln(n!) - (ln(x1!)+ln(x2!)+...+ln(xk!))] + [x1*ln(p1)+x2*ln(p2)+...+xk*ln(pk)]     

    multinomial_coeff_denom= 0
    prod_probs = 0
    for x in range(0,len(obs)): # loop through state counts in each observation
        multinomial_coeff_denom = multinomial_coeff_denom + math.log(math.factorial(obs[x]))
        prod_probs = prod_probs + obs[x]*math.log(probs[x])

multinomial_coeff = math.log(math.factorial(sum(obs))) -  multinomial_coeff_denom
likelihood = multinomial_coeff + prod_probs
return likelihood

# 1st:  Coin B, {HTTTHHTHTH}, 5H,5T
# 2nd:  Coin A, {HHHHTHHHHH}, 9H,1T
# 3rd:  Coin A, {HTHHHHHTHH}, 8H,2T
# 4th:  Coin B, {HTHTTTHHTT}, 4H,6T
# 5th:  Coin A, {THHHTHHHTH}, 7H,3T
# so, from MLE: pA(heads) = 0.80 and pB(heads)=0.45

# represent the experiments
head_counts = np.array([5,9,8,4,7])
tail_counts = 10-head_counts
experiments = zip(head_counts,tail_counts)

# initialise the pA(heads) and pB(heads)
pA_heads = np.zeros(100); pA_heads[0] = 0.60
pB_heads = np.zeros(100); pB_heads[0] = 0.50

# E-M begins!
delta = 0.001  
j = 0 # iteration counter
improvement = float('inf')
while (improvement>delta):
    expectation_A = np.zeros((5,2), dtype=float) 
    expectation_B = np.zeros((5,2), dtype=float)
    for i in range(0,len(experiments)):
        e = experiments[i] # i'th experiment
        ll_A = get_mn_log_likelihood(e,np.array([pA_heads[j],1-pA_heads[j]])) # loglikelihood of e given coin A
        ll_B = get_mn_log_likelihood(e,np.array([pB_heads[j],1-pB_heads[j]])) # loglikelihood of e given coin B

        weightA = math.exp(ll_A) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of A proportional to likelihood of A 
        weightB = math.exp(ll_B) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of B proportional to likelihood of B                            

        expectation_A[i] = np.dot(weightA, e) 
        expectation_B[i] = np.dot(weightB, e)

    pA_heads[j+1] = sum(expectation_A)[0] / sum(sum(expectation_A)); 
    pB_heads[j+1] = sum(expectation_B)[0] / sum(sum(expectation_B)); 

    improvement = max( abs(np.array([pA_heads[j+1],pB_heads[j+1]]) - np.array([pA_heads[j],pB_heads[j]]) ))
    j = j+1

1
理解这个的关键是知道辅助变量是什么,使得估计变得微不足道。我将快速解释第一个示例,第二个示例遵循类似的模式。
将每个正面/反面序列与两个二进制变量相结合,指示是否使用硬币1或硬币2。现在我们的数据如下所示:
c_11 c_12 c_21 c_22 c_31 c_32 ...
对于每个i,要么c_i1=1,要么c_i2=1,另一个为0。如果我们知道这些变量在样本中的取值,参数的估计将是微不足道的:p1将是在c_i1=1的样本中头的比例,同样适用于c_i2,而\ lambda将是c_i1s的平均值。
然而,我们不知道这些二进制变量的值。因此,我们基本上是猜测它们(实际上,取它们的期望值),然后在假设我们的猜测是正确的情况下更新模型中的参数。因此,E步骤是取c_i1和c_i2的期望值。M步骤是在给定这些c的情况下最大化似然估计p_1、p_2和\ lambda。

这样是否更容易理解一些?如果您愿意,我可以写出E步骤和M步骤的更新。EM算法保证了在迭代增加的过程中,通过遵循此过程,似然性永远不会降低。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接