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

首先，我使用“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)

头(x)

[1] 11 9 10 13 6 8

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

以下是参数的真实值:

均值(x)＃平均值

[1] 10.986

均值(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(开始= c(0.1,0.1)，下= 0，上= Inf，log12)

$ par

[1] 0.1 0.1

$目的

[1] Inf

$收敛

[1] 0

$次

[1] 1

$评估

功能梯度

2 4

$消息
[1]“X收敛(3)”

无疑，估计参数0.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
## ...

关于r - 双泊松分布的最大似然，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/50203075/