最大似然分类python 最大似然分类算法

文章目录

• 参考资料
• 1. 最大似然估计

• 1.1 原理
• 1.2 示例

• 2. EM算法

• 2.1 原理
• 2.2 示例

参考资料

• 统计计算中的优化问题

1. 最大似然估计

1.1 原理

统计中许多问题的计算最终都归结为一个最优化问题， 典型代表是最大似然估计(MLE)、各种拟似然估计方法、 非线性回归、惩罚函数方法（如svm、lasso）等。

最大似然估计经常需要用最优化算法计算， 最大似然估计问题有自身的特点， 可以直接用一般优化方法进行最大似然估计的计算， 但是利用最大似然估计的特点可以得到更有效的算法。

设总体 $\boldsymbol{X}$ 有概率密度（连续型随机变量）或概率分布（离散型随机变量） $p(\boldsymbol{x} \mid \boldsymbol{\theta}), \boldsymbol{\theta}$ 为 $m$ 维的分布参数。有了一组样本 $\boldsymbol{X}_1, \boldsymbol{X}_2, \ldots, \boldsymbol{X}_n$ 后，似然函数 为
$L(\boldsymbol{\theta})=\prod_{i=1}^n p\left(\boldsymbol{X}_i \mid \boldsymbol{\theta}\right)$
对数似然函数为
$l(\boldsymbol{\theta})=\sum_{i=1}^n \ln p\left(\boldsymbol{X}_i \mid \boldsymbol{\theta}\right)$

最大化的步骤通过对$l(\boldsymbol{\theta})$求导等于0来解得。

1.2 示例

• 例1：设 $X \sim b(1, p), X_1, \ldots, X_n$ 是来自 $X$ 的一个样本， $x_1, \ldots, x_n$ 为观察值。试求参数 $p$ 的最大似然 估计。
  解：可知 $X$ 的分布律为:
  $P\{X=x\}=p^x(1-p)^{1-x}$
  故似然函数为：
  $\begin{aligned} &L(p \mid \boldsymbol{x})=\prod_{i=1}^n p^{x_i}(1-p)^{1-x_i}=p^{\sum_{i=1}^n x_i}(1-p)^{n-\sum_{i=1}^n x_i} \\ &\ln L(p \mid \boldsymbol{x})=\sum_{i=1}^n x_i \ln p+\left(n-\sum_{i=1}^n x_i\right) \ln (1-\mathrm{p}) \end{aligned}$
  对待估参数求导为 0 , 有：
  $\frac{d}{d p} \ln L(p \mid \boldsymbol{x})=\frac{\sum_{i=1}^n x_i}{p}-\frac{n-\sum_{i=1}^n x_i}{1-p}=0$
  解得 $p$ 的最大似然估计值为 $\hat{p}=\bar{x}$ ，最大估计量为 $\hat{p}=\bar{X}$.
• 例2： 设 $X \sim N\left(\mu, \sigma^2\right), \mu, \sigma^2$ 为未知参数, $x_1, \ldots, x_n$ 为来自 $X$ 的一组观察值, 求 $\mu, \sigma^2$ 的最大似然估计量。
  解：$X$ 的概率密度为
  $f\left(x ; \mu, \sigma^2\right)=\frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{(x-\mu)^2}{2 \sigma^2}}$
  似然函数为:
  $\begin{gathered} L(\boldsymbol{\mu, \sigma} \mid \boldsymbol{x})=\prod_{i=1}^n \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(x_i-\mu\right)^2}{2 \sigma^2}} \\ \ln L(\boldsymbol{\mu, \sigma} \mid \boldsymbol{x})=-\operatorname{nln} \sqrt{2 \pi} \sigma-\frac{1}{2 \sigma^2} \sum_{i=1}^n\left(x_i-\mu\right)^2 \\ \left\{\begin{array}{c} \frac{d}{d \mu} \ln L(\boldsymbol{\mu, \sigma} \mid \boldsymbol{x})=\frac{1}{2 \sigma^2} 2 \sum_{i=1}^n\left(x_i-\mu\right)=0 \\ \frac{d}{d \sigma} \ln L(\boldsymbol{\mu, \sigma} \mid \boldsymbol{x})=-\frac{n}{\sigma}+\frac{1}{\sigma^3} 2 \sum_{i=1}^n\left(x_i-\mu\right)^2=0 \end{array}\right. \end{gathered}$
  最后解得：
  $\left\{\begin{array}{c} \mu=\bar{x} \\ \sigma^2=\frac{1}{n} \sum_{i=1}^n\left(x_i-\mu\right)^2 \end{array}\right.$
  因此得到 $\mu, \sigma^2$ 的最大似然估计量分别为
  $\hat{\mu}=\bar{X}, \hat{\sigma}^2=\frac{1}{n} \sum_{i=1}^n\left(X_i-\bar{X}\right)^2$
• 例3 (多项分布)： 设 $X$ 取值于 $\{1,2,3,4\}$, 分布概率为 $p(x \mid \theta) \triangleq \pi_x(\theta):$
  $\begin{aligned} \pi_1(\theta) &=\frac{1}{4}(2+\theta), & \pi_2(\theta) &=\frac{1}{4}(1-\theta) \\ \pi_3(\theta) &=\frac{1}{4}(1-\theta), & \pi_4(\theta) &=\frac{1}{4} \theta \end{aligned}$
  设 $n$ 次试验得到的 $X$ 值有 $n_j$ 个 $j(j=1,2,3,4)$ ，求参数 $\theta$
  解：观测数据的对数似然函数为(去掉了与参数无关的加性常数)
  $l(\theta)=n_1 \ln (2+\theta)+\left(n_2+n_3\right) \ln (1-\theta)+n_4 \ln \theta .$
  令 $l^{\prime}(\theta)=0$ 得到一个关于 $\theta$ 的二次方程，由此写出 $\operatorname{argmax} l(\theta)$ 的解析表达式:
  $l^{\prime}(\theta)=\frac{n_1}{2+\theta}-\frac{n_2+n_3}{1-\theta}+\frac{n_4}{\theta}$
  令 $l^{\prime}(\theta)=0$ ，即求解二次方程
  $n \theta^2+\left(-n_1+2 n_2+2 n_3+n_4\right) \theta-2 n_4=0 .$
  得最大似然估计为
  $\hat{\theta}=\frac{-b+\sqrt{b^2+8 n n_4}}{2 n}$
  其中 $b=-n_1+2 n_2+2 n_3+n_4$

2. EM算法

2.1 原理

EM算法最初用于缺失数据模型参数估计，现在已经用在许多优化问题中。设模型中包含 $\boldsymbol{X}_{\text {obs }} 和 \boldsymbol{X}_{\mathrm{mis}}$ 两个 随机成分， 有联合密度函数或概率函数 $f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right), \boldsymbol{\theta}$ 为未知参数。称 $f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right)$ 为完全数据的 密度，一般具有简单的形式。实际上我们只有 $\boldsymbol{X}_{\mathrm{obs}}$ 的观测数据 $\boldsymbol{X}_{\mathrm{obs}}=\boldsymbol{x}_{\mathrm{obs}} ， \boldsymbol{X}_{\mathrm{mis}}$ 不能观测得到， 这一 部分可能是缺失观测数据，也可能是潜在影响因素。所以实际的似然函数为
$\tag{1} L(\boldsymbol{\theta})=f\left(\boldsymbol{x}_{\mathrm{obs}} \mid \boldsymbol{\theta}\right)=\int f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right) d \boldsymbol{x}_{\mathrm{mis}},$
这个似然函数通常比完全数据的似然函数复杂得多，所以很难直接从 $L(\boldsymbol{\theta})$

EM算法的想法是，已经有了参数的近似估计值 $\boldsymbol{\theta}^{(t)}$ 后， 假设 $\left(\boldsymbol{X}_{\mathrm{obs}}, \boldsymbol{X}_{\mathrm{mis}}\right)$ 近似服从完全密度 $f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}^{(t)}\right)$, 这里 $\boldsymbol{X}_{\mathrm{obs}}=\boldsymbol{x}_{\mathrm{obs}}$ 已知，所以认为 $\boldsymbol{X}_{\mathrm{mis}}$ 近似服从由 $f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}^{(t)}\right)$ 导出的条件分 布
$\tag{2} f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right)=\frac{f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}^{(t)}\right)}{f\left(\boldsymbol{x}_{\mathrm{obs}} \mid \boldsymbol{\theta}^{(t)}\right)}$
其中 $f\left(\boldsymbol{x}_{\mathrm{obs}} \mid \boldsymbol{\theta}^{(t)}\right)$ 是由 $f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}^{(t)}\right)$ 决定的边缘密度。据此近似条件分布，在完全数据对数似然函数 $\ln f\left(\boldsymbol{X}_{\mathrm{obs}}, \boldsymbol{X}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right)$ 中， 把 $\boldsymbol{X}_{\mathrm{obs}}=\boldsymbol{x}_{\mathrm{obs}}$ 看成已知， 关于未知部分 $\boldsymbol{X}_{\mathrm{mis}}$ 按密度 $f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\boldsymbol { x }}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right)$ 求期 望，得到 $\boldsymbol{\theta}$ 的函数 $Q_t(\boldsymbol{\theta})$ ，再求 $Q_t(\boldsymbol{\theta})$ 的最大值点作为下一个 $\boldsymbol{\theta}^{(t+1)}$

EM算法每次迭代有如下的E步（期望步）和M步（最大化步）：

• E步: 计算完全数据对数似然函数的期望 $Q_t(\boldsymbol{\theta})=E\left\{\ln f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{X}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right)\right\}$, 其中期望针对随机变量 $\boldsymbol{X}_{\mathrm{mis}}$ ， 求期望时假定 $\boldsymbol{X}_{\mathrm{mis}}$ 服从条件密度 $f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right)$
• M步: 求 $Q_t(\boldsymbol{\theta})$ 的最大值点，记为 $\boldsymbol{\theta}^{(t+1)}$

定理1： EM算法得到的估计序列 $\boldsymbol{\theta}^{(t)}$ 使得公式(1)中的似然函数值 $L\left(\boldsymbol{\theta}^{(t)}\right)$ 单调不减。
证明: 对任意参数 $\boldsymbol{\theta}$ ，有
$\begin{aligned} \ln L(\theta)=& \ln f\left(\boldsymbol{x}_{\mathrm{obs}} \mid \theta\right) \cdot \int f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \\ =&\ln\frac{f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right)}{f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}\right)} \int f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \\ =& \int\left[\ln f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right)-\ln f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}\right)\right] f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \\ =& \int \ln f\left(\boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{\theta}\right) f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \\ &-\int \ln f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}\right) f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \\ =& Q_t(\boldsymbol{\theta})-\int \ln f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}\right) f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \end{aligned}$
由信息不等式知
$\begin{aligned} & \int \ln f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}\right) f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \\ \leq & \int \ln f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \end{aligned}$
又EM迭代使得 $Q_t\left(\boldsymbol{\theta}^{(t+1)}\right) \geq Q_t\left(\boldsymbol{\theta}^{(t)}\right)$, 所以
$\begin{aligned} \ln L\left(\boldsymbol{\theta}^{(t+1)}\right) & \geq Q_t\left(\boldsymbol{\theta}^{(t)}\right)-\int \ln f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) f\left(\boldsymbol{x}_{\mathrm{mis}} \mid \boldsymbol{x}_{\mathrm{obs}}, \boldsymbol{\theta}^{(t)}\right) d \boldsymbol{x}_{\mathrm{mis}} \\ &=\ln L\left(\boldsymbol{\theta}^{(t)}\right) . \end{aligned}$

定理证毕。

在适当正则性条件下， EM算法的迭代序列依概率收敛到$\theta^{(t)}$的最大值点。 但是， 定理（1）仅保证EM算法最终能收敛， 但不能保证EM算法会收敛到似然函数的全局最大值点， 算法也可能收敛到局部极大值点或者鞍点。

在实际问题中， 往往E步和M步都比较简单， 有时E步和M步都有解析表达式， 这时EM算法实现很简单。 EM算法优点是计算稳定， 可以保持原有的参数约束， 缺点是收敛可能很慢， 尤其是接近最大值点时可能收敛更慢。 如果公式（1）中的似然函数不是凸函数， 算法可能收敛不到全局最大值点， 遇到这样的问题可以多取不同初值比较， 用矩估计等合适的近似值作为初值。

原理暂时看不懂没关系，结合后面例题看就更容易懂了。

2.2 示例

(混合分布) EM算法可以用来估计混合分布的参数。 设随机变量 $Y_1 \sim \mathrm{N}\left(\mu_1, \delta_1\right), Y_2 \sim \mathrm{N}\left(\mu_2, \delta_2\right)$, $Y_1, Y_2$ 独立。记 $N(\mu, \delta)$ 的密度为 $f(x \mid \mu, \delta)$ 。设随机变量 $W \sim \mathrm{b}(1, \lambda), 0<\lambda<1, W$ 与 $Y_1, Y_2$ 独立，令
$X=(1-W) Y_1+W Y_2,$
则 $W=0$ 条件下 $X \sim \mathrm{N}\left(\mu_1, \delta_1\right), W=1$ 条件下 $X \sim \mathrm{N}\left(\mu_2, \delta_2\right)$, 但 $X$ 的边缘密度为
$f(x \mid \boldsymbol{\theta})=(1-\lambda) f\left(x \mid \mu_1, \delta_1\right)+\lambda f\left(x \mid \mu_2, \delta_2\right),$
其中 $\boldsymbol{\theta}=\left(\mu_1, \delta_1, \mu_2, \delta_2, \lambda\right)$

设 $X$ 有样本 $\boldsymbol{X}=\left(X_1, \ldots, X_n\right)$, 样本值为 $\boldsymbol{x}$ ， 实际观测数据的似然函数为
$L(\boldsymbol{\theta})=\prod_{i=1}^n f\left(x_i \mid \boldsymbol{\theta}\right)$
这个函数是光滑函数但是形状很复杂， 直接求极值很容易停留在局部极值点。
用EM算法，以 $\boldsymbol{W}=\left(W_1, \ldots, W_n\right)$ 为没有观测到的部分, 完全数据的似然函数和对数似然函数为
$\begin{aligned} \tilde{L}(\boldsymbol{\theta} \mid \boldsymbol{x}, \boldsymbol{W})=& \prod_{W_i=0} f\left(x_i \mid \mu_1, \delta_1\right) \prod_{W_i=1} f\left(x_i \mid \mu_2, \delta_2\right) \lambda^{\sum_{i=1}^n W_i}(1-\lambda)^{n-\sum_{i=1}^n W_i} \\ \tilde{l}(\boldsymbol{\theta} \mid \boldsymbol{x}, \boldsymbol{W})=& \sum_{i=1}^n\left[\left(1-W_i\right) \log f\left(x_i \mid \mu_1, \delta_1\right)+W_i \log f\left(x_i \mid \mu_2, \delta_2\right)\right] \\ &+\left(\sum_{i=1}^n W_i\right) \log \lambda+\left(n-\sum_{i=1}^n W_i\right) \log (1-\lambda) \end{aligned}$
在E步，设已有 $\boldsymbol{\theta}$ 的近似值 $\boldsymbol{\theta}^{(t)}=\left(\mu_1^{(t)}, \delta_1^{(t)}, \mu_2^{(t)}, \delta_2^{(t)}, \lambda^{(t)}\right)$, 以 $\boldsymbol{\theta}^{(t)}$ 为分布参数，在 $\boldsymbol{X}=\boldsymbol{x}$ 条件下， $\boldsymbol{W}_i$ 的 条件分布为
$\begin{aligned} \gamma_i^{(t)} & \triangleq P\left(W_i=1 \mid \boldsymbol{x}, \boldsymbol{\theta}^{(t)}\right)=P\left(W_i=1 \mid X_i=x_i, \boldsymbol{\theta}^{(t)}\right) \\ &=\frac{\lambda^{(t)} f\left(x_i \mid \mu_2^{(t)}, \delta_2^{(t)}\right)}{\left(1-\lambda^{(t)}\right) f\left(x_i \mid \mu_1^{(t)}, \delta_1^{(t)}\right)+\lambda^{(t)} f\left(x_i \mid \mu_2^{(t)}, \delta_2^{(t)}\right)} . \end{aligned}$
这里的推导类似于逆概率公式。利用 $W_i$ 的条件分布求完全数据对数似然的期望，得
$\tag{3} \begin{aligned} Q_t(\boldsymbol{\theta})=& \sum_{i=1}^n\left[\left(1-\gamma_i^{(t)}\right) \log f\left(x_i \mid \mu_1, \delta_1\right)+\gamma_i^{(t)} \log f\left(x_i \mid \mu_2, \delta_2\right)\right] \\ &+\left(\sum_{i=1}^n \gamma_i^{(t)}\right) \log \lambda+\left(n-\sum_{i=1}^n \gamma_i^{(t)}\right) \log (1-\lambda) . \end{aligned}$
令 $\nabla Q_t(\boldsymbol{\theta})=\mathbf{0}$ ，求得 $Q_t(\boldsymbol{\theta})$ 的最大值点 $\boldsymbol{\theta}^{(t+1)}$ 为
$\tag{4} \left\{\begin{array}{l} \mu_1^{(t+1)}=\frac{\sum_{i=1}^n\left(1-\gamma_i^{(t)}\right) x_i}{\sum_{i=1}^n\left(1-\gamma_i^{(t)}\right)} \\ \delta_1^{(t+1)}=\frac{\sum_{i=1}^n\left(1-\gamma_i^{(t)}\right)\left(x_i-\mu_1^{(t+1)}\right)^2}{\sum_{i=1}^n\left(1-\gamma_i^{(t)}\right)} \\ \mu_2^{(t+1)}=\frac{\sum_{i=1}^n \gamma_i^{(t)} x_i}{\sum_{i=1}^n \gamma_i^{(t)}} \\ \delta_2^{(t+1)}=\frac{\sum_{i=1}^n \gamma_i^{(t)}\left(x_i-\mu_2^{(t+1)}\right)^2}{\sum_{i=1}^n \gamma_i^{(t)}} \\ \lambda^{(t+1)}=\frac{1}{n} \sum_{i=1}^n \gamma_i^{(t)} \end{array}\right.$
适当选取初值 $\boldsymbol{\theta}^{(0)}$ 用公式(3)和(4)迭代就可以计算 $\boldsymbol{\theta}$

• (多项分布)： 设 $X$ 取值于 $\{1,2,3,4\}$, 分布概率为 $p(x \mid \theta) \triangleq \pi_x(\theta):$
  $\begin{aligned} \pi_1(\theta) &=\frac{1}{4}(2+\theta), & \pi_2(\theta) &=\frac{1}{4}(1-\theta) \\ \pi_3(\theta) &=\frac{1}{4}(1-\theta), & \pi_4(\theta) &=\frac{1}{4} \theta \end{aligned}$
  设 $n$ 次试验得到的 $X$ 值有 $n_j$ 个 $j(j=1,2,3,4)$ ，求参数 $\theta$
  $P(X=12)=\frac{1}{4} \theta_{\text {。 }}$ 这时 $Z_2+n_4$ 代表结果12和结果4的出现次数， 这两种结果出现概率为 $\frac{1}{2} \theta$ ，其它结 果 $(11,2,3)$ 的出现概率为 $1-\frac{1}{2} \theta_{\text {。 }}$ 令 $Y=Z_2+n_4$ ，则 $Y \sim \mathrm{B}\left(n, \frac{1}{2} \theta\right)$ 。
  数据 $\left(Z_1, Z_2, n_2, n_3, n_4\right)$ 的全似然函数为
  $L_c(\theta) \propto\left(\frac{1}{2}\right)^{Z_1}\left(\frac{\theta}{4}\right)^{Z_2}\left(\frac{1}{4}-\frac{\theta}{4}\right)^{n_2+n_3}\left(\frac{\theta}{4}\right)^{n_4}$
  对数似然函数 (差一个与 $\theta$ 无关的常数项) 为
  $l_c(\theta)=\ln L_c(\theta)=\left(Z_2+n_4\right) \ln \theta+\left(n_2+n_3\right) \ln (1-\theta)$
  在EM迭代中， 假设已经得到的参数 $\theta$ 近似值为 $\theta^{(t)}$ ， 设 $\theta=\theta^{(t)}$ ， 在给定 $n, n_1, n_2, n_3, n_4$ 条件下求 $l_c(\theta)$ 的 条件期望， 这时 $Z_2$ 的条件分布为
  $\mathrm{B}\left(n_1, \frac{\theta^{(t)}}{2+\theta^{(t)}}\right),$
  于是
  $Z_2^{(t)}=E\left(Z_2 \mid \theta^{(t)}, n, n_1, n_2, n_3, n_4\right)=n_1 \frac{\theta^{(t)}}{2+\theta^{(t)}},$
  从而完全对数似然函数的条件期望为
  $Q_t(\theta)=E\left(\ln L_c(\theta) \mid \theta^{(t)}, n, n_1, n_2, n_3, n_4\right)=\left(Z_2^{(t)}+n_4\right) \ln \theta+\left(n_2+n_3\right) \ln (1-\theta) .$
  求解 $Q_t(\theta)$ 的最大值，令
  $\frac{d}{d \theta} Q_t(\theta)=\frac{n_4+Z_2^{(t)}}{\theta}-\frac{n_2+n_3}{1-\theta}=0$
  得下一个参数近似值为
  $\theta^{(t+1)}=\frac{n_4+Z_2^{(t)}}{n_2+n_3+n_4+Z_2^{(t)}} .$
  于是, EM迭代步骤从某个 $\theta^{(0)}$ 出发，比如 $\theta^{(0)}=\frac{1}{2}$ ， 在第 $t$ 步计算
  $Z_2^{(t)}=n_1 \frac{\theta^{(t)}}{2+\theta^{(t)}}, \quad \theta^{(t+1)}=\frac{n_4+Z_2^{(t)}}{n_2+n_3+n_4+Z_2^{(t)}}$
  迭代到两次的近似参数值变化小于 $\epsilon=10^{-6}$