R语言学习笔记（八） -极大似然估计

欢迎大家来到IT世界,在知识的湖畔探索吧!

导语：首先继续跟大家一起学习一份R语言资料：

https://web.stanford.edu/class/bios221/book/index.html
。前两期跟大家介绍了常见的离散分布，本期跟大家一起学习极大似然估计。

1. 极大似然思想

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

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

欢迎大家来到IT世界,在知识的湖畔探索吧!

这个概率看起来很低，但如果我们推测池塘里鲤鱼和草鱼比例为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))) + }

欢迎大家来到IT世界,在知识的湖畔探索吧!

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

欢迎大家来到IT世界,在知识的湖畔探索吧!> 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二项分布的极大似然估计

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

欢迎大家来到IT世界,在知识的湖畔探索吧!> cb = c(rep(0, 110), rep(1, 10)) > table(cb) cb 0 1 110 10 > mean(cb) [1] 0.0

#设置一系列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的均值：

欢迎大家来到IT世界,在知识的湖畔探索吧!> cbmax [1] 0.083 > mean(cb) [1] 0.0

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的极大似然估计。

欢迎大家来到IT世界,在知识的湖畔探索吧!> 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.

4. 小结

本期给大家介绍了极大似然估计的思想和计算方法。简单来说就是写出似然函数，取对数求导，令导数等于0，最后求出参数的极大似然估计。

R语言学习笔记（二） ​​​R语言学习笔记（五）——曼哈顿图 R语言学习笔记（七） -离散型数据的模型预测2 R语言学习笔记（四）—pheatmap R语言学习笔记（三）——实用的内置函数 R语言学习笔记（六） -离散型数据的模型预测1