Beta分布和抗性基因频率的推导

工作以后发现学的更多的是实践方面的，有一定的局限性，更喜欢记笔记而不是写博客了。习惯了latex后，也越来越用不惯博客如此麻烦的调格式等，markdown没有完全支持latex吧，至少我的/section命令并没有得到相应解释。看了两篇关于抗性基因频率估计的论文，顺便写篇总结吧。

参数估计

参数估计有两种方法:
第一种是频率学派，也就是通过抽样调查，计算似然函数来选择特定参数值，如对期望的无偏估计。
第二种是贝叶斯学派，是利用贝叶斯推断来计算参数，认为参数不是固定值也是变量，给参数赋予先验分布，并使得先验与似然共轭，通过求后验均值来得到参数的估计。
由于贝叶斯学派的计算极其依赖于先验概率的假设，而先前的理论和技术的限制使得当时的人们不能很好的取得先验概率，因而贝叶斯的理论虽然出现的极早，但并没有得到人们的重视，也未能在生产实践中发挥重大作用。当进入20世纪中期，计算机技术逐步应用于数学分析领域，使得人们有能力通过大量的迭代计算来得到较好的先验概率，所以近几十年贝叶斯理论重新进入人们的视野，并在计算机应用领域做出了重要贡献。

后验分布函数的推导

以下x是试验结果中的正有效值，例如抽检样品中的不合格产品，抗性检测中的抗性个体数目等，θ是需要估计的参数，例如与前面例子对应的真实的产品不合格概率，抗性个体的概率等。
π(θ)是θ的概率函数，p(x)是事件x的似然函数，也就是在θ一定的情况下，使得事件发生概率最大的函数。

p(θ|x)=p(x,θ)p(x)=p(x|θ)π(θ)∫p(x|θ)π(θ)dθ

确定参数分布函数

在计算和实践中，会选取一定的似然函数，使得参数的先验分布和后验分布函数形式一致，这个时候又可以称此参数的先验分布函数是似然函数的共轭函数，例如beta分布函数（参数的先验函数）是二项式分布（似然函数)的共轭函数，此例以后证明。

事件似然函数和参数先验分布函数的对应关系

1.二项分布 — beta分布
2.Gaussian分布 — Gaussian分布
3.指数分布 — Gamma分布

例:Beta分布的推导

下面，利用上述理论实现Beta分布的推导。
首先,假设事件发生t次，事件发生的概率是θ,并且认为事件的发生服从二项分布，由上述推得θ的后验概率函数是:

π(θ|x)=p(x,θ)π(x)==p(x|θ)π(θ)∫p(x|θ)π(θ)dθ=(nr)θt(1−θ)n−tπ(θ)∫(nr)θt(1−θ)n−tπ(θ)dθ

化简:积分项中，二项分布系数(nr)与θ无关，可以消掉，同时注意分母是p(x)的全概率公式，结果是p(x)，也与θ无关，所以可以推得如下结果:

π(θ|x)=θt(1−θ)n−tπ(θ)p(x)∝θt(1−θ)n−tπ(θ)

为了使得后验分布和先验分布具有相同的函数形式，那么必须有如下的形式:

π(θ)∝θa−1(1−θ)b−1

这时候，容易想到参照标准beta分布：

f(x;α,β)=xα−1(1−x)β−1∫10uα−1(1−u)β−1du=Γ(α+β)Γ(α)Γ(β)xα−1(1−x)β−1∝xα−1(1−x)β−1

所以可以假设先验分布服从beta分布，由此得到:

π(θ)=Γ(a+b)Γ(a)Γ(b)θa−1(1−θ)b−1

将上述函数带入后验分布函数，又由于Γ(a+b)Γ(a)Γ(b)是常数，可得到如下关系:

π(θ|x)∝θt+a−1(1−θ)n+b−t−1

实践中，为了简便起见，设a=b=1,所以得到

π(θ|x)∝θt(1−θ)n−t

这时，可以得到后验分布的函数:

π(θ|x)=Γ(α+β)Γ(α)Γ(β)θα−1(1−θ)β−1

由上述两个方程得到系数的关系如下:

{t=α−1n−t=β−1⇒{α=t+1β=n−t+1

根据beta分布获取估计参数的期望和方差

⎧⎩⎨⎪⎪⎪⎪E(x)=αα+βVar(θ)=αβ(α+β)2(α+β+1)⇒⎧⎩⎨⎪⎪⎪⎪E(θ)=αα+β=t+1n+2Var(θ)=E(θ)(1−E(θ))n+3

实践:估计抗性基因频率

设培养的lines数目是n,检测出的有抗性基因的line数目是t。设抗性line出现的概率是ψ根据上述推导，可以计算出抗性line频率的期望是E(ψ)=t+1n+2.
注意的是由于祖父代每四个等位基因有一个抗性基因，那么才可能在F2代检测出来抗性,设祖父代抗性基因出现的频率是θ,有如下关系:
ψ=4θ
所以得到:

⎧⎩⎨⎪⎪⎪⎪E(ψ)=αα+β=t+1n+2Var(ψ)=E(ψ)(1−E(ψ))n+3⇒⎧⎩⎨⎪⎪⎪⎪E(4θ)=αα+β=t+1n+2⇒E(θ)=t+14(n+2)Var(4θ)=E(4θ)(1−E(4θ))n+3⇒Var(θ)=E(θ)(1−4E(θ))4(n+3)