冀望最大化（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 # number of flips in experimentn = 5 # number of experimentsp_1 = 0.8p_2 = 0.3xs = [] # number of heads in each experimentzs = [0,0,1,0,1] # which coin to flipnp.random.seed(5)for i in zs: if i==0: xs.append(stats.binom(n=m, p=p_1).rvs()) # flip coin 1 elif i==1: xs.append(stats.binom(n=m, p=p_2).rvs()) # flip coin 2 xs = np.array(xs)print(xs)exp1 = [0,1,3] # indexes of experiments with coin 1exp2 = [2,4] # indexes of experiments with coin 2print('Analytical solutions:')print('p1: ', xs[exp1].sum() / (m * len(exp1)))print('p2: ', xs[exp2].sum() / (m * len(exp2)))这些实现都比较简单。咱们只是实现下面计算的公式并插入数字。运行此代码的后果如下所示。
解决这个问题的另一种办法是应用数值求解器。咱们须要找到一个最大化对数似然函数的解决方案，当应用数值求解器时，不须要计算导数并手动求解最大化对数似然函数的参数。只需实现一个咱们想要最大化的函数并将其传递给数值求解器。
因为 Python 中的大多数求解器旨在最小化给定函数，因而咱们实现了一个计算负对数似然函数的函数（因为最小化负对数似然函数与最大化对数似然函数雷同）。
代码和后果如下所示。
ef neg_log_likelihood(probs, m, xs, zs): ''' compute negative log-likelihood ''' ll = 0 for x,z in zip(xs,zs): ll += stats.binom(p=probs[z], n=m).logpmf(x) return -ll res = optimize.minimize(neg_log_likelihood, [0.5,0.5], bounds=[(0,1),(0,1)], args=(m,xs,zs), method='tnc')print('Numerical solution:')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）：计算残缺对数似然函数绝对于 Z 给定数据 X 的以后条件散布和以后参数估计 theta 的条件期望
- 最大化步骤（M-step）：找到最大化该冀望的参数 theta 的值
能够应用贝叶斯定理在给定 X_i 和 theta 的状况下找到 Z_i 的条件散布：
当初定义齐全对数似然的条件期望如下：
插入残缺的对数似然函数并重新排列：
这样就实现了 E-step。在 M-step中，咱们依据参数 theta 最大化下面计算的函数：
在这个的示例中，能够通过剖析找到参数（通过求导并将其设置为零）。这里就不演示整个过程，只是提供答案。
在上面的实现中将应用与论文中雷同数据来查看是否取得了雷同的后果。Python 代码如下
m = 10 # number of flips in each samplen = 5 # number of samplesxs = np.array([5,9,8,4,7])theta = [0.6, 0.5] # initial guess p_1, p_2for i in range(10): p_1,p_2 = theta T_1s = [] T_2s = [] # E-step for x in xs: 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-step 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。
咱们反复与之前雷同的步骤。定义齐全似然函数：
在后面的例子中，给定 theta 的 Z 的条件概率是恒定的，当初它是带有参数 tau 的伯努利散布的 PMF。
计算残缺的对数似然函数：
求给定 X 和 theta 的暗藏变量 Z 的条件散布：
计算对数似然的条件期望：
剩下的就是最大化对于参数 theta 的条件期望。上式中波及概率 p_j 的项与之前的齐全一样，所以 p_1 和 p_2 的解是雷同的。最大化对于 tau 的条件期望，能够失去：
Python 中实现这个算法的代码如下
# model parametersp_1 = 0.1p_2 = 0.8tau_1 = 0.3tau_2 = 1-tau_1m = 10 # number of flips in each samplen = 10 # number of samples# generate datanp.random.seed(123)dists = [stats.binom(n=m, p=p_1), stats.binom(n=m, p=p_2)]xs = [dists[x].rvs() for x in np.random.choice([0,1], size=n, p=[tau_1,tau_2])]# random initial guessnp.random.seed(123)theta = [np.random.rand() for _ in range(3)]last_ll = 0max_iter = 100for i in range(max_iter): tau,p_1,p_2 = theta T_1s = [] T_2s = [] # E-step lls = [] for x in xs: 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 * stats.binom(n=m,p=p_1).pmf(x) / denom T_2 = (1-tau) * stats.binom(n=m,p=p_2).pmf(x) / denom T_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-step T_1s = np.array(T_1s) T_2s = np.array(T_2s) tau = np.sum(T_1s) / n 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}, tau={tau}, p1={p_1:.2f}, p2={p_2:.2f}, log_likelihood={sum(lls):.2f}') theta = [tau, p_1, p_2] # stop when likelihood doesn't improve anymore if abs(sum(lls) - last_ll) < 0.001: break else: last_ll=sum(lls)在后面的示例中，将算法运行了 10 次迭代，复制论文中的后果。然而在实践中，咱们运行算法直到收敛，这意味着当对数似然进行改良时就进行算法。
算法的输入如上所示。因为问题的对称性，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
作者：Alexander Pavlov