如何解决在 R 中使用 optim 进行伽马分布的最大似然估计
我正在尝试使用 R 中的 optim 函数获取此数据的形状和比例参数。
incomeData = data.frame(L = c(850,rep(1000,24),rep(2001,112),rep(3001,267),rep(4001,598),rep(5001,1146)),U = c(999,rep(2000,rep(3000,rep(4000,rep(5000,rep(10000,Interval = c(1,rep(2,rep(3,rep(4,rep(5,rep(6,1146)))
初始参数采用矩量法计算
incomeData$middle = (incomeData$U+incomeData$L)/2 # middle point of the interval
middlePointMean = mean(incomeData$middle) # mean of the middle points
middlePointvar = var(incomeData$middle) # variance of the middle points
initialPar1 = middlePointvar/(middlePointMean^2) # initial shape parameter (this was suggested)
initialPar2 = initialPar1/middlePointMean # initial scale parameter
这是我用来运行优化的代码
# The likelihood function for this problem is defined by the product of the difference between the
# cumulative gamma evaluated in the upper bound of the interval - the cumulative gamma evaluated in
# the lower bound of the interval.
logLikelihood = function(par){
ub = incomeData$U
lb = incomeData$L
# I'm applying sum instead of prod since the log of a product would be the sum
logLike = sum(pgamma(ub,shape = par[1],scale = par[2]) -
pgamma(lb,scale = par[2]))
return(-logLike)
}
optim(par = c(initialPar1,initialPar2),fn = logLikelihood,method = "L-BFGS-B",lower = 0.00001,upper=.99999)
我得到了这些结果
$par
[1] 1.014180e-01 1.737418e-05
$value
[1] 0
$counts
function gradient
1 1
$convergence
[1] 0
$message
[1] "CONVERGENCE: norM OF PROJECTED GRADIENT <= PGTOL"
当我用这些参数测试结果时,这些值太低,我无法绘制分布或似然函数,这对我来说没有意义。这应该给出在特定收入区间下降的概率。我做错了什么?
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。