Подгонка кривой распределения Вейбулла в R с помощью NLS

Question

Я пытаюсь соответствовать эти данные в Вейбулла распределения:

Мои x и y переменные:

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) 

Сюжет выглядит следующим образом:

Я ищу что-то вроде этого: fitted plot

и я хочу подгонять к нему кривую Weibull. Я использую функцию NLS в R, как это:

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

Эта функция всегда подбрасывает ошибку говоря:

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 

Так первый я пробовал различные начальные значения без какого-либо успеха. Я не могу понять, как сделать «хорошее» предположение о стартовых значениях. Затем я пошел с функцией SSweibull(x, Asym, Drop, lrc, pwr), которая является функцией самонастройки. Теперь функция SSWeibull ожидает значения для Asym, Drop, lrc и pwr, и у меня нет никаких подсказок относительно того, какими могут быть эти значения.

Буду признателен, если кто-то поможет мне разобраться, как действовать.

Фон данных: Я взял некоторые данные из bugzilla, а моя переменная «y» - количество ошибок, сообщенных за определенный месяц, а переменная «x» - это номер месяца после выпуска.

Answer 1

Возможно, вы захотите изменить формулу, чтобы лучше соответствовать данным. Например, вы можете добавить перехват (так как ваши данные плоские линии в 1 вместо 0, которые модель хочет сделать) и скалярный множитель, так как вы на самом деле не устанавливаете плотность.

Всегда стоит потратить некоторое время на размышления о том, какие параметры имеют смысл, поскольку процедуры подгонки моделей часто весьма чувствительны к первоначальным оценкам. Вы также можете выполнить поиск по сетке, где вы найдете диапазоны возможных параметров, и попробуйте подгонять модель с помощью различных комбинаций, используя функции ловушки ошибок. nls2 имеет возможность сделать это за вас.

Так, например,

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

Не идеально, но это начало.