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

21

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.

Amanda
źródło
Sugeruję rozpoczęcie odszyfrowywania, przeszukując naszą stronę w poszukiwaniu komunikatu o błędzie .
whuber
3
Właściwie to zrobiłem i moje wyszukiwanie pełnego błędu przyniosło na wpół wypalone pytanie z trzema punktami danych i bez odpowiedzi. Ale Twoje bardziej szczegółowe wyszukiwanie przynosi pewne wyniki. Być może dlatego, że masz tutaj więcej doświadczenia i wiesz, które warunki wyróżniają się jako odpowiednie.
Amanda,
Jedną z rzeczy, które znalazłem na temat błędów oprogramowania, jest to, że wyszukiwanie określonego komunikatu o błędzie (zwykle w cudzysłowie) jest najpewniejszym sposobem, aby dowiedzieć się, czy zostało to omówione wcześniej. (Dotyczy to całego Internetu, nie tylko witryn SE). Jak głosi nasz komunikat „wstrzymane”, jeśli dodatkowe badania nie rozwiążą problemu, proszę wróć i odeślij się do nas trochę: to pytanie jest na przecięcie statystyki i informatyki i może ujawnić niektóre kwestie o dużym znaczeniu tutaj.
whuber
1
Dopasowanie wartości początkowych jest bardzo dalekie od danych; porównaj exp(50)i exp(95)do wartości y przy x = 50 i x = 95. Jeśli ustawisz c=0i weź dziennik y (tworząc relację liniową), możesz użyć regresji, aby uzyskać wstępne szacunki dla dziennika ( ) ib, które będą wystarczające dla twoich danych (lub jeśli dopasujesz linię przez początek, możesz wyjść a na 1 i po prostu użyj szacunku dla b ; to również wystarcza dla twoich danych). Jeśli b jest znacznie poza dość wąskim przedziałem wokół tych dwóch wartości, napotkasz pewne problemy. [Alternatywnie spróbuj innego algorytmu]ababb
Glen_b
1
Dzięki @Glen_b. Miałem nadzieję, że mógłbym użyć R zamiast kalkulatora graficznego do pracy z podręcznikiem wprowadzającym statystyki (i przeskoczyć sam kurs), więc zaczynam od najciekawszego wglądu statystycznego, ale z dużym doświadczeniem w robieniu innych krojów i krojenia w R .
Amanda

Odpowiedzi:

38

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ę

mi(Y)=zaexp(bx)+do

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 ).za,b,dodozadoYYzadoY

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 jakdodo0yja

mi(Y)-do0zaexp(bx).

Że możemy wziąć dziennik:

log(mi(Y)-do0)log(za)+bx.

log(za)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)

Postać

Działa dobrze!

yza<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 .

Whuber
źródło
2

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))
joshring
źródło
Ta funkcja wydaje się być nls.lmteraz wywoływana .
Matt
-1

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")
Amanda
źródło
2
To nie odpowiada na zadane pytanie - zmienia to na coś innego (a raczej mniej interesującego, IMHO).
whuber
6
Chociaż może to rozwiązać problem dopasowania funkcji do reprezentowania danych, zaakceptowane odpowiedzi nie są oczekiwaniem na twoje pytanie. Pan @whuber zapewnił ci doskonałe wyjaśnienie i zasługujesz na zaakceptowaną odpowiedź.
Lourenco,