Anpassung einer 3-Parameter-Weibull-Verteilung in R

Question

Ich habe einige Datenanalysen in R durchgeführt und versuche herauszufinden, wie ich meine Daten an eine 3-Parameter-Weibull-Verteilung anpassen kann. Ich habe herausgefunden, wie man es mit einem 2-Parameter-Weibull macht, bin aber bei der Suche nach einem 3-Parameter zu kurz gekommen.

Hier ist, wie ich die Daten unter Verwendung des fitdistr() -Funktion aus der Masse Paket passen:

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

x[[6]] ist eine Teilmenge meiner Daten und Y ist, wo ich das Ergebnis des Beschlags am speichert.

Answer 1

Zuerst möchten Sie vielleicht FAdist package betrachten. Das ist jedoch nicht so schwer rweibull3-rweibull zu gehen:

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

und in ähnlicher Weise von dweibull3 zu dweibull

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

so haben wir diese

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

Edit: Als Im Kommentar erwähnt, gibt es verschiedene Warnungen beim Versuch, den Verteiler anzupassen Ion auf diese Weise

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 

Für mich anfangs war es nur NaNs produced, und das ist nicht das erste Mal, wenn ich es sehe so dachte ich, dass es nicht so sinnvoll ist, da Schätzungen waren gut. Nach einigem Suchen schien es ein recht beliebtes Problem zu sein und ich konnte weder Ursache noch Lösung finden. Eine Alternative könnte stats4 Paket und mle() Funktion sein, aber es schien auch einige Probleme zu haben. Aber ich kann Sie biete eine modifizierte Version von code von danielmedic zu verwenden, die ich habe ein paar Mal geprüft:

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