双序列比对的理论基础之建造替换矩阵的合理性证明

双序列比对的理论基础之建造替换矩阵的合理性证明

 前言：如果对最大似然估计没有概念的话，可以看看我之前写的《似然，似然，似是而然》
 结合前几篇文章我们大致的了解了计分矩阵的流程：对某以蛋白质家族进行多序列对比，然后按某一阈值(等同残基比）进行聚类，之后将匹配的无空位的区域划分为block,然后统计各个block中残基之间的联配的频率，用归一化的频率估计概率，进行最大似然估计，估计出在自然界中各残基联配的概率（即匹配模型M的参数）。
  《双序列比对的基础（2）之替换（计分）矩阵系列》提出了疑问：怎么没进行最大似然估计啊？没有列似然函数啊，没有求极值点啊。从样本数据得到的频率怎么直接估计为总体的参数呢？ 所以怀疑建造的替换矩阵是不合理的！
 那么本文就探讨探讨这个问题。
 首先，之前我们说过我们必须将生物学的问题抽象成数学模型。而残基对之间的联配可看做是多项分布。多项分布是二项分布的推广。
 二项分布就是我抛出一枚硬币，它的结果不是正面就是反面。我们抛出10次，计算出现5次正面的概率很简单。而多项分布则是我扔色子会有六种状态，现在我扔了十次，我想知道出现事件A（A事件=点数为1出现2次和点数为2出现4次和点数为3出现1次和点数为4出现1次和点数为5出现1次和点数为6出1次。）的概率是多少？注意在多项分布中，每个状态出现的频率必须大于0！
  而残基对之间的联配情况就可看做多项分布，只不过有210中状态（20种氨基酸的两两组合），并且在blocks中每种状态都出现过。所以我们就可以按照多项分布式列出似然函数，为进行最大似然估计奠定了基础。
 符号说明：
n={$n_1$、$n_2$、$n_3$、、、$n_{210}$}
$n_i$:第i种状态被观测到的频率。如 $n_1$可表示鸟氨酸与亮氨酸在blocks中联配的频率。
N:blocks中出现的残基对总数=$n_1$+ $n_2$+$n_3$+…+$n_{210}$
$\theta$={ ${\theta }_1$、 ${\theta }_2$、 ${\theta }_3$、、、 ${\theta }_{210}$}
${\theta }_i$=第i种状态出现的概率参数。
${\theta }^{ML}$=n/N（用归一化的频率表示概率）
 所以现在我们需要证明当$\theta$=${\theta }^{ML}$ 时，似然函数L（$n_1$,$n_2$,$n_3$,$n_{210}$|$\theta$)取得最大值。即可将 ${\theta }^{ML}$估计为在自然界中残基对联配和残基与空位联配的概率。
 下面我们将提供两种证明方法。分别用到了相对熵的性质和拉格朗日乘数法。
  第一种：此证明等价于任意$\theta$!=${\theta }^{ML}$,$\log \frac{L（n_1,n_2,n_3,,,n_{210}|{\theta }^{ML})}{L\left(n_1,n_2,n_3,,,n_{210}|\theta \right)}$ >0
 根据多向分布的概率函数可以轻易的写出对应的似然函数。（两函数形式上一样，只是看待的角度不一样） 因此
L($n_1$,$n_2$,$n_3$,$n_{210}$|$\theta$)=$\frac{N!}{{\prod }_{i=1}^{210}{n}_i!}$*${\prod }_{i=1}^{210}{\theta }_i^{n_i}$

从而，$\log \frac{L（n_1,n_2,n_3,,,n_{210}|{\theta }^{ML})}{L\left(n_1,n_2,n_3,,,n_{210}|\theta \right)}$=$\log \frac{{\prod }_{i=1}^{210}{\left({\theta }_i^{ML}\right)}^{n_i}}{{\prod }_{i=1}^{210}{\theta }_i^{n_i}}$=${\sum }_{i=1}^{210}{n}_i\log \frac{{\theta }_i^{ML}}{{\theta }_i}$
又因为
$n_i$/N=${\theta }_i^{ML}$
所以将$n_i$换掉，
${\sum }_{i=1}^{210}{n}_i\log \frac{{\theta }_i^{ML}}{{\theta }_i}$ = $N{\sum }_{i=1}^{210}{\theta }_i^{ML}\log \frac{{\theta }_i^{ML}}{{\theta }_i}$
我们可以发现：
${\sum }_{i=1}^{210}{\theta }_i^{ML}\log \frac{{\theta }_i^{ML}}{{\theta }_i}$是相对熵的形式，而根据相对熵大于等于0的性质，且N大于0，所以我们得证。
ps:

相对熵：$H(P\Vert Q)={\sum }_{i}P\left(x_i\right)\log \frac{P\left(x_i\right)}{Q\left(x_i\right)}$

（相对熵大于等于0的证明将放在另一篇博客。）
 第二种证明方法：
 最大似然估计本身就是多元函数求极值点。所以我们应该求出并证明极值点为${\theta }^{ML}$=n\N
  我们现在所已知的是多元函数（即上文的似然函数）和变量的限制条件：${\sum }_{i=1}^{210}{\theta }_i$=1
  多元函数、极值、约束条件，看到这些我们可以很快的联想到可以使用拉格朗日乘数法。(另一篇博客将介绍此法的背景应用)
  证明此法如下：
似然函数：L($n_1$,$n_2$,$n_3$,$n_{210}$|$\theta$)=$\frac{N!}{{\prod }_{i=1}^{210}{n}_i!}$*${\prod }_{i=1}^{210}{\theta }_i^{n_i}$
 构造Lagrange函数:
Lagrange(${\theta }_1$,${\theta }_2$,${\theta }_3$,${\theta }_{210}$,$\lambda$)=L($n_1$,$n_2$,$n_3$,$n_{210}$|$\theta$)-$\lambda$(${\sum }_{i=1}^{210}{\theta }_i$ -1)
 对似然函数取对数：
Lagrange(${\theta }_1$,${\theta }_2$,${\theta }_3$,${\theta }_{210}$,$\lambda$)=$\log N!-{\sum }_{i=1}^{210}{n}_i!+{\sum }_{i=1}^{210}{n}_i\log {\theta }_i$-$\lambda$(${\sum }_{i=1}^{210}{\theta }_i$ -1)
  然后对Lagrange函数求关于${\theta }_i$的偏导：
$\frac{\partial Lagrange\left({\theta }_1,\dots ,{\theta }_{210},\lambda \right)}{\partial {\theta }_i}$=$\frac{n_i}{{\theta }_i}-\lambda$
  取极值时，导数为0：
$\frac{n_i}{{\theta }_i}-\lambda$=0
  故 $\frac{n_i}{\lambda }={\theta }_i$
又因为 ${\sum }_{i=1}^{210}{\theta }_i={\sum }_{i=1}^{210}\frac{n_i}{\lambda }=1$
 故$\lambda$=N
  故${\theta }_i=\frac{n_i}{\lambda }=\frac{n_i}{N}$
  故可证：${\theta }_i=\frac{n_i}{N}={\theta }_i^{ML}$
 以上就是两种证明方法。写到这里，我自己明白了并最大似然估计出看此文的读者们明白了，233333333333。