Uzyskiwanie właściwych wartości początkowych dla modelu nls w języku R

Próbuję dopasować prosty model prawa mocy do zestawu danych, który jest następujący:

mydf:

rev     weeks
17906.4 1
5303.72 2
2700.58 3
1696.77 4
947.53  5
362.03  6

Celem jest przepuszczenie linii energetycznej i wykorzystanie jej do przewidywania revwartości na przyszłe tygodnie. Kilka badań doprowadziło mnie do tej nlsfunkcji, którą wdrożyłem w następujący sposób.

newMod <- nls(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))
predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))

Chociaż działa to w przypadku lmmodelu, pojawia się singular gradientbłąd, który, jak rozumiem, ma związek z moimi wartościami początkowymi ai b. Próbowałem różnych wartości, nawet posunąłem się do wykreślenia tego w Excelu, zaliczenia samotnego, uzyskania równania, a następnie użycia wartości z równania, ale wciąż pojawia się błąd. Spojrzałem na pęczek odpowiedzi jak ten i próbował drugą odpowiedź (nie mógł zrozumieć pierwszy), ale bez rezultatu.

Naprawdę przydałaby mi się pomoc w znalezieniu właściwych wartości początkowych. Lub alternatywnie, jakiej innej funkcji mogę użyć zamiast nls.

Jeśli chcesz mydfz łatwością odtworzyć :

mydf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6)) 

Jest to powszechny problem z nieliniowymi modelami najmniejszych kwadratów; jeśli twoje wartości początkowe są bardzo dalekie od optymalnego, algorytm może się nie zbiegać, nawet jeśli może być dobrze zachowany w pobliżu optymalnego.

Jeśli zaczniesz biorąc dzienniki obu stronach i dopasować model liniowy, masz szacunki i jako nachylenie i przecięcia (9.947 i -2.011) (edit: to logarytm naturalny)$\log(a)$$b$

Jeśli używasz tych kierować wartości wyjściowych dla i wszystko wydaje się działać dobrze:$a$$b$

 newMod <- nls(rev ~ a*weeks^b, data=mydf, start = list(a=exp(9.947),b=-2.011))
 predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))
 [1] 17919.2138  5280.7001  2584.0109  1556.1951  1050.1230   761.4947   580.3091   458.6027
 [9]   372.6231   309.4658

Próbować

 newMod <- nls(rev ~ a*weeks^b, data=mydf, startlist(a=17919.2127344,b=-1.76270557120))

Poproszono mnie o rozszerzenie tej odpowiedzi. Ten problem jest tak prosty, że jestem trochę zaskoczony, że nls mu się nie udaje. Prawdziwym problemem jest jednak całe podejście R i filozofia dopasowywania modeli nieliniowych. W świecie rzeczywistym przeskalibyśmy x, aby leżeć w przedziale od -1 do 1, a yy y w przedziale od 0 do 1 (y = ax ^ b). Prawdopodobnie wystarczyłoby to, aby nls się zjednoczył. Oczywiście, jak wskazuje Glen, można dopasować odpowiedni model logarytmiczno-liniowy. Opiera się to na tym, że istnieje prosta transformacja, która linearyzuje model. Często tak nie jest. Problem z procedurami R, takimi jak nls, polega na tym, że nie oferują one pomocy w ponownej parametryzacji modelu. W tym przypadku ponowna parametryzacja jest prosta, wystarczy przeskalować / zaktualizować xiy. Jednak po dopasowaniu modelu użytkownik będzie miał inne parametry aib niż oryginalne. Chociaż łatwo jest obliczyć te pierwotne, drugą trudnością jest to, że uzyskanie oszacowanych odchyleń standardowych dla tych oszacowań parametrów nie jest ogólnie takie proste. Odbywa się to metodą delta, która obejmuje Hesjan prawdopodobieństwa logarytmicznego i niektóre pochodne. Oprogramowanie do szacowania parametrów nieliniowych powinno automatycznie przeprowadzać te obliczenia, aby łatwo można było przeprowadzić ponowną parametryzację modelu. Inną rzeczą, którą oprogramowanie powinno obsługiwać, jest pojęcie faz. Możesz pomyśleć o pierwszym dopasowaniu modelu do wersji Glena jako fazie 1. „Prawdziwy” model pasuje do etapu 2. Inną trudnością jest to, że uzyskanie oszacowanych odchyleń standardowych dla tych oszacowań parametrów nie jest ogólnie takie proste. Odbywa się to metodą delta, która obejmuje Hesjan prawdopodobieństwa logarytmicznego i niektóre pochodne. Oprogramowanie do szacowania parametrów nieliniowych powinno automatycznie przeprowadzać te obliczenia, aby łatwo można było przeprowadzić ponowną parametryzację modelu. Inną rzeczą, którą oprogramowanie powinno obsługiwać, jest pojęcie faz. Możesz pomyśleć o pierwszym dopasowaniu modelu do wersji Glena jako fazie 1. „Prawdziwy” model pasuje do etapu 2. Inną trudnością jest to, że uzyskanie oszacowanych odchyleń standardowych dla tych oszacowań parametrów nie jest ogólnie takie proste. Odbywa się to metodą delta, która obejmuje Hesjan prawdopodobieństwa logarytmicznego i niektóre pochodne. Oprogramowanie do szacowania parametrów nieliniowych powinno automatycznie przeprowadzać te obliczenia, aby łatwo można było przeprowadzić ponowną parametryzację modelu. Inną rzeczą, którą oprogramowanie powinno obsługiwać, jest pojęcie faz. Możesz pomyśleć o pierwszym dopasowaniu modelu do wersji Glena jako fazie 1. „Prawdziwy” model pasuje do etapu 2. Oprogramowanie do szacowania parametrów nieliniowych powinno automatycznie przeprowadzać te obliczenia, aby łatwo można było przeprowadzić ponowną parametryzację modelu. Inną rzeczą, którą oprogramowanie powinno obsługiwać, jest pojęcie faz. Możesz pomyśleć o pierwszym dopasowaniu modelu do wersji Glena jako fazie 1. „Prawdziwy” model pasuje do etapu 2. Oprogramowanie do szacowania parametrów nieliniowych powinno automatycznie przeprowadzać te obliczenia, aby łatwo można było przeprowadzić ponowną parametryzację modelu. Inną rzeczą, którą oprogramowanie powinno obsługiwać, jest pojęcie faz. Możesz pomyśleć o pierwszym dopasowaniu modelu do wersji Glena jako fazie 1. „Prawdziwy” model pasuje do etapu 2.

Dopasowuję twój model do AD Model Builder, który w naturalny sposób obsługuje fazy. W pierwszej fazie oszacowano tylko a. To sprawi, że Twój model trafi na boisko. Szacuje się, że w drugiej fazie a i b uzyskają rozwiązanie. Narzędzie AD Model Builder automatycznie oblicza odchylenia standardowe dla dowolnej funkcji parametrów modelu za pomocą metody delta, aby zachęcić do stabilnej ponownej parametryzacji modelu.

Algorytm Levenberga-Marquardta może pomóc:

modeldf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6))

require(minpack.lm)
fit <- nlsLM(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))

require(broom)
fit_data <- augment(fit)

plot(.fitted~rev, data=fit_data)

Z mojego doświadczenia wynika, że ​​dobrym sposobem na znalezienie wartości początkowych dla parametrów modeli NLR jest zastosowanie algorytmu ewolucyjnego. Z początkowej populacji (100) losowych danych szacunkowych (rodziców) w polu wyszukiwania wybierz najlepsze 20 (potomstwo) i użyj ich, aby pomóc w zdefiniowaniu wyszukiwania w kolejnej populacji. Powtarzaj aż do konwergencji. Nie potrzeba gradientów ani hessianów, tylko oceny SSE. Jeśli nie jesteś zbyt chciwy, to bardzo często działa. Problemy, które często mają ludzie, polegają na tym, że korzystają z wyszukiwania lokalnego (Newton-Raphson) do wykonywania wyszukiwania globalnego. Jak zawsze chodzi o użycie odpowiedniego narzędzia do danego zadania. Bardziej sensowne jest użycie globalnego wyszukiwania EA, aby znaleźć wartości początkowe dla lokalnego wyszukiwania Newton, a następnie pozwolić, aby sprowadzało się to do minimum. Ale, jak w przypadku wszystkich rzeczy, diabeł tkwi w szczegółach.