1
votes

Has anyone had problems with the Weibull distribution using the ExtDist Package?

From the documentation:

Parameter Estimation for a distribution with unknown shape parameters Example from: Rinne (2009) Dataset p.338 and example pp.418-419 Parameter estimates are given as shape = 99.2079 and scale = 2.5957. The log-likelihood for this data and Rinne's parameter estimates is -1163.278.

data <- c(35,38,42,56,58,61,63,76,81,83,86,90,99,104,113,114,117,119,141,183)
est.par <- eWeibull(X=data, method="numerical.MLE"); est.par
plot(est.par)

However when I run this I get the following output:

Parameters for the Weibull distribution.
(found using the  numerical.MLE method.)

 Parameter  Type   Estimate       S.E.
     shape shape 5.82976007 1.79326460
     scale scale 0.06628166 0.02129258

This is clearly wrong but I am not sure if I have made a mistake or if there is a bug in the package?

1
did you read the next bit of the example? "# Estimates calculated by eWeibull differ from those given by Rinne(2009). However, eWeibull's parameter estimates appear to be an improvement, due to a larger log-likelihood of -99.09037 (as given by lWeibull below)." - Ben Bolker
Hi Ben, yes I did but the answer I got from the code of 5.82976007 is definitely not even close, the output should be in the region of 99 - James

1 Answers

1
votes

It seems to me it's a bug in the package. I did my own independent MLE and got the same answer as Rinne:

library(bbmle)
m1 <- mle2(y~dweibull(shape=exp(lshape),scale=exp(lscale)),
     data=data.frame(y=data),
     start=list(lshape=0,lscale=0))

Then I dug in and looked at the source of the dWeibull function:

function (x, shape = 2, scale = 2, params = list(shape = 2, scale = 2)) 
{
    if (!missing(params)) {
        shape <- params$shape
        scale <- params$scale
    }
    out = stats::dgamma(x, shape, scale)
    return(out)
}

It seems that out should be set to the result of dweibull(...) rather than dgamma(...) ... ?? Looking at the rest of the weibull code, this error seems to be repeated -- maybe this is just a sloppy cut-and-paste? I would definitely contact the maintainer (maintainer("ExtDist")).

PS. If I fit a Gamma distribution using my alternative method I get exactly the same answers as the ExtDist package:

m1g <- mle2(y~dgamma(shape=exp(lshape),rate=exp(lrate)),
     data=data.frame(y=data),
     start=list(lshape=0,lrate=0))
exp(coef(m1g))
##     lshape      lrate 
## 5.82976007 0.06628166