这是一个通过实用且（在我看来）非常直观的“投币”示例来学习 EM 的秘诀：

1. 阅读 Do 和 Batzoglou 的这篇简短的EM 教程论文。这是解释抛硬币示例的模式：

   

2. 您的脑海中可能会有问号，尤其是关于期望步骤中的概率来自何处。请查看此数学堆栈交换页面上的说明。

3. 查看/运行我在 Python 中编写的代码，它模拟了项目 1 的 EM 教程论文中抛硬币问题的解决方案：

   import numpy as np
   import math
   import matplotlib.pyplot as plt
   
   ## E-M Coin Toss Example as given in the EM tutorial paper by Do and Batzoglou* ##
   
   def get_binomial_log_likelihood(obs,probs):
       """ Return the (log)likelihood of obs, given the probs"""
       # Binomial Distribution Log PDF
       # ln (pdf)      = Binomial Coeff * product of probabilities
       # ln[f(x|n, p)] =   comb(N,k)    * num_heads*ln(pH) + (N-num_heads) * ln(1-pH)
   
       N = sum(obs);#number of trials  
       k = obs[0] # number of heads
       binomial_coeff = math.factorial(N) / (math.factorial(N-k) * math.factorial(k))
       prod_probs = obs[0]*math.log(probs[0]) + obs[1]*math.log(1-probs[0])
       log_lik = binomial_coeff + prod_probs
   
       return log_lik
   
   # 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((len(experiments),2), dtype=float) 
       expectation_B = np.zeros((len(experiments),2), dtype=float)
       for i in range(0,len(experiments)):
           e = experiments[i] # i'th experiment
             # loglikelihood of e given coin A:
           ll_A = get_binomial_log_likelihood(e,np.array([pA_heads[j],1-pA_heads[j]])) 
             # loglikelihood of e given coin B
           ll_B = get_binomial_log_likelihood(e,np.array([pB_heads[j],1-pB_heads[j]])) 
   
             # corresponding weight of A proportional to likelihood of A 
           weightA = math.exp(ll_A) / ( math.exp(ll_A) + math.exp(ll_B) ) 
   
             # corresponding weight of B proportional to likelihood of B
           weightB = math.exp(ll_B) / ( math.exp(ll_A) + math.exp(ll_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
   
   plt.figure();
   plt.plot(range(0,j),pA_heads[0:j], 'r--')
   plt.plot(range(0,j),pB_heads[0:j])
   plt.show()
   

听起来您的问题有两个部分：基本思想和具体示例。我将从基本思想开始，然后链接到底部的示例。

EM 在您似乎需要知道的 Catch-22 情况下很有用$A$在你计算之前$B$你需要知道$B$在你计算之前$A$.

人们处理的最常见的情况可能是混合分布。对于我们的示例，让我们看一个简单的高斯混合模型：

您有两个具有不同均值和单位方差的不同单变量高斯分布。

你有一堆数据点，但你不确定哪些点来自哪个分布，你也不确定这两个分布的均值。

现在你被困住了：

• 如果你知道真正的手段，你可以弄清楚哪些数据点来自哪个高斯。例如，如果一个数据点的值非常高，它可能来自具有较高均值的分布。但是你不知道手段是什么，所以这行不通。

• 如果您知道每个点来自哪个分布，那么您可以使用相关点的样本均值来估计两个分布的均值。但是您实际上并不知道将哪些点分配给哪个分布，所以这也行不通。

因此，这两种方法似乎都行不通：您需要先知道答案，然后才能找到答案，而您却被卡住了。

EM 让您做的是在这两个易于处理的步骤之间交替，而不是一次处理整个过程。

您需要从对这两种方法的猜测开始（尽管您的猜测不一定非常准确，但您确实需要从某个地方开始）。

如果您对方法的猜测是准确的，那么您将有足够的信息来执行我上面第一个要点中的步骤，并且您可以（概率地）将每个数据点分配给两个高斯中的一个。即使我们知道我们的猜测是错误的，我们还是试试这个。然后，给定每个点的分配分布，您可以使用第二个要点获得均值的新估计值。事实证明，每次循环执行这两个步骤时，都会提高模型可能性的下限。

这已经很酷了：尽管上面要点中的两个建议看起来不像单独工作，但您仍然可以将它们一起使用来改进模型。EM的真正魔力在于，经过足够多的迭代后，下限将如此之高，以至于它与局部最大值之间不会有任何空间。结果，您已经在本地优化了可能性。

因此，您不仅改进了模型，还找到了可以通过增量更新找到的最佳模型。

维基百科的这个页面显示了一个稍微复杂的例子（二维高斯和未知协方差），但基本思想是相同的。它还包括R用于实现示例的注释良好的代码。

在代码中，“期望”步骤（E-step）对应于我的第一个要点：在给定每个高斯的当前参数的情况下，找出哪个高斯负责每个数据点。给定这些分配，“最大化”步骤（M-step）更新均值和协方差，如我的第二个要点。

正如您在动画中看到的那样，这些更新很快使算法从一组糟糕的估计变为一组非常好的估计：似乎确实有两团以 EM 发现的两个高斯分布为中心的点云。

这是用于估计均值和标准差的期望最大化 (EM) 示例。代码用 Python 编写，但即使您不熟悉该语言，也应该很容易理解。

EM的动机

下面显示的红点和蓝点来自两个不同的正态分布，每个正态分布都有特定的平均值和标准差：

为了计算红色分布的“真实”均值和标准差参数的合理近似值，我们可以很容易地查看红色点并记录每个点的位置，然后使用熟悉的公式（蓝色组也是如此） .

现在考虑这样一种情况，我们知道有两组点，但我们看不到哪个点属于哪个组。换句话说，颜色是隐藏的：

如何将点分成两组并不明显。我们现在不能只查看位置并计算红色分布或蓝色分布的参数的估计值。

这就是EM可以用来解决问题的地方。

使用 EM 估计参数

这是用于生成上述点的代码。您可以看到从中提取点的正态分布的实际均值和标准差。变量red和blue分别保存红色组和蓝色组中每个点的位置：

import numpy as np
from scipy import stats

np.random.seed(110) # for reproducible random results

# set parameters
red_mean = 3
red_std = 0.8

blue_mean = 7
blue_std = 2

# draw 20 samples from normal distributions with red/blue parameters
red = np.random.normal(red_mean, red_std, size=20)
blue = np.random.normal(blue_mean, blue_std, size=20)

both_colours = np.sort(np.concatenate((red, blue)))

如果我们能看到每个点的颜色，我们将尝试使用库函数恢复均值和标准差：

>>> np.mean(red)
2.802
>>> np.std(red)
0.871
>>> np.mean(blue)
6.932
>>> np.std(blue)
2.195

但由于颜色对我们隐藏，我们将开始 EM 过程......

首先，我们只是猜测每组参数的值（步骤 1）。这些猜测不一定是好的：

# estimates for the mean
red_mean_guess = 1.1
blue_mean_guess = 9

# estimates for the standard deviation
red_std_guess = 2
blue_std_guess = 1.7

非常糟糕的猜测——这些方法看起来离一组点的任何“中间”都很远。

为了继续使用 EM 并改进这些猜测，我们计算每个数据点（无论其秘密颜色）出现在这些猜测下的均值和标准差的可能性（步骤 2）。

该变量both_colours保存每个数据点。该函数stats.norm使用给定参数计算点在正态分布下的概率：

likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)
likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)

例如，这告诉我们，根据我们目前的猜测，1.761 处的数据点更可能是红色（0.189）而不是蓝色（0.00003）。

我们可以将这两个似然值转换为权重（步骤 3），使它们总和为 1，如下所示：

likelihood_total = likelihood_of_red + likelihood_of_blue

red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total

使用我们当前的估计和新计算的权重，我们现在可以计算新的、可能更好的参数估计（步骤 4）。我们需要一个均值函数和一个标准差函数：

def estimate_mean(data, weight):
    return np.sum(data * weight) / np.sum(weight)

def estimate_std(data, weight, mean):
    variance = np.sum(weight * (data - mean)**2) / np.sum(weight)
    return np.sqrt(variance)

这些看起来与数据平均值和标准差的常用函数非常相似。不同之处在于使用了weight为每个数据点分配权重的参数。

这种加权是 EM 的关键。数据点上颜色的权重越大，数据点对该颜色参数的下一次估计的影响就越大。最终，这具有将每个参数拉向正确方向的效果。

使用以下函数计算新的猜测：

# new estimates for standard deviation
blue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)
red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)

# new estimates for mean
red_mean_guess = estimate_mean(both_colours, red_weight)
blue_mean_guess = estimate_mean(both_colours, blue_weight)

然后用这些新的猜测从第 2 步开始重复 EM 过程。我们可以对给定的迭代次数（比如 20 次）重复这些步骤，或者直到我们看到参数收敛。

经过五次迭代后，我们看到我们最初的错误猜测开始变得更好：

经过 20 次迭代后，EM 过程或多或少地收敛了：

为了比较，这里是 EM 过程的结果与未隐藏颜色信息的计算值的比较：

          | EM guess | Actual 
----------+----------+--------
Red mean  |    2.910 |   2.802
Red std   |    0.854 |   0.871
Blue mean |    6.838 |   6.932
Blue std  |    2.227 |   2.195

注意：这个答案改编自我在 Stack Overflow 上的回答。

按照 Zhubarb 的回答，我在 GNU R 中实现了 Do 和 Batzoglou“抛硬币”EM 示例。请注意，我使用了包的mle功能stats4- 这有助于我更清楚地了解 EM 和 MLE 的关系。

require("stats4");

## sample data from Do and Batzoglou
ds<-data.frame(heads=c(5,9,8,4,7),n=c(10,10,10,10,10),
    coin=c("B","A","A","B","A"),weight_A=1:5*0)

## "baby likelihood" for a single observation
llf <- function(heads, n, theta) {
  comb <- function(n, x) { #nCr function
    return(factorial(n) / (factorial(x) * factorial(n-x)))
  }
  if (theta<0 || theta >1) { # probabilities should be in [0,1]
    return(-Inf);
  }
  z<-comb(n,heads)* theta^heads * (1-theta)^(n-heads);
  return (log(z))
}

## the "E-M" likelihood function
em <- function(theta_A,theta_B) {
  # expectation step: given current parameters, what is the likelihood
  # an observation is the result of tossing coin A (vs coin B)?
  ds$weight_A <<- by(ds, 1:nrow(ds), function(row) {
    llf_A <- llf(row$heads,row$n, theta_A);
    llf_B <- llf(row$heads,row$n, theta_B);

    return(exp(llf_A)/(exp(llf_A)+exp(llf_B)));
  })

  # maximisation step: given params and weights, calculate likelihood of the sample
  return(- sum(by(ds, 1:nrow(ds), function(row) {
    llf_A <- llf(row$heads,row$n, theta_A);
    llf_B <- llf(row$heads,row$n, theta_B);

    return(row$weight_A*llf_A + (1-row$weight_A)*llf_B);
  })))
}

est<-mle(em,start = list(theta_A=0.6,theta_B=0.5), nobs=NROW(ds))