restimation

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"
    )
    

    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()
    

    enter image description here