爱奇百科 手机版
您的位置: 首页 > 常识 >

哈迪温伯格定律是什么(R语言学习笔记 -极大似然估计)

100次浏览     发布时间:2024-11-13 10:19:47    


1. 极大似然思想


还是举个例子说明,假如一个池塘里面有鲤鱼和草鱼,想知道这两种鱼的比例,所以打捞了100条鱼,发现有80条鲤鱼、20条草鱼。现在如果让你估计池塘里两种鱼的比例,你会很自然认为有80%的鲤鱼,20%的草鱼。那么这样推理的逻辑是什么呢?其实这背后就是极大似然的思想。


对上面的例子我们可以再深入思考一下,根据我们认为池塘里的鲤鱼和草鱼比例为4:1,那么我们打捞100条鱼,出现80条鲤鱼、20条草鱼的概率是多少呢?这个其实是我们之前提到的二项分布的例子,其概率等于:



这个概率看起来很低,但如果我们推测池塘里鲤鱼和草鱼比例为1:1呢,再计算一下出现80条鲤鱼、20条草鱼的概率:



看出差别了吧,当我们推测鲤鱼和草鱼比例为4:1时,出现80条鲤鱼、20条草鱼的概率其实就是最大的。


现在总结一下这个过程,一个随机试验如有若干个可能的结果A,B,C,…。若在一次试验中,结果A出现,则一般认为试验条件对A出现有利,也即A出现的概率很大,这也就是极大似然的思想。


2. 极大似然估计


求极大似然估计值的一般步骤:


(1)写出似然函数;


我们通常假定是独立同分布的样本,因此观察到这些样本的概率就是各样本概率的乘积,即L(θ) =


(2)对似然函数取对数,并整理;


这一步其实是为了简化运算,将乘积变为求和形式。


(3)求导数,解似然方程。


对似然函数求一阶导数,一般使一阶导数为0的θ值,就是我们要估计的极大似然估计值。


下面分别以泊松分布和二项分布为例,用R代码计算。


2.1泊松分布的极大似然估计


现在我们有一个泊松分布的样本e100,要估计泊松分布的参数λ。按照上面的步骤,首先写出似然函数:



用R写出这个函数,并取对数:


> loglikelihood  =  function(lambda, data = e100) {
+   sum(log(dpois(data, lambda)))
+ }


现在我们设置一系列λ值,分别计算似然函数值,并画出似然函数与λ的图像,观察λ取值多少时,似然函数最大:


> lambdas = seq(0.05, 0.95, length = 100)   #设置一系列λ值
> loglik = vapply(lambdas, loglikelihood, numeric(1))  #计算似然函数值
> plot(lambdas, loglik, type = "l", col = "red",ylab = "", lwd = 2,
+      xlab = expression(lambda))
> m0 = lambdas[which.max(loglik)]       #似然函数取最大值时的λ值
> abline(v = m0, col = "blue", lwd = 2)
> abline(h = loglikelihood(m0), col = "purple", lwd = 2)
> m0
[1] 0.55



从上面的计算中我们可以看到λ取0.55时,似然函数取最大值。我们再计算一下e100数据的样本均值:


> mean(e100)
[1] 0.55


竟然也是0.55!这是巧合吗?当然不是!因为泊松分布λ的极大似然估计量就是样本均值。在这里就不给出严格的证明了,但通过上面的例子我们可以清楚看到λ取样本均值时,似然函数取最大值。


2.2二项分布的极大似然估计


还是以一个简单的数据为例说明:


> cb  =  c(rep(0, 110), rep(1, 10))
> table(cb)
cb
  0  1
110  10
> mean(cb)
[1] 0.08333333


#设置一系列p值,分别计算似然函数值


> probs = seq(0, 0.3, by = 0.001)
> likelihood = dbinom(sum(cb), prob = probs, size = length(cb))
> plot(probs, likelihood, pch = 16, xlab = "probability ofsuccess",
+      ylab = "likelihood",cex=0.6)
> cbmax <- probs[which.max(likelihood)]
> abline(v = cbmax, col = "red", lwd = 2)



二项分布参数p的极大似然估计也是样本均值,可以看到使似然函数取极大值的点就是数据集cb的均值:


> cbmax
[1] 0.083
> mean(cb)
[1] 0.08333333


3. 多项式分布的极大似然估计


这里介绍一个数量遗传学中的多项式分布,哈迪-温伯格平衡(hardy-weinberg equilibrium)。该定律是说在一个大且随机交配的种群中,基因频率和基因型频率在没有迁移、突变和选择的条件下会保持不变。此时若等位基因M的频率为p,等位基因N的频率为q = 1-p,则该基因位点三种基因型频率为:


pMM = p2,pMN = 2pq, pNN = q2


nMM、nMN、nNN分别表示三种基因型的个体数,S = nMM + nMN + nNN表示总个体数,现在我们要计算p的极大似然估计。


还是先写出似然函数:



对该函数取对数,计算极值,得到的就是p的极大似然估计:



下面用代码说明,这里用到了HardyWeinberg包中的Mourant数据集,该数据集包含216个群体MN血型系统的基因型频率。


> install.packages("HardyWeinberg")       #安装HardyWeinberg包
> library("HardyWeinberg")
> data("Mourant")
> Mourant[214:216,]
    Population    Country Total  MM MN  NN
214    Oceania Micronesia   962 228 436 298
215    Oceania Micronesia   678 36 229 413
216    Oceania     Tahiti  580 188 296  96


可以看到Mourant数据集有5列,前2列是群体信息,后4列是基因型频率。这里以第216个群体为例,计算等位基因M频率p的极大似然估计。


> nMM = Mourant$MM[216]
> nMN = Mourant$MN[216]
> nNN = Mourant$NN[216]
> loglik = function(p, q = 1 - p) {             #似然函数取对数
+   2 * nMM * log(p) + nMN *log(2*p*q) + 2 * nNN * log(q)
+ }
> xv = seq(0.01, 0.99, by = 0.0001)
> yv = loglik(xv)
> #ggplot画图
> library(ggplot2)
> plot_lik <- data.frame(xv, yv)
> p <- ggplot(data = plot_lik, aes(x = xv, y = yv))
> p+geom_point()+
+   geom_vline(xintercept =xv[imax],lwd = 1, col = "blue")+
+   geom_hline(yintercept =yv[imax],lwd = 1, col = "red")



我们验证一下我们上面的结论:


> xv[imax]
[1] 0.5793
> (nMM+nMN/2)/(nMM+nMN+nNN)
[1] 0.5793103

相关文章

  • 沭阳太阳能售后故障维修客服服务热线是多少实时反馈-今-日-汇-总

    沭阳太阳能24小时售后客服热线:400-883-2086沭阳太阳能全国统一售后维修电话:400-883-2086沭阳太阳能的技术团队是一支由业界精英组成的强大阵容,他们在安全机构的研发设计方面有着深厚的造诣,多项创新技术已申请专利,为公司的持续发展注入源源不断的动力。更重要的是,沭阳太阳

    2025-05-09 12:42:56
  • 意大利优尼卡壁挂炉全国售后服务网点号码实时反馈-今-日-汇-总

    意大利优尼卡壁挂炉产品使用介绍:便捷安全,操作无忧意大利优尼卡壁挂炉作为领域的佼佼者,凭借精湛工艺与前沿科技,为用户珍贵财物保驾护航。当您选择意大利优尼卡壁挂炉,全国售后服务点热线号码 400-883-2086 随时待命,无论使用中遇到任何疑问,都能获得专业解答。下面为您详细介绍意大利优尼卡壁挂

    2025-05-09 12:42:05
  • 森太热水器售后号码多少实时反馈-今-日-更-新

    森太热水器售后服务电话400-883-2086越来越多的人意识到家财安全的重要性,为了保护贵重物品的安全,森太热水器成为了很多家庭的必备品。而森太热水器作为国内知名品牌,其专业服务备受广大消费者的赞誉。森太热水器专业服务电话(以产品说明书或保修卡电话为准),全方位解决您的家财安全问题。专

    2025-05-09 12:40:37
  • 巴特利锅炉全国各市售后热线号码实时反馈-今-日-更-新

    巴特利锅炉:保障您的财产安全在现代社会,人们对于财产安全的需求也日益增加,在人们的生活中扮演着不可或缺的角色。巴特利锅炉作为国内领先的品牌,致力于为广大消费者提供高品质、可靠的产品和服务。巴特利锅炉24小时售后客服热线:400-883-2086了解巴特利锅炉:拨打全国统一专业维修

    2025-05-09 12:39:13
  • 萨弘壁挂炉服务热线号码各区24小时维修实时反馈-今-日-汇-总

    萨弘壁挂炉售后服务电话:400-883-2086萨弘壁挂炉24小时维修热线:400-883-2086萨弘壁挂炉,作为业内的践行者,不仅专注于高品质的研发与制造,更将客户服务体验提升至战略高度。我们的维修团队由一群技术精湛、素质卓越的专业人员组成,他们不仅拥有丰富的现场处理经验和扎实的理论

    2025-05-09 12:39:02

热门文章

最新文章