椒花蛾(peppered moth)的色彩已经被确认是通过某单个基因决定的, 该基因有三个可能的等位基因: C, I和T. 其中C对I是显性的, T对I是隐性的. 因此基因型CC, CI和CT导致(深色)黑色的表型. 基因型TT导致浅色的表型. II和IT导致一种中间的表型, 外观上变化很广泛, 但通常以中间色彩杂色而成. 因此有6种
基因型, 3种表现型.
在英国和北美, 在燃煤的工业区, 浅色的蛾子几乎以及被深色的蛾子所替代. 这种等位基因频率在种群内的变化被认为是在人类时间刻度下可以观察到生物微进化的一个例子. (被实验证实的)理论结果是"鸟类在不同反射背景下对蛾子的捕食程度明显不同". 在燃煤的工业区, 污染减弱了黑色表型蛾子栖息在树皮表面的反射程度, 因而对其有利. 当环境改善后, 浅色表型增加, 黑色表型减少, 这并不令人奇怪.因此, 监视这种蛾子的等位基因C,I,T的变化, 除了可以研究微进化过程外, 也可以用于环境污染监控.
假设Hardy-Weinberg平衡律成立, 记各个等位基因在种群中的频率为 \(p_{C}, p_I, p_T, (p_{C}+p_I+p_T=1)\), 则基因型CC, CI, CT, II, IT, TT的
频率分别为 \(p^2_{C}\), \(2p_{C} p_{I}\), \(2p_{C} p_{T}\), \(p_{I}^2\), \(2p_{I}p_{T}\), \(p_{T}^2\).
假设我们随机捕了 \(n\) 只蛾子, 其中黑色表型的有 \(n_C=85\) 只, 中间表型的有 \(n_I=196\) 只, 浅色表型的有 \(n_T=341\) 只. 求等位基因的概率 \(p_{C},p_{I},p_{T}\)的最大似然估计.
由于各个基因型的频数不可观测: $$ \begin{array}{|c|ccc|cc|c|} \hline \text{基因型}& \mathrm{CC} & \mathrm{CI} & \mathrm{CT} & \text { II } & \text { IT } & \text { TT } \\ \hline \text { 不可观测 } & n_{C C} & n_{C I} & n_{C T} & n_{I I} & n_{I T} & n_{T T} \\ \hline\text { 观测 } & & n_C & & & n_I & n_T \\ & &\text { 黑色 } & & & \text { 中间 } & \text { 浅色 } \\ \hline \end{array} $$ 在左侧输入观测数和初始概率, EM迭代过程结果如下:
- 迭代次数 t=0,1,...,10
- 第t步迭代等位基因C的概率 \(pct=p_C^{(t)}\)
- 第t步迭代等位基因I的概率 \(pit=p_I^{(t)}\)
- 第t步迭代相对变化量 \(rcc=\|p^{(t)}-p^{(t-1)}\|/\|p^{(t-1)}\|, p=(p_C, p_I)\), 可用于收敛标准
- 第t步迭代相对t-1步迭代变化量 \(d1=(p_C^{(t)}-p_C^{(10)})/(p_C^{(t-1)}-p_C^{(10)})\), 验证收敛速度为线性
- 第t步迭代相对t-1步迭代变化量 \(d2=(p_I^{(t)}-p_I^{(10)})/(p_I^{(t-1)}-p_I^{(10)})\)为平稳点, 验证收敛速度为线性
此时观测数据为 \(x=(n_{ C},n_{ I},n_{ T})\), 而完全数据为 \(y=(n_{ CC},n_{ CI},n_{ CT}\), \(n_{ II},n_{ IT},n_{ TT})\),
这里 \(n_{ CC}+n_{ CI}+n_{ CT}=n_{ C}\), \(n_{ II}+n_{ IT}=n_{ I}\), 以及 \(n_{ TT}=n_{ T}\). 我们希望由此数据
来估计各个等位基因的概率 \(p_{ C},p_{ I},p_{ T}\).
从而完全数据的似然函数为
$$
ln f(y|p)=n_{ CC}\ln p_{ C}^2+n_{ CI}\ln (2p_{ C}p_{ I})+n_{ CT}\ln (2p_{ C}p_{ T})\\
+n_{ II}\ln (p_{ I}^2)+n_{ IT}\ln (2p_{ I}p_{ T})+n_{ TT}\ln (p_{ TT}^2)\\
+\ln \frac{n!}{n_{ CC}!n_{ CI}!n_{ CT}!n_{ II}!n_{ IT}!n_{ TT}!}
$$
由于完全数据不可观测, 若记 \(N_{ CC},N_{ CI},N_{ CT}, N_{ II},N_{ IT},N_{ TT}\)分别为各个基因
型的个数, \(p=(p_{ CC},p_{ CI},p_{ CT},p_{ II},p_{ IT},p_{ TT})\), 则有
$$
\begin{aligned}
1.& N_{ CC},N_{ CI},N_{ CT}|n_{ C},n_{ I},n_{ T},p \sim
MN\left(n_{ C},\frac{(p_{ C})^2}{1-(p_{ I}+p_{ T})^2},
\frac{2p_{ C}p_{ I}}{1-(p_{ I}+p_{ T})^2},\frac{2p_{ C}p_{ T}}{1-(p_{ I}+p_{ T})^2}\right)\\
2.& N_{ II},N_{ IT}|n_{ C},n_{ I},n_{ T},p \sim
MN\left(n_{ I},\frac{(p_{ I})^2}{2p_{ C}p_{ I}+(p_{ T})^2},
\frac{2p_{ C}p_{ I}}{2p_{ C}p_{ I}+(p_{ T})^2}\right)\\
3. & N_{ TT}|n_{ C},n_{ I},n_{ T},p \sim B(n_{ T},p_{ T}^2)
\end{aligned}
$$
从而在EM算法的第k步中, 对完全似然函数取期望得到
$$
Q(p|p^{(k)})=n_{ CC}^{(k)}\ln (p_{ C}^2)+n_{ CI}^{(k)}\ln (\ln (2p_{ C}p_{ I})+n_{ CT}^{(k)}\ln (2p_{ C}p_{ T})\\
+n_{ II}^{(k)}\ln (p_{ I}^2)+n_{ IT}^{(k)}\ln (2p_{ I}p_{ T})+n_{ TT}^{(k)}\ln (p_{ TT}^2)\\
+k(n_{ C},n_{ I},n_{ T},p^{(k)}),
$$
其中
$$
n_{ CC}^{(k)}=E\{N_{ CC}|n_{ C},n_{ I},n_{ T},p^{(k)}\}=n_{ C}\frac{(p_{ C}^{(k)})^2}{1-(p_{ I}^{(k)}+p_{ T}^{(k)})^2},\\
n_{ CI}^{(k)}=E\{N_{ CC}|n_{ C},n_{ I},n_{ T},p^{(k)}\}=n_{ C}\frac{2p_{ C}^{(k)}p_{ I}^{(k)}}{1-(p_{ I}^{(k)}+p_{ T}^{(k)})^2},\\
n_{ CT}^{(k)}=E\{N_{ CC}|n_{ C},n_{ I},n_{ T},p^{(k)}\}=n_{ C}\frac{2p_{ C}^{(k)}p_{ T}^{(k)}}{1-(p_{ I}^{(k)}+p_{ T}^{(k)})^2},\\
n_{ II}^{(k)}=E\{N_{ II}|n_{ C},n_{ I},n_{ T},p^{(k)}\}=n_{ I}\frac{(p_{ I}^{(k)})^2}{2p^{(k)}_Cp^{(k)}_I+(p^{(k)}_T)^2}\\
n_{ IT}^{(k)}=E\{N_{ IT}|n_{ C},n_{ I},n_{ T},p^{(k)}\}=n_{ I}\frac{2p_{ I}^{(k)}p_{ T}^{(k)}}{2p^{(k)}_Cp^{(k)}_I+(p^{(k)}_T)^2}.
$$
最后一项
$$
k(n_{ C},n_{ I},n_{ T},p^{(k)})
=E\Big\{\ln \frac{n!}{n_{ CC}!n_{ CI}!n_{ CT}!n_{ II}!n_{ IT}!n_{ TT}!}|n_{ C},n_{ I},n_{ T},p^{(k)}\Big\}
$$
与 \(p\)无关.
下面对 \(Q(p|p^{(k)})\)进行最大化, 注意 \(p_{ C}+p_{ I}+p_{ T}=1\), 于是关于 \(p_{ C},p_{ I}\)求导
$$
\frac{\partial Q(p|p^{(k)})}{\partial p_{ C}}=\frac{2n_{ CC}^{(k)}+n_{ CI}^{(k)}+n_{ CT}^{(k)}}{p_{ C}}-
\frac{2n_{ TT}^{(k)}+n_{ CT}^{(k)}+n_{ IT}^{(k)}}{1-p_{ C}-p_{ I}},\\
\frac{\partial Q(p|p^{(k)})}{\partial p_{ I}}=\frac{2n_{ II}^{(k)}+n_{ II}^{(k)}+n_{ CI}^{(k)}}{p_{ I}}-
\frac{2n_{ TT}^{(k)}+n_{ CT}^{(k)}+n_{ IT}^{(k)}}{1-p_{ C}-p_{ I}}
$$
令导数为零,得到
$$
p_{ C}^{(k+1)}=\frac{2n_{ CC}^{(k)}+n_{ CT}^{(k)}+n_{ CI}^{(k)}}{2n},\\
p_{ I}^{(k+1)}=\frac{2n_{ II}^{(k)}+n_{ IT}^{(k)}+n_{ CI}^{(k)}}{2n},\\
p_{ T}^{(k+1)}=\frac{2n_{ TT}^{(k)}+n_{ CT}^{(k)}+n_{ IT}^{(k)}}{2n}.
$$
代入数据计算得到 \(\hat{p}_{ C}=0.071, \hat{p}_{ I}=0.189, \hat{p}_{ T}=0.74\).