établir une courbe de distribution de Weibull en R à l'aide nls

Question

Je suis en train de s'adapter ces données à une distribution de Weibull:

Mes x et y variables sont:

y <- c(1, 1, 1, 4, 7, 20, 7, 14, 19, 15, 18, 3, 4, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1) 
x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24) 

L'intrigue ressemble à ceci:

Je cherche quelque chose comme ceci:

et je veux y insérer une courbe de Weibull. J'utilise la fonction nls dans R comme ceci:

nls(y ~ ((a/b) * ((x/b)^(a-1)) * exp(- (x/b)^a))) 

Cette fonction lance toujours une erreur en disant:

Error in numericDeriv(form[[3L]], names(ind), env) : 
    Missing value or an infinity produced when evaluating the model 
In addition: Warning message: 
In nls(y ~ ((a/b) * ((x/b)^(a - 1)) * exp(-(x/b)^a))) : 
    No starting values specified for some parameters. 
Initializing ‘a’, ‘b’ to '1.'. 
Consider specifying 'start' or using a selfStart model 

Alors d'abord j'ai essayé différentes valeurs de départ, sans succès. Je n'arrive pas à comprendre comment faire une "bonne" estimation des valeurs de départ. Ensuite, je suis allé avec la fonction SSweibull(x, Asym, Drop, lrc, pwr) qui est une fonction de démarrage automatique. Maintenant, la fonction SSWeibull attend des valeurs pour Asym, Drop, lrc et pwr et je n'ai aucune idée de ce que ces valeurs pourraient être.

J'apprécierais que quelqu'un puisse m'aider à comprendre comment procéder. Contexte des données: J'ai pris quelques données de bugzilla et ma variable "y" est le nombre de bogues signalés dans un mois particulier et "x" variable est le numéro du mois après la publication.

Answer 1

Vous pouvez envisager de modifier votre formule pour mieux l'adapter aux données. Par exemple, vous pouvez ajouter une interception (puisque vos données flatlines à 1 au lieu de 0, ce que le modèle veut faire) et un multiplicateur scalaire puisque vous n'effectuez pas de densité.

Il est toujours utile de passer du temps à réfléchir aux paramètres qui ont du sens, car les procédures d'ajustement des modèles sont souvent assez sensibles aux estimations initiales. Vous pouvez également faire une recherche de grille où vous arrivez avec des gammes de paramètres possibles et essayez d'ajuster le modèle avec les différentes combinaisons en utilisant des fonctions de capture d'erreur. nls2 a une option pour le faire pour vous.

Ainsi, par exemple,

## Put the data in a data.frame 
dat <- data.frame(x=x, y=y) 

## Try some possible parameter combinations 
library(nls2) 
pars <- expand.grid(a=seq(0.1, 100, len=10), 
        b=seq(0.1, 100, len=10), 
        c=1, 
        d=seq(10, 100, len=10)) 

## brute-force them 
## note the model has changed slightly 
res <- nls2(y ~ d*((a/b) * ((x/b)^(a-1)) * exp(- (x/b)^a)) + c, data=dat, 
      start=pars, algorithm='brute-force') 

## use the results with another algorithm 
res1 <- nls(y ~ d*((a/b) * ((x/b)^(a-1)) * exp(- (x/b)^a)) + c, data=dat, 
      start=as.list(coef(res))) 

## See how it looks 
plot(dat, col="steelblue", pch=16) 
points(dat$x, predict(res), col="salmon", type="l", lwd=2) 

Pas parfait, mais il est un début.