用R语言做数据分析——点估计之极大似然法

最大似然估计提供了一种给定观察数据来评估模型参数的方法，即：“模型已定，参数未知”。简单而言，假设我们要统计全国人口的身高，首先假设这个身高服从服从正态分布，但是该分布的均值与方差未知。我们没有人力与物力去统计全国每个人的身高，但是可以通过采样，获取部分人的身高，然后通过最大似然估计来获取上述假设中的正态分布的均值与方差。

最大似然估计中采样需满足一个很重要的假设，就是所有的采样都是独立同分布的。下面我们具体描述一下最大似然估计：

首先，假设X1,X2,...,Xn为独立同分布的采样，θ为模型参数,f为我们所使用的模型，遵循我们上述的独立同分布假设。参数为θ的模型f产生上述采样可表示为

回到上面的“模型已定，参数未知”的说法，此时，我们已知的为X1,X2,...,Xn，未知为θ，故似然定义为:

在实际应用中常用的是两边取对数，得到公式如下：

其中左侧公式被称为对数似然，右侧公式称为平均对数似然。而我们平时所称的最大似然为最大的对数平均似然，即：

接下来我们以R语言“MASS”包中的geyser数据，用极大似然法拟合模型曲线。

>#引入MASS包， 使用geyser数据

> library(MASS)

> attach(geyser)

>#绘制数据的频率分布直方图

> hist(waiting)

从图中可以看出，其分布是两个正态分布的混合。可以用如下的分布函数来描述该数据：

该函数中有5个参数[p、μ1、σ1、μ2、σ2]需要确定。上述分布函数的对数极大似然函数为：

在R语言中定义对数似然函数，保存文件为“log.R”：

LL<-function(params,data){

#参数"params"是一个向量，依次包含了五个参数：p,mu1,sigma1,mu2,sigma2.

#参数"data"，是观测数据。

t1<-dnorm(data,params[2],params[3])

t2<-dnorm(data,params[4],params[5])

#这里的dnorm()函数是用来生成正态密度函数的

f<-params[1]*t1+(1-params[1])*t2

#混合密度函数

ll<-sum(log(f))

return(-ll)

}

引入对数似然函数，进行参数估计：

> source('log.R')

> #用hist函数找出初始值

> hist(waiting,freq = F)

> lines(density(waiting))

>#拟合函数

>geyser.res<-nlminb(c(0.5,50,10,80,10),LL,data=waiting,lower=c(0.0001,-Inf,0.0001,-Inf,-Inf,0.0001),upper = c(0.9999,Inf,Inf,Inf,Inf))

> #初始值为p=0.5,mu1=50,sigma1=10,mu2=80,sigma2=10

> #LL是被最小化的函数。

> #data是拟合用的数据

> #lower和upper分别指定参数的上界和下界。

下一步进行结果的估计

> #查看拟合的参数

> geyser.res$par

[1] 0.3075937 54.2026518 4.9520026 80.3603085 7.5076330

> #拟合的效果

> X<-seq(40,120,length=100)

> #读出估计的参数

> p<-geyser.res$par[1]

> mu1<-geyser.res$par[2]

> sig1<-geyser.res$par[3]

> mu2<-geyser.res$par[4]

> sig2<-geyser.res$par[5]

> #将估计的参数函数代入原密度函数。

> f<-p*dnorm(X,mu1,sig1)+(1-p)*dnorm(X,mu2,sig2)

> #作出数据的直方图

> hist(waiting,probability=T,col=0,ylab="Density",

+ ylim=c(0,0.04),xlab="Eruption waiting times")

> #画出拟合的曲线

> lines(X,f)

>#释放数据

> detach()

由上可知最大似然估计的一般求解过程：

1. 写出似然函数；

2. 对似然函数取对数，并整理；

3. 求导数 ；

4. 解似然方程