Gompertz-Makeham parameter estimation

I would like estimate the parameters of the Gompert-Makeham distribution, but I haven't got a result. I would like a method in R, like this Weibull parameter estimation code:

weibull_loglik <- function(parm){
   gamma <- parm[1]
   lambda <- parm[2]
   loglik <- sum(dweibull(vec, shape=gamma, scale=lambda, log=TRUE))
   return(-loglik)
}
weibull <- nlm(weibull_loglik,parm<-c(1,1), hessian = TRUE, iterlim=100)

weibull$estimate
c=weibull$estimate[1];b=weibull$estimate[2]

My data:

  [1]  872   52   31   26   22   17   11   17   17    8   20   12   25   14   17
 [16]   20   17   23   32   37   28   24   43   40   34   29   26   32   34   51
 [31]   50   67   84   70   71  137  123  137  172  189  212  251  248  272  314
 [46]  374  345  411  494  461  505  506  565  590  535  639  710  733  795  786
 [61]  894  963 1019 1149 1185 1356 1354 1460 1622 1783 1843 2049 2262 2316 2591
 [76] 2730 2972 3187 3432 3438 3959 3140 3612 3820 3478 4054 3587 3433 3150 2881
 [91] 2639 2250 1850 1546 1236  966  729  532  375  256  168  107   65   39   22
[106]   12    6    3    2    1    1

summary(vec)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0    32.0   314.0   900.9  1355.0  4054.0

Solution

It would be nice to have a reproducible example, but something like:

library(bbmle)
library(eha)
set.seed(101)
vec <- rmakeham(1000, shape = c(2,3), scale = 2)
dmwrap <- function(x, shape1, shape2, scale, log) {
   res <- try(dmakeham(x, c(shape1, shape2), scale, log = log), silent = TRUE)
   if (inherits(res, "try-error")) return(NA)
   res
}
m1 <- mle2(y ~ dmwrap(shape1, shape2, scale), 
          start = list(shape1=1,shape2=1, scale=1),
       data = data.frame(y = vec),
       method = "Nelder-Mead"
)

Define a wrapper that (1) takes shape parameters as separate values; (2) returns NA rather than throwing an error when e.g. parameters are negative
Use Nelder-Mead rather than default BFGS for robustness
the fitdistrplus package might help too
if you're going to do a lot of this it may help to fit parameters on the log scale (i.e. use parameters logshape1, etc., and use exp(logshape1) etc. in the fitting formula)

I had to work a little harder to fit your data; I scaled the variable by 1000 (and found that I could only compute the log-likelihood; the likelihood gave an error that I didn't bother trying to track down). Unfortunately, it doesn't look like a great fit (too many small values).

x <- scan(text = "872   52   31   26   22   17   11   17   17    8   20   12   25   14   17
   20   17   23   32   37   28   24   43   40   34   29   26   32   34   51
   50   67   84   70   71  137  123  137  172  189  212  251  248  272  314
  374  345  411  494  461  505  506  565  590  535  639  710  733  795  786
  894  963 1019 1149 1185 1356 1354 1460 1622 1783 1843 2049 2262 2316 2591
 2730 2972 3187 3432 3438 3959 3140 3612 3820 3478 4054 3587 3433 3150 2881
2639 2250 1850 1546 1236  966  729  532  375  256  168  107   65   39   22
12    6    3    2    1    1")
m1 <- mle2(y ~ dmwrap(shape1, shape2, scale), 
          start = list(shape1=1,shape2=1, scale=10000),
       data = data.frame(y = x/1000),
       method = "Nelder-Mead"
)

cc <- as.list(coef(m1))
png("gm.png")
hist(x,breaks = 25, freq=FALSE)
with(cc,
     curve(exp(dmwrap(x/1000, shape1, shape2, scale, log = TRUE))/1000, add = TRUE)
     )
dev.off()