r - 双泊松分布中的最大似然

Question

首先我使用“rmutil”包来模拟双泊松分布数据。泊松和双泊松的区别在于，双泊松允许过分散和欠分散，其中均值和方差不必相等。

此链接显示了双泊松分布的函数:
http://ugrad.stat.ubc.ca/R/library/rmutil/html/DoublePoisson.html

我模拟了一组大小为 500 的数据。

set.seed(10)
library("rmutil")

nn = 500 #size of data
gam = 0.7 #dispersion parameter
mu = 11

x <- rdoublepois(nn, mu, gam)

head(x)

[1] 11 9 10 13 6 8

 mean(x) #mean
 mean(x)/var(x) #dispersion

以下是参数的真实值:

mean(x) #mean

[1] 10.986

mean(x)/var(x) #dispersion

[1] 0.695784

为了通过 MLE 获取参数，我使用了 nlminb 函数来最大化对数似然函数。对数似然函数由“rmutil”包中的双重分布的密度函数形成。

logl <- function(par) {
  mu.new <- par[1]
  gam.new <- par[2]

  -sum(ddoublepois(x, mu.new, gam.new, log=TRUE))
 }
nlminb(start = c(0.1,0.1), lower = 0, upper = Inf, logl)

结果出现错误:

Error in ddoublepois(x, mu.new, gam.new) : s must be positive

所以我再次尝试，我输入了双泊松密度函数的方程。

logl2 <- function(par) {
  mu.new <- vector() #mean
  gam.new <- vector() #dispersion
  ddpoi <- vector()


for (i in 1:nn){    
    ddpoi[i] <- 0.5*log(gam.new[i])-gam.new[i]*mu.new[i]
    +x[i]*(log(x[i])-1)-log(factorial(x[i]))
    +(gam.new[i])*x[i]*(1+log(mu.new[i]/x[i]))
  }
  -sum(ddpoi)
 }
nlminb(start = c(0.1,0.1), lower = 0, upper = Inf, logl2)

输出:

nlminb(start = c(0.1,0.1), lower = 0, upper = Inf, logl2)

$par

[1] 0.1 0.1

$objective

[1] Inf

$convergence

[1] 0

$iterations

[1] 1

$evaluations

function gradient

2 4

$message [1] "X-convergence (3)"

当然，估计参数为 0.1(与初始值相同)表明此代码失败。

任何人都可以告诉我如何对双泊松分布进行正确的最大似然估计？

提前致谢。

Answer 1

你的问题是 nlminb正在尝试评估边界上的函数(即 s 完全等于 0)。

解决这个问题的一种方法是修改 logl包括调试语句:

logl <- function(par,debug=FALSE) {
    mu.new <- par[1]
    gam.new <- par[2]
    if (debug) cat(mu.new,gam.new," ")
    r <- -sum(ddoublepois(x, m=mu.new, s=gam.new,log=TRUE))
    if (debug) cat(r,"\n")
    return(r)
}
nlminb(start = c(0.1,0.1), lower = 0, upper = Inf, logl, debug=TRUE)
## 0.1 0.1 3403.035 
## 0.1 0.1 3403.035 
## 0.1 0.1 3403.035 
## 1.022365 0 Error in ddoublepois(x, m = mu.new, s = gam.new, log = TRUE) : 
## s must be positive

现在尝试将边界从零稍微偏移:

nlminb(start = c(0.1,0.1), lower = 1e-5, upper = Inf, logl)

给出合理的答案

## $par
## [1] 10.9921451  0.7183259
## ...