Ajuste de una distribución Weibull de 3 parámetros en R

Question

He estado haciendo algunos análisis de datos en R y estoy tratando de encontrar la manera de adaptar mis datos a una distribución Weibull de 3 parámetros. Encontré cómo hacerlo con un Weibull de 2 parámetros pero me quedé corto al encontrar cómo hacerlo con un parámetro de 3.

Aquí es cómo adaptarse a la función() del paquete de datos utilizando el MASA fitdistr:

y <- fitdistr(x[[6]], 'weibull') 

x[[6]] es un subconjunto de los datos de mi ey es donde estoy almacenando el resultado de la conexión.

Answer 1

En primer lugar, es posible que desee consultar FAdist package. Sin embargo, eso no es tan difícil de pasar de rweibull3 a rweibull:

> rweibull3 
function (n, shape, scale = 1, thres = 0) 
thres + rweibull(n, shape, scale) 
<environment: namespace:FAdist> 

y de manera similar de dweibull3 a dweibull

> dweibull3 
function (x, shape, scale = 1, thres = 0, log = FALSE) 
dweibull(x - thres, shape, scale, log) 
<environment: namespace:FAdist> 

así que tenemos esta

> x <- rweibull3(200, shape = 3, scale = 1, thres = 100) 
> fitdistr(x, function(x, shape, scale, thres) 
     dweibull(x-thres, shape, scale), list(shape = 0.1, scale = 1, thres = 0)) 
     shape   scale   thres  
    2.42498383  0.85074556 100.12372297 
( 0.26380861) ( 0.07235804) ( 0.06020083) 

Editar: Como mencionado en el comentario, aparecen varias advertencias al intentar ajustar la distribución ión de esta manera

Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573, : 
    non-finite finite-difference value [3] 
There were 20 warnings (use warnings() to see them) 
Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573, : 
    L-BFGS-B needs finite values of 'fn' 
In dweibull(x, shape, scale, log) : NaNs produced 

para mí al principio fue sólo NaNs produced, y que no es la primera vez que lo veo, así que pensé que no es tan significativo ya que las estimaciones eran buenas. Después de algunas búsquedas, parecía ser un problema bastante popular y no pude encontrar causa ni solución. Una alternativa podría ser usar el paquete stats4 y la función mle(), pero parecía tener algunos problemas también. Pero puedo ofrecer a usar una versión modificada de code por danielmedic que he comprobado varias veces:

thres <- 60 
x <- rweibull(200, 3, 1) + thres 

EPS = sqrt(.Machine$double.eps) # "epsilon" for very small numbers 

llik.weibull <- function(shape, scale, thres, x) 
{ 
    sum(dweibull(x - thres, shape, scale, log=T)) 
} 

thetahat.weibull <- function(x) 
{ 
    if(any(x <= 0)) stop("x values must be positive") 

    toptim <- function(theta) -llik.weibull(theta[1], theta[2], theta[3], x) 

    mu = mean(log(x)) 
    sigma2 = var(log(x)) 
    shape.guess = 1.2/sqrt(sigma2) 
    scale.guess = exp(mu + (0.572/shape.guess)) 
    thres.guess = 1 

    res = nlminb(c(shape.guess, scale.guess, thres.guess), toptim, lower=EPS) 

    c(shape=res$par[1], scale=res$par[2], thres=res$par[3]) 
} 

thetahat.weibull(x) 
    shape  scale  thres 
3.325556 1.021171 59.975470