期望最大化抛硬币示例

Question

我最近一直在自学期望最大化，并在这个过程中抓住了一些简单的例子：

http://cs.dartmouth.edu/~cs104/CS104_11.04.22.pdf 有 3 个硬币 0、1 和 2，投掷时有 P0、P1 和 P2 概率落在正面。掷硬币 0，如果结果是正面，则掷硬币 1 3 次，否则掷硬币 2 3 次。硬币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-Step 期间，我看不到它们如何最大化任何东西。似乎他们正在重新计算参数，并且不知何故，新参数比旧参数更好。而且，这两个E-Step甚至看起来都不相似，更不用说原始理论的E-Step了。

那么这些示例究竟是如何工作的呢？

Answer 1

第二个 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 在下一阶段是否被抛了 3 次，这将允许您估计硬币 1 和硬币 2 出现正面的概率。这些估计是合理的，因为它们最大化了观察数据的可能性，这不仅包括给你的结果，还包括你不是的隐藏数据——硬币 0 的结果。 A 正面和 B 反面，您会发现 A 正面概率的最大可能性是 A/(A+B) - 这可能值得您详细计算，因为它是 M 步骤的构建块。

在 EM 算法中，您说虽然您不知道隐藏的数据，但您会得到概率估计，从而可以为它写下概率分布。对于隐藏数据的每个可能值，您可以找到可以优化包括隐藏数据在内的数据的对数似然的参数值，这几乎总是意味着计算某种加权平均值（如果它不是 EM步骤可能太难而不实用）。

EM 算法要求您做的是找到使所有可能的隐藏数据值给出的对数似然的加权和最大化的参数，其中权重由给定观察值的相关隐藏数据的概率给出，使用EM 步骤开始时的参数。这就是几乎所有人，包括维基百科算法，都称之为 Q 函数。维基百科文章中给出的 EM 算法背后的证明说，如果你改变参数以增加 Q 函数（这只是达到目的的一种手段），你也会改变它们以增加观察到的数据的可能性（您确实关心）。您在实践中往往会发现，如果您知道隐藏数据，您可以使用隐藏数据的概率来最大化 Q 函数，但使用隐藏数据的概率，给定 EM 开始时的估计步骤，以某种方式对观察结果进行加权。

在您的示例中，这意味着将每个硬币产生的正面和反面的数量相加。在 PDF 中，他们计算出 P(Y=H|X=) = 0.6967。这意味着您在 Y=H 的情况下使用权重 0.6967，这意味着您将 Y=H 的计数增加 0.6967，并将硬币 1 中 X=H 的计数增加 3*0.6967，并且您将 Y 的计数增加=T 增加 0.3033 并将硬币 2 中 X=H 的计数增加 3*0.3033。如果您对为什么 A/(A+B) 是标准情况下硬币概率的最大似然有详细的理由，您应该准备好将其转化为为什么这种加权更新方案最大化 Q 函数的理由。

最后，观察数据的对数似然（您正在最大化的东西）为您提供了一个非常有用的检查。它应该随着每个 EM 步骤而增加，至少直到您非常接近收敛以至于出现舍入误差，在这种情况下，您可能会有非常小的减少，表明收敛。如果它急剧下降，则说明您的程序或数学存在错误。

Answer 2

理解这一点的关键是了解使估计变得微不足道的辅助变量。我将快速解释第一个示例，第二个遵循类似的模式。

用两个二进制变量增加每个正面/反面序列，这表明使用的是硬币 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_i1s和c_i2s的期望值。 M 步骤是在给定这些 cs 的情况下对 p_1、p_2 和 λ 进行最大似然估计。

这是否更有意义？如果您愿意，我可以写出 E 和 M 步骤的更新。然后 EM 只是保证通过遵循这个过程，可能性永远不会随着迭代的增加而降低。

Answer 3

幸运的是，我最近也一直在为这种材料而苦苦挣扎。以下是我的想法：

考虑一种相关但不同的算法，称为分类最大化算法，我们可以将其用作混合模型问题的解决技术。混合模型问题是这样一个问题：我们有可能由 N 个不同过程中的任何一个产生的数据序列，我们知道这些数据的一般形式（例如，高斯），但我们不知道过程的参数（例如，均值和/或方差），甚至可能不知道过程的相对可能性。 （通常我们至少知道过程的数量。没有它，我们就进入了所谓的“非参数”领域。）从某种意义上说，生成每个数据的过程是“丢失”或“隐藏”数据的问题。

现在，这个相关的分类最大化算法所做的是从对过程参数的一些任意猜测开始。根据这些参数过程中的每一个对每个数据点进行评估，并生成一组概率——数据点由第一个过程、第二个过程等生成的概率，直到最后的第 N 个过程。然后根据最可能的过程对每个数据点进行分类。

此时，我们将数据分为 N 个不同的类。因此，对于每个类数据，我们可以通过一些相对简单的演算，使用最大似然技术优化该集群的参数。 （如果我们试图在分类之前对整个数据集执行此操作，通常在分析上是难以处理的。）

然后我们更新我们的参数猜测，重新分类，更新我们的参数，重新分类等，直到收敛。

期望最大化算法的作用类似，但更通用：我们现在使用软分类，而不是将数据点硬分类为第 1 类、第 2 类、... 到第 N 类，点以一定的概率属于每个过程。 （显然，每个点的概率总和需要为 1，因此需要进行一些标准化。）我认为我们也可以将其视为每个过程/猜测对每个数据都有一定的“解释力”点。

所以现在，我们不是针对绝对属于每个类别的点优化猜测（忽略绝对不属于的点），而是在那些软分类或解释能力的上下文中重新优化猜测。碰巧的是，如果您以正确的方式编写表达式，那么您要最大化 的是一个函数，它的形式是一个期望。

话虽如此，但有一些警告：

1) 这听起来很简单。至少对我来说不是。文献中充斥着各种特殊技巧和技术——使用似然表达式而不是概率表达式，转换为对数似然，使用指示变量，将它们置于基向量形式并将它们置于指数中，等等。

一旦您有了大致的想法，这些可能会更有帮助，但它们也会混淆核心想法。

2) 无论您对问题有什么限制，都很难将其整合到框架中。特别是，如果您知道每个过程的概率，您可能处于良好状态。如果不是，您也在估计这些，并且过程的概率之和必须为 1；他们必须生活在概率单纯形上。如何保持这些约束完好无损并不总是显而易见的。

3) 这是一种足够通用的技术，我不知道如何编写通用代码。这些应用程序远远超出了简单的聚类范围，并扩展到了许多实际丢失数据的情况，或丢失数据的假设可能对您有所帮助。对于许多应用程序来说，这里有一种可怕的独创性。

4) 这种技术被证明是收敛的，但收敛不一定是全局最大值；小心点。

我发现以下链接有助于提出上述解释：Statistical learning slides

以下文章详细介绍了一些令人痛苦的数学细节：Michael Collins' write-up

Answer 4

我用 Python 编写了以下代码，解释了 Do 和 Batzoglou 在您的 second example paper 中给出的示例。

我建议您先阅读此link，以清楚地了解以下代码中的'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