I'm trying to get the shape and scale parameters for this data using the optim function in R.
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)))
The initial parameters were calculated using the method of moments
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
This is the code I used to run the optimization
# 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)
I get these results
$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"
When I test the results with those parameters the values are too low and I can't plot the distribution nor the likelihood function and it doesn't make sense to me. This is supposed to give the proability of falling in a particular income interval. What I'm doing wrong?