Dlaczego nls () podaje mi błędy „pojedynczej macierzy gradientu przy początkowych oszacowaniach parametrów”?

Mam podstawowe dane na temat redukcji emisji i kosztu na samochód:

q24 <- read.table(text = "reductions  cost.per.car
    50  45
    55  55
    60  62
    65  70
    70  80
    75  90
    80  100
    85  200
    90  375
    95  600
    ",header = TRUE, sep = "")

Wiem, że jest to funkcja wykładnicza, dlatego spodziewam się, że uda mi się znaleźć model pasujący do:

    model <- nls(cost.per.car ~ a * exp(b * reductions) + c, 
         data = q24, 
         start = list(a=1, b=1, c=0))

ale pojawia się błąd:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Przeczytałem mnóstwo pytań na temat błędu, który widzę i rozumiem, że problemem jest prawdopodobnie to, że potrzebuję lepszych / różnych startwartości ( initial parameter estimatesma to trochę więcej sensu), ale nie jestem pewien, biorąc pod uwagę dane, które mam, jak poszedłbym na temat szacowania lepszych parametrów.

Automatyczne znajdowanie dobrych wartości początkowych dla modelu nieliniowego jest sztuką. (Jest to stosunkowo łatwe w przypadku jednorazowych zestawów danych, gdy można po prostu wykreślić dane i dokonać pewnych domysłów wizualnie.) Jednym z podejść jest linearyzacja modelu i użycie oszacowań metodą najmniejszych kwadratów.

W takim przypadku model ma formę

dla nieznanych parametrów . Obecność wykładnicza zachęca nas do korzystania z logarytmów - ale dodanie c utrudnia to. Zauważmy jednak, że jeśli jest dodatnie, to C jest mniejsza od najmniejszej wartości oczekiwanej Y --and zatem może być nieco mniejsza niż najmniejszy obserwowanych wartości Y . (Jeśli a może być ujemne, musisz również wziąć pod uwagę wartość c, która jest nieco większa niż największa obserwowana wartość Y ).$a,b,c$$c$$a$$c$$Y$$Y$$a$$c$$Y$

Zajmijmy się zatem , stosując jako wstępne oszacowanie c 0 około połowy minimum obserwacji y i . Model można teraz przepisać bez tego kolczastego dodatku jak$c$$c_0$$y_i$

Że możemy wziąć dziennik:

$\log(a)$$b$

Oto poprawiony kod:

c.0 <- min(q24$cost.per.car) * 0.5
model.0 <- lm(log(cost.per.car - c.0) ~ reductions, data=q24)
start <- list(a=exp(coef(model.0)[1]), b=coef(model.0)[2], c=c.0)
model <- nls(cost.per.car ~ a * exp(b * reductions) + c, data = q24, start = start)

Jego dane wyjściowe (dla przykładowych danych) to

Nonlinear regression model
  model: cost.per.car ~ a * exp(b * reductions) + c
   data: q24
        a         b         c 
 0.003289  0.126805 48.487386 
 residual sum-of-squares: 2243

Number of iterations to convergence: 38 
Achieved convergence tolerance: 1.374e-06

Konwergencja wygląda dobrze. Nakreślmy to:

plot(q24)
p <- coef(model)
curve(p["a"] * exp(p["b"] * x) + p["c"], lwd=2, col="Red", add=TRUE)

Działa dobrze!

$y$$a\lt 0$

Inna metoda szacowania wartości początkowych polega na zrozumieniu, co one oznaczają, co może być oparte na doświadczeniu, teorii fizycznej itp . Rozszerzony przykład (umiarkowanie trudnego) dopasowania nieliniowego, którego wartości początkowe można określić w ten sposób, opisano w mojej odpowiedzi na /stats//a/15769 .

Analiza wizualna wykresu rozrzutu (w celu określenia wstępnych oszacowań parametrów) została opisana i zilustrowana na stronie /stats//a/32832 .

W niektórych okolicznościach powstaje sekwencja nieliniowych dopasowań, w której można oczekiwać, że rozwiązania będą się zmieniać powoli. W takim przypadku często wygodne (i szybkie) jest użycie poprzednich rozwiązań jako wstępnych szacunków dla następnych . Pamiętam, używając tej techniki (bez komentarza) na /stats//a/63169 .

Ta biblioteka była w stanie rozwiązać mój problem z nls singular gradient: http://www.r-bloggers.com/a-better-nls/ Przykład:

library(minpack.lm)
nlsLM(function, start=list(variable=2,variable2=12))

Więc ... Myślę, że źle odczytałem to jako funkcję wykładniczą. Potrzebowałem tylkopoly()

model <- lm(cost.per.car ~ poly(reductions, 3), data=q24)
new.data <- data.frame(reductions = c(91,92,93,94))
predict(model, new.data)

plot(q24)
lines(q24$reductions, predict(model, list(reductions = q24$reductions)))

Lub używając lattice:

xyplot(cost.per.car ~ reductions, data = q24,
       panel = function(x, y) {
         panel.xyplot(x, y)
         panel.lines(x, predict(model,list(reductions = x) ))
       }, 
       xlab = "Reductions", 
       ylab = "Cost per car")