I'm trying to estimate the Weibull-Gamma Distribution parameters, but I'm encountering the following error:
"the function mle failed to estimate the parameters, with the error code 7"
What do I do?
The Weibull-Gamma Distribution
Density Function
dWeibullGamma <- function(x, alpha, beta, lambda)
{
((alpha*beta)/(lambda))*(x^(alpha-1))*(1+(1/lambda)*x^(alpha))^(-(beta+1))
}
Cumulative Distribution Function
pWeibullGamma <- function(x, alpha, beta, lambda)
{
1-(1+(1/lambda)*x^(alpha))^(-(beta))
}
Hazard Function
hWeibullGamma <- function(x, alpha, beta, lambda)
{
((alpha*beta)/(lambda))*(x^(alpha-1))*(1+(1/lambda)*x^(alpha))^(-(beta+1))/(1+(1/lambda)*x^(alpha))^(-(beta))
}
Survival Function
sWeibullGamma <- function(x,alpha,beta,lambda)
{
(1+(1/lambda)*x^(alpha))^(-(beta))
}
Estimation
paramWG = fitdist(data = dadosp, distr = 'WeibullGamma', start = c(alpha=1.5,beta=1,lambda=1.5), lower= c(0, 0))
summary(paramWG)
Sample:
dadosp = c(240.3,71.9,271.3, 186.3,241,253,287.4,138.3,206.9,176,270.4,73.3,118.9,203.1,139.7,31,269.6,140.2,205.1,133.2,107,354.6,277,27.6,186,260.9,350.4,242.6,292.5, 112.3,242.8,310.7,309.9,53.1,326.5,145.7,271.5, 117.5,264.7,243.9,182,136.7,103.8,188.3,236,419.8,338.6,357.7)
For your sample, the algorithm does not converge when estimating the ML. Fitting a Weibull-Gamma distribution to this data would require an extremely high lambda
value. You can solve this problem by estimating log10(lambda)
instead of lambda
.
You can add lambda <- 10^lambda
inside your 4 functions, e.g.
dWeibullGamma <- function(x, alpha, beta, lambda)
{
lambda <- 10^lambda
((alpha*beta)/(lambda))*(x^(alpha-1))*(1+(1/lambda)*x^(alpha))^(-(beta+1))
}
Then, the algorithm seems to converge:
library(fitdistrplus)
paramWG = fitdist(data = data, distr = 'WeibullGamma',
start = list(alpha=1, beta=1, lambda=1), lower = c(0, 0, 0))
summary(paramWG)$estimate
Output:
alpha beta lambda
2.432939 799.631852 8.680802
We see that the estimate of lambda is 10^8.68
, hence the convergence problem when not taking the log.
You can also have a look at the fit as follows:
newx <- 0:500
pars <- summary(paramWG)$estimate
pred <- dWeibullGamma(newx, pars["alpha"], pars["beta"], pars["lambda"])
hist(data, freq = FALSE)
lines(newx, pred, lwd = 2)
Note: maybe fitting another distribution would make more sense?