高山流水分享 http://blog.sciencenet.cn/u/friendpine 走在科学路上的一位无名侠客,只是静静的走着.........

博文

R语言中的最大似然估计

已有 22466 次阅读 2011-5-12 16:37 |个人分类:统计学与R语言学习|系统分类:科研笔记

对于最大似然估计,都是依赖于似然函数的,因此关键在于写出似然函数,然后对该似然函数进行优化。似然函数依个人的问题而定,比如说下面的是正态分布的似然函数:

normal <- function(theta,x){
  mu <- theta[1]
  sigma2 <- theta[2]
  n <- length(x)
  logL <- -0.5*n*log(2*pi)-0.5*n*log(sigma2)-(1/(2*sigma2))*sum((x-mu)**2)
  return (-logL)
}
上面的theta是指需要估计的正态分布的均值和方差,logL是似然值,之所以返回负数是因为后面用到的优化似然值的函数是最小化函数。
写出来了似然函数之后,下一步就是如何优化似然函数的值。给定数据x是一个符合正态分布的数据,函数optim()可以进行后续的优化。该函数的一般形式如下:
optim(initial valurs of theta, likelihood function,data)

因此,我们这里的优化结果就可以通过下面的语句给出来:

x <- rnorm(100)

result <- optim(c(0,1),normal,x=x)
result中包含对于theta的估计值,以及优化之后的似然值,其他返回值可以看该函数的文档

也可以采用其他的函数进行估计,比如maxLik包是一个专门用来进行最大似然估计的,其中的maxLik()函数也可以对于给定的似然函数进行优化,比如说上面的正态分布似然函数normal(),可以直接输入到maxLik()函数中进行估计,不过需要注意的是该函数默认是正的似然函数值,而且不是在函数中输入数据,因此normal()函数需要做一些修改如下:
normal <- function(theta){
   mu <- theta[1]
   sigma <- theta[2]
   logL <- -0.5*N*log(2*pi) - N*log(sigma) - sum(0.5*(x - mu)^2/sigma^2)
   return (logL)
}
对于给定的数据 x <- rnorm(100,1,2),N <- length(x),然后就可以使用下面的语句进行最大似然估计了:
result <- maxLik(normal,start=c(0,1))
得到的结果如下:
print(result)
Maximum Likelihood estimation
Newton-Raphson maximisation, 8 iterations
Return code 1: gradient close to zero
Log-Likelihood: -2117.389 (2 free parameter(s))
Estimate(s): 1.007240 2.010635
可见,虽然在maxLik()函数中给出的初始值是错误的,但是最后得到的估计值非常接近真值。


http://blog.sciencenet.cn/blog-54276-443431.html

上一篇:关于回归的一些问题
下一篇:对于流行病学建模的粗浅认识-1

2 田灿荣 黄俊

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-2-25 18:05

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部