双泊松分布中的最大似然估计

Question

双泊松分布中的最大似然估计

4

首先，我使用“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) #平均数 [1] 10.986 mean(x)/var(x) #离散度 [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)

出现错误：

ddoublepois(x, mu.new, gam.new)中的错误：s必须为正数

因此，我尝试了另一种方法，输入了双泊松分布函数的方程。

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（与初始值相同）表明此代码失败了。

有谁能告诉我如何正确进行双泊松分布的最大似然估计？

提前感谢。

- Miyazaki

在你的logl函数中，你应该在某个时候取对数。 - Dason

嗨Dason，我也尝试了这个-sum(ddoublepois(x, mu.new, gam.new, log=TRUE))，但依然返回相同的错误信息。 - Miyazaki

请注意，这不是“替代方案”。替代方案将是取非对数版本的乘积。但如果您不取对数，就不能简单地求和。 - Dason

谢谢您的注意，Dason。 - Miyazaki

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

你的问题是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
## ...