r - 使用 R 中的 optim 对伽马分布进行最大似然估计

Question

我正在尝试使用 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,24),rep(3000,112),rep(4000,267),rep(5000,598),rep(10000,1146)),
       Interval = c(1,rep(2,24),rep(3,112),rep(4,267),rep(5,598),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,shape = par[1],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"

当我用这些参数测试结果时，值太低，我无法绘制分布或似然函数，这对我来说没有意义。这应该给出在特定收入区间下降的可能性。我做错了什么？

标签： rstatistics