双泊松分布中的极大似然
双泊松分布中的极大似然
提问于 2018-05-06 18:18:18
首先，使用"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) 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(start = c(0.1,0.1)，Inf= 0，Inf= Inf，logl2) $par 1 0.1 0.1 $objective 1 Inf $convergence 1%0 $iterations 1 1 $evaluations 函数梯度 2 4 $message 1“X-收敛(3)”
当然，估计参数0.1 (与初始值相同)显示此代码失败。
谁能教我如何对双泊松分布做正确的最大似然估计？
提前谢谢。
回答 1
Stack Overflow用户
回答已采纳
发布于 2018-05-06 19:27:29
您的问题是，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
## ...页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/50203075
复制相关文章
相似问题
腾讯云开发者
