3

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)
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
user55546
  • 37
  • 1
  • 15
  • I cannot reproduce the error, could you provide a sample with which you get the error? Also, can you confirm you are using the function `fitdist`from the `fitdistrplus` package? – Alexandre Oct 09 '18 at 13:26
  • Many thanks for the reply. Yes, I can, I'll send it to you !! Yes I am using fitdist from the fitdistrplus package. My database contains the monthly rainfall level of the city of Maringá in the State of Paraná, located in Brazil. I'm just sending you all the January observations. (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) – user55546 Oct 09 '18 at 14:31
  • Breno, read in your data into your session, `data <- read.csv("some_data, etc)`, then `dput(data)`, this will print it all out to your screen, copy everything between `structure(...), and class = "data.frame")`, and paste it above to have a working example of your data. – Chris Oct 09 '18 at 14:52
  • my database has 500+ observations, it does not fit. – user55546 Oct 09 '18 at 20:06
  • You can just provide a subset of your data as you did in you comment. Just edit your post so that it is easier for us to reproduce your problem and find a solution. For example, add something like `library(fitdistrplus) dadosp <- c(240.3, 71.9, 271.3, ...)` in the beginning of your example. See also [the guidelines to provide a reproducible example](https://stackoverflow.com/help/mcve) – Alexandre Oct 09 '18 at 20:22
  • IT'S OK! See an issue – user55546 Oct 10 '18 at 01:15

1 Answers1

4

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)

fit

Note: maybe fitting another distribution would make more sense?

Alexandre
  • 452
  • 2
  • 10
  • Alexandre, could you help me with this other question? https://stackoverflow.com/questions/52731375/modified-weibull-error-statistics-aic-what-could-be-happening A guy helped me on the density function (see the question answers) but now I have the error code 100. What do I do? – user55546 Oct 10 '18 at 17:24
  • It's ok! How can I start my own help? Regarding the new question, I tried to use your tip, but it did not work :( – user55546 Oct 10 '18 at 18:49
  • I understood, I already verified 3 times the functions, I do not know the reason of the error. I'll post the link to the article in which you have a role: https://www.sciencedirect.com/science/article/pii/S0951832012002396 – user55546 Oct 10 '18 at 19:11
  • 1
    Ok, I checked the functions and they are correct. The problem is actually the same because of unusual parameter values. Another dirty solution is for example to add `alpha <- alpha/1e5`, `beta <- beta/1e5`, etc. in the beginning of each function in a similar way. Then I managed to estimate your parameters. Otherwise, you would have to implement yourself an optimization algorithm or use another one. You could also consider in these cases transforming your original data. – Alexandre Oct 10 '18 at 19:48
  • Could you estimate? What initial kick did you use in the fitdist ?? I've set up the functions as you suggested – user55546 Oct 10 '18 at 20:01
  • 1
    Yes, I used `paramYuan = fitdist(data = dadosp, distr = 'nmw', start = list(alpha = 0.05, beta = 1, gama = 1.25, theta = 1, lambda = 1), lower = c(0,0,0,0,0))`. And just to be sure, I divided by `1e5` all five parameters in each function. – Alexandre Oct 10 '18 at 20:04