张金龙的博客分享 http://blog.sciencenet.cn/u/zjlcas 物种适应性、分布与进化

博文

二项分布的对数似然面:算法与R语言实现

已有 12985 次阅读 2011-11-19 00:47 |个人分类:统计分析|系统分类:科研笔记| 二项分布, 极大似然, 对数似然函数

二项分布的对数似然面:算法与R语言实现
jinlongzhang01@gmail.com

对数似然面是极大似然估计的结果之一,这里介绍在R语言中如何绘制对数似然面(Log Likelihood Surface),这里用到的对数似然函数是二项分布的,对于泊松分布,正态分布,负二项分布等,只需对对数似然函数进行适当修改即可。

 
绘制二项分布的对数似然面举例:
主要过程
(1) 写出对数似然函数
(2) 用nlm求对数似然函数的估计值
(3) 若绘制似然曲线,则固定一个参数值,并变化另外一个参数,

    计算Loglikelihood,用lines添加曲线。

(4) 若绘制趋势面,则先要计算参数组合矩阵,并求各栅格处的

      LogLikelihood,用persp, image, contour等函数即可。也

      可以用lattice函数包的wireframe()

以下的例子给出全部的R代码和相应说明。
 
该例子引自:
http://www.unc.edu/courses/2010fall/ecol/563/001/docs/lectures/lecture10.htm#Rcode
 

## 实例: 假设有以下数据, 记录每个枝条上蚜虫的数量,该数据出自

##Krebs, Charles. 1999. Ecological Methodology.

## Menlo Park, CA: Addison-Wesley. 130

###########################################################

#######################原始数据############################

###########################################################

### 数据导入

num.stems <- c(6,8,9,6,6,2,5,3,1,4)

aphids_count <- c(0,1,2,3,4,5,6,7,8,9)

aphids <- data.frame(aphids_count, num.stems)

### Raw Data,整理出二项分布拟合所需要的格式

aphids_raw <- rep(aphids_count, num.stems)

 

###########################################################

####################写出对数似然函数#######################

###########################################################

 

## 二项分布的对数似然函数,用来绘制对数似然曲线和趋势面

nnominb.loglik <- function(p) {

    res <- sum(log(dnbinom(aphids_raw, mu = p[1], size = p[2])))

    return(res)

}

## 二项分布的负对数似然函数,nlm()函数只能求最小值,因此极大

## 似然值时必须用到此函数

neg.nnominb.loglik <- function(p) {

    res <- sum(log(dnbinom(aphids_raw, mu = p[1], size = p[2])))

    return(-res)

}

###########################################################

###################估计对数似然函数极值####################

###########################################################

#利用nlm估计极值

out.NB <- nlm(neg.nnominb.loglik, c(4,4))

# estimate 是对两个参数的估计,即对neg.nnominb.loglik函数参数p的估计

out.NB

###########################################################

################### 1 ##绘制似然曲线#####################

###########################################################

 

### x 序列

x <- seq(0.1, 10, by = 0.1)

### 计算y

y_max <- sapply(x, function(x) nnominb.loglik(p = c(out.NB$estimate[1], x)))

plot(y_max ~ x, type = "l", ylab = "Log Likelihood", xlab = "Size")

### mu = 2

y_2 <- sapply(x, function(x) nnominb.loglik(p = c(2, x)))

lines(y_2 ~ x, col = 2)

### mu = 5

y_5 <- sapply(x, function(x) nnominb.loglik(p = c(5, x)))

lines(y_5 ~ x, col = 3)

### 添加图例

legend('bottomright', c(expression(mu=='MLE (3.46)'),

        expression(mu==2), expression(mu==5)),

        col=c(1,2,3), lty=c(1,1,1), cex=.9, bty='n')

###########################################################

################## 2#绘制栅格图##########################

 

### 绘制栅格图

## 按照值设色

image(seq(2,5,.1), seq(1,4,.1), z.matrix, col = cm.colors(30),

      xlab = expression(mu), ylab = "Number of Trials")

## 添加等高线

contour(seq(2,5,.1), seq(1,4,.1), z.matrix, add = TRUE)

## 极大似然值

points(x = out.NB$estimate[1], y = out.NB$estimate[2], pch = 3)

###########################################################

############### 3绘制3D 对数似然面#######################

###########################################################

 

### 计算趋势面原始数据

### 构建网格

grids <- expand.grid(mu=seq(2,5,.1), theta=seq(1,4,.1))

### 计算z

grids$z <- apply(grids, 1, function(x) nnominb.loglik(p = c(x[1], x[2])))

### 趋势面原始数据

z.matrix <- matrix(grids$z, nrow=length(seq(2,5,.1)))

 

### 绘制三维趋势面

persp(seq(2,5,.1), seq(1,4,.1), z.matrix,

      xlab=expression(mu), ylab=expression(theta), zlab='Log Likelihood',

      ticktype='detailed', theta=30, phi=20, d = 4)

###########################################################

########## 4##绘制对数似然面的三维设色图#################

 

library(lattice)       

wireframe(z~mu*theta, data=grids, xlab=expression(mu),

          ylab=expression(theta), zlab="Log-nlikelihood",

          scales = list(arrows = FALSE), drape=TRUE,

          par.settings=list(par.zlab.text=list(cex=.75),

          axis.text=list(cex=.75)))

 




https://wap.sciencenet.cn/blog-255662-509501.html

上一篇:用R语言绘制似然面Likelihood Surface
下一篇:一个简单的马尔可夫链
收藏 IP: 213.21.84.*| 热度|

0

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

数据加载中...
扫一扫,分享此博文

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

GMT+8, 2024-12-27 08:19

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部