双泊松分布的最大似然

Question

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

这个link展示了双泊松分布的函数： 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
## ...

双泊松分布的最大似然

maximum likelihood in double poisson distribution

r

mle