期望最大化（EM）算法被广泛用于估计不同统计模型的参数。它是一种迭代算法，可以将一个困难的优化问题分解为几个更简单的优化问题。本文通过几个简单的例子解释了它是如何工作的。这个算法最流行的例子（也是互联网上讨论最多的）可能来自这篇论文（http://www.nature.com/nbt/jou...）。这是一个非常简单的示例，所以让我们也从这里开始。假设我们有两枚硬币（硬币1和硬币2），正面朝上的概率不同。我们选择其中一枚硬币，翻转m=10并记录正面朝上的次数。假设我们重复这个实验n=5次。我们的任务是确定每枚硬币正面朝上的概率。我们有：首先假设我们知道在每个实验中使用了哪种硬币。在这种情况下，在信息完整的情况下，可以使用最大似然估计(MLE)技术轻松求解p_1和p_2。首先计算似然函数并取其对数（因为更容易最大化对数似然函数）。由于我们有n个独立的实验，似然函数只是在x_i处评估的各个概率质量函数(PMF)的乘积（实验i中的数字为正数）。现在我们需要最大化关于概率p_1和p_2的对数似然函数。它可以通过数值或分析来完成。我将演示这两种方法。我们先尝试求解，可以分别求解p_1和p_2。对p_1取对数似然函数的导数，将其设置为零并求解p_1。在对数似然函数求微分时，p_2项的导数会等于0。所以我们只使用硬币1的实验数据。得到的答案很直观：就是我们得到的硬币正面朝上的总数1个实验除以硬币1个实验中的翻转总数。p_2的计算将类似。现在我将在Python中实现这个解决方案。m=10#experimentn=5#experimentsnumbersp_1=0.8p_2=0.3xs=[]#eachexperimentzs=[0,0,1,0,1]#flipnp.random的硬币数量.seed(5)foriinzs:ifi==0:xs.append(stats.binom(n=m,p=p_1).rvs())#flipcoin1elifi==1:xs.append(stats.binom(n=m,p=p_2).rvs())#flipcoin2xs=np.array(xs)print(xs)exp1=[0,1,3]#coin1exp2的实验索引=[2,4]#硬币实验索引2print('解析解:')print('p1:',xs[exp1].sum()/(m*len(exp1)))print('p2:',xs[exp2].sum()/(m*len(exp2)))这些实现比较简单。我们只是实现上面计算的公式并插入数字。运行此代码的结果如下所示。解决此问题的另一种方法是使用数值求解器。我们需要找到最大化对数似然函数的解决方案，当使用数值求解器时，无需计算导数和手动求解最大化对数似然函数的参数。只需实现一个我们想要最大化的函数并将其传递给数值求解器即可。由于Python中的大多数求解器旨在最小化给定函数，因此我们实现了一个计算负对数似然的函数（因为最小化负对数似然与最大化对数似然相同）。代码和结果如下所示。efneg_log_likelihood(probs,m,xs,zs):'''计算负对数似然'''ll=0forx,zinzip(xs,zs):ll+=stats.binom(p=probs[z],n=m).logpmf(x)return-llres=optimize.minimize(neg_log_likelihood,[0.5,0.5],bounds=[(0,1),(0,1)],args=(m,xs,zs),method='tnc')print('数值解:')print('p1:',res.x[0])print('p2:',res.x[1])和上面一个方法产生与解析解相同的结果。现在假设我们不知道每个实验中使用的是哪种硬币。在这种情况下，未观察到变量Z（称为潜在或隐藏变量）并且数据集变得不完整。估计概率p_1和p_2现在更加困难，但仍然可以借助EM算法完成。如果您知道选择硬币1或硬币2的概率，则可以使用贝叶斯定理来估计每个硬币的偏差。如果您知道每个硬币的偏差，则可以估计在给定实验中使用硬币1或硬币2的概率。在EM算法中，我们对这些概率进行初始猜测，然后在两个步骤之间迭代（在给定偏差的情况下估计使用每个硬币的概率和在给定偏差的情况下估计使用每个硬币的概率）直到收敛。可以从数学上证明，这种算法收敛于（似然函数的）局部最大值。现在尝试复制论文中提供的示例。问题设置如下。完整数据集由两个变量X和Z组成，但仅观察到X。由于随机选择了两个硬币之一，因此选择每个硬币的概率等于0.5。为了求解theta，我们需要最大化以下似然函数：上述函数称为不完全似然函数。如果我们计算它的对数，我们得到：在对数下求和使得最大化问题变得困难。让我们在似然函数中包含隐变量Z以获得全似然：全似然函数的对数为：这样就没有对数内求和，更容易解决这个函数的最大化问题。由于未观察到Z值，因此无法直接计算和最大化此函数。但是如果我们知道Z的分布，我们可以计算它的期望值并用它来最大化似然（或对数似然）函数。这里的问题是要找到Z的分布需要知道参数theta，而要找到参数theta需要知道Z的分布。EM算法可以让我们解决这个问题。我们从对参数theta的随机猜测开始，然后迭代以下步骤：期望步骤（E-step）：在给定当前条件分布数据X和当前参数估计theta步骤（M-step）：找到最大化此期望的参数theta的值贝叶斯定理可用于找到给定X_i和theta的Z_i的条件分布：现在定义条件期望完整的对数似然如下：插入完整的对数似然函数并重新排列：这完成了E步。在M步中，我们最大化上面计算的关于参数theta的函数：在这个例子中，可以通过分析找到参数（通过取导数并将其设置为零）。此处不演示整个过程，但提供了答案。在下面的实现中，我们将使用与论文中相同的数据来检查是否获得了相同的结果。Python代码如下m=10#每个样本的翻转次数=5#样本数xs=np.array([5,9,8,4,7])theta=[0.6,0.5]#初始猜测p_1,p_2foriinrange(10):p_1,p_2=thetaT_1s=[]T_2s=[]#e-stepforxinxs:T_1=stats.binom(n=m,p=theta[0]).pmf(x)/(stats.binom(n=m,p=theta[0]).pmf(x)+stats.binom(n=m,p=theta[1]).pmf(x))T_2=stats.binom(n=m,p=theta[1]).pmf(x)/(stats.binom(n=m,p=theta[0]).pmf(x)+stats.binom(n=m,p=theta[1]).pmf(x))T_1s.append(T_1)T_2s.append(T_2)#M步T_1s=np.array(T_1s)T_2s=np.array(T_2s)p_1=np.dot(T_1s,xs)/(m*np.sum(T_1s))p_2=np.dot(T_2s,xs)/(m*np.sum(T_2s))print(f'step:{i},p1={p_1:.2f},p2={p_2:.2f}')theta=[p_1,p_2]所有参数与论文中完全一致。您可以在下面看到论文中的图表和上述代码的输出。结果与论文完全一致，说明EM算法的实现是正确的。也可以使用数值求解器来最大化完整对数似然函数的条件期望，但在这种情况下使用解析解更容易。现在让我们试着让事情变得更复杂一点。假设选择每个硬币的概率是未知的。然后我们将混合两个二项式，其中选择第一个硬币的概率为tau，选择第二个硬币的概率为1-tau。我们重复与之前相同的步骤。定义全似然函数：在前面的例子中，Z给定theta的条件概率是常数，现在它是带有参数tau的伯努利分布的PMF。计算完整的对数似然函数：找到给定X和theta的隐藏变量Z的条件分布：计算对数似然的条件期望：剩下的就是最大化关于参数theta的条件期望。上式中涉及概率p_j的各项与前面完全相同，所以p_1和p_2的解是一样的。最大化tau的条件期望可以得到：用Python实现这个算法的代码如下#模型参数sp_1=0.1p_2=0.8tau_1=0.3tau_2=1-tau_1m=10#每个样本的翻转次数=10#number样本生成数据np.random.seed(123)dists=[stats.binom(n=m,p=p_1),stats.binom(n=m,p=p_2)]xs=[dists[x].rvs()forxinnp.random.choice([0,1],size=n,p=[tau_1,tau_2])]#randominitialguessnp.random.seed(123)theta=[np.random.rand()for_inrange(3)]last_ll=0max_iter=100foriinrange(max_iter):tau,p_1,p_2=thetaT_1s=[]T_2s=[]#E-steplls=[]forxinxs:denom=(tau*stats.binom(n=m,p=p_1).pmf(x)+(1-tau)*stats.binom(n=m,p=p_2).pmf(x))T_1=tau*统计数据。binom(n=m,p=p_1).pmf(x)/denomT_2=(1-tau)*stats.binom(n=m,p=p_2).pmf(x)/denomT_1s.append(T_1)T_2s.append(T_2)lls.append(T_1*np.log(tau*stats.binom(n=m,p=p_1).pmf(x))+T_2*np.log(tau*stats.binom(n=m,p=p_2).pmf(x)))#M步T_1s=np.array(T_1s)T_2s=np.array(T_2s)tau=np.sum(T_1s)/np_1=np.dot(T_1s,xs)/(m*np.sum(T_1s))p_2=np.dot(T_2s,xs)/(m*np.sum(T_2s))print(f'step:{i},tau={tau},p1={p_1:.2f},p2={p_2:.2f},log_likelihood={sum(lls):.2f}')theta=[tau,p_1,p_2]#如果abs(sum(lls)-last_ll)<0.001:breakelse:last_ll=sum(lls)在前面的例子中，算法运行了10次迭代，复制了论文中的结果，但实际上我们运行算法直到收敛，这意味着当Stop当数字似然停止改进时的算法。算法的输出如上所示。由于问题的对称性，交换了p_1和p_2的概率，但结果仍然正确且接近参数的真实值：我们有p=0.85概率为0.8，p=0.1概率为0.2（真实值：p=0.8，概率为0.7，p=0.1，概率为0.3）。本文的目的是通过一些简单的例子来演示EM算法的基础知识。本文的完整代码可以在这里找到：https://avoid.overfit.cn/post/f618e5e3c5304fceb36abbdec8816107作者：AlexanderPavlov
