Inserimento di una distribuzione Weibull a 3 parametri in R

Question

Ho eseguito alcune analisi dei dati in R e sto cercando di capire come adattare i miei dati a una distribuzione Weibull a 3 parametri. Ho trovato come farlo con un Weibull a 2 parametri, ma sono venuto meno nel trovare come farlo con un parametro 3.

Ecco come mi adatto i dati utilizzando la funzione fitdistr() dal pacchetto MASS:

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

x[[6]] è un sottoinsieme dei miei dati e Y è dove mi trovo memorizzare il risultato del raccordo.

Answer 1

In primo luogo, si potrebbe voler guardare FAdist package. Tuttavia, non è così difficile andare rweibull3-rweibull:

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

e simile, da dweibull3 a dweibull

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

così abbiamo questo

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

Modifica: Come menzionato nel commento, appaiono vari avvisi quando si cerca di adattare la distribuzione ione in questo modo

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 

Per me in un primo momento è stato solo NaNs produced, e che non è la prima volta quando la vedo così ho pensato che non è così significativo dal momento che le stime erano buone. Dopo alcune ricerche sembrava essere un problema abbastanza popolare e non riuscivo a trovare né la causa né la soluzione. Un'alternativa potrebbe essere l'utilizzo del pacchetto stats4 e della funzione mle(), ma sembra che abbia anche qualche problema. Ma io posso offrire di utilizzare una versione modificata di code da danielmedic che ho controllato un paio di volte:

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