Erhalten der richtigen Startwerte für ein nls-Modell in R.

12

Ich versuche, ein einfaches Potenzgesetzmodell an einen Datensatz anzupassen, der wie folgt lautet:

mydf::

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

Das Ziel ist es, die Stromleitung durchzuleiten und damit revWerte für zukünftige Wochen vorherzusagen . Eine Reihe von Recherchen hat mich zu der nlsFunktion geführt, die ich wie folgt implementiert habe.

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

Während dies für ein lmModell funktioniert , erhalte ich einen singular gradientFehler, der meines Wissens mit meinen Startwerten aund zu tun hat b. Ich habe verschiedene Werte ausprobiert und bin sogar so weit gegangen, dies in Excel zu zeichnen, einen einzelnen zu übergeben, eine Gleichung zu erhalten und dann die Werte aus der Gleichung zu verwenden, aber ich erhalte immer noch den Fehler. Ich sah mir eine Reihe von Antworten wie diese an und versuchte die zweite Antwort (konnte die erste nicht verstehen), aber ohne Ergebnis.

Ich könnte hier wirklich Hilfe gebrauchen, um die richtigen Startwerte zu finden. Oder alternativ, welche andere Funktion kann ich anstelle von nls verwenden.

Für den Fall, dass Sie mydfmit Leichtigkeit neu erstellen möchten :

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)) 
NeonBlueHair
quelle
1
Obwohl in R angegeben (es muss wirklich in einer Sprache angegeben werden), ist es statistisch genug, geeignete Startwerte für eine nichtlineare Modellanpassung zu finden, um hier zum Thema zu kommen, IMO. Es ist nicht wirklich ein Programmier-Q, z.
Gung - Reinstate Monica

Antworten:

13

Dies ist ein häufiges Problem bei nichtlinearen Modellen der kleinsten Quadrate. Wenn Ihre Startwerte sehr weit vom Optimum entfernt sind, konvergiert der Algorithmus möglicherweise nicht, obwohl er sich in der Nähe des Optimums gut verhält.

Log(ein)b

einb

 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
Glen_b - Monica neu starten
quelle
Das ist sehr hilfreich, vielen Dank! Ich habe allerdings eine Frage, wie Sie hier zu Ihrem "a" -Wert gekommen sind. Ich habe versucht, lm (log10 (rev) ~ log10 (Wochen)) auszuführen und dann die Funktion "summary" zu verwenden. Während ich den gleichen "b" -Wert erhalte, wird mein "a" -Wert auf 4.3201 ausgegeben. Was haben Sie anders gemacht, um zu a = 9.947 zu gelangen?
NeonBlueHair
exploge
Ah, du hast vollkommen recht. Amateurfehler meinerseits. Denken Sie an die mathematische Notation und erwarten Sie, dass "log" log log 10 und "ln" für natürliches log bedeutet. Danke für die Klarstellung.
NeonBlueHair
1
Für viele Mathematiker (und viele Statistiker) ist ein schmuckloses "Protokoll" das natürliche Protokoll, ähnlich wie ein schmuckloses Argument für eine Sündenfunktion im Bogenmaß. [Kollisionskonventionen können leider zu Verwirrung führen, aber als ich zum Beispiel anfing, R zu verwenden, habe ich nicht zweimal über die Verwendung der Protokollfunktion
nachgedacht,
4

Versuchen

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

Ich wurde gebeten, diese Antwort etwas zu erweitern. Dieses Problem ist so einfach, dass ich überrascht bin, dass nls daran scheitert. Das eigentliche Problem ist jedoch der gesamte R-Ansatz und die Philosophie der nichtlinearen Modellanpassung. In der realen Welt würde man x so skalieren, dass es zwischen -1 und 1 liegt und y und y zwischen 0 und 1 liegen (y = ax ^ b). Das würde wahrscheinlich ausreichen, um nls zur Konvergenz zu bringen. Natürlich können Sie, wie Glen betont, das entsprechende logarithmische lineare Modell anpassen. Dies beruht auf der Tatsache, dass es eine einfache Transformation gibt, die das Modell linearisiert. Das ist oft nicht der Fall. Das Problem bei R-Routinen wie nls ist, dass sie keine Unterstützung für die Neuparametrisierung des Modells bieten. In diesem Fall ist die Neuparametrisierung einfach. Skalieren Sie einfach x / y neu. Wenn der Benutzer jedoch zum Modell passt, hat er andere Parameter a und b als die ursprünglichen. Während es einfach ist, die ursprünglichen daraus zu berechnen, besteht die andere Schwierigkeit darin, dass es im Allgemeinen nicht so einfach ist, die geschätzten Standardabweichungen für diese Parameterschätzungen zu erhalten. Dies erfolgt nach der Delta-Methode, bei der der Hessische Wert der logarithmischen Wahrscheinlichkeit und einige Derivate berücksichtigt werden. Nichtlineare Parameterschätzungssoftware sollte diese Berechnungen automatisch bereitstellen, damit die Neuparametrisierung des Modells leicht unterstützt wird. Eine andere Sache, die Software unterstützen sollte, ist der Begriff der Phasen. Sie können sich vorstellen, das Modell zuerst mit Glen's Version als Phase 1 auszurüsten. Das "echte" Modell wird in Stufe 2 angepasst. Die andere Schwierigkeit besteht darin, dass es im Allgemeinen nicht so einfach ist, die geschätzten Standardabweichungen für diese Parameterschätzungen zu erhalten. Dies erfolgt nach der Delta-Methode, bei der der Hessische Wert der logarithmischen Wahrscheinlichkeit und einige Derivate berücksichtigt werden. Nichtlineare Parameterschätzungssoftware sollte diese Berechnungen automatisch bereitstellen, damit die Neuparametrisierung des Modells leicht unterstützt wird. Eine andere Sache, die Software unterstützen sollte, ist der Begriff der Phasen. Sie können sich vorstellen, das Modell zuerst mit Glen's Version als Phase 1 auszurüsten. Das "echte" Modell wird in Stufe 2 angepasst. Die andere Schwierigkeit besteht darin, dass es im Allgemeinen nicht so einfach ist, die geschätzten Standardabweichungen für diese Parameterschätzungen zu erhalten. Dies erfolgt nach der Delta-Methode, bei der der Hessische Wert der logarithmischen Wahrscheinlichkeit und einige Derivate berücksichtigt werden. Nichtlineare Parameterschätzungssoftware sollte diese Berechnungen automatisch bereitstellen, damit die Neuparametrisierung des Modells leicht unterstützt wird. Eine andere Sache, die Software unterstützen sollte, ist der Begriff der Phasen. Sie können sich vorstellen, das Modell zuerst mit Glen's Version als Phase 1 auszurüsten. Das "echte" Modell wird in Stufe 2 angepasst. Nichtlineare Parameterschätzungssoftware sollte diese Berechnungen automatisch bereitstellen, damit die Neuparametrisierung des Modells leicht unterstützt wird. Eine andere Sache, die Software unterstützen sollte, ist der Begriff der Phasen. Sie können sich vorstellen, das Modell zuerst mit Glen's Version als Phase 1 auszurüsten. Das "echte" Modell wird in Stufe 2 angepasst. Nichtlineare Parameterschätzungssoftware sollte diese Berechnungen automatisch bereitstellen, damit die Neuparametrisierung des Modells leicht unterstützt wird. Eine andere Sache, die Software unterstützen sollte, ist der Begriff der Phasen. Sie können sich vorstellen, das Modell zuerst mit Glen's Version als Phase 1 auszurüsten. Das "echte" Modell wird in Stufe 2 angepasst.

Ich passe Ihr Modell mit AD Model Builder an, der Phasen auf natürliche Weise unterstützt. In der ersten Phase wurde nur a geschätzt. Dies bringt Ihr Modell in den Ballpark. In der zweiten Phase werden a und b geschätzt, um die Lösung zu erhalten. AD Model Builder berechnet automatisch die Standardabweichungen für jede Funktion der Modellparameter über die Delta-Methode, um eine stabile Neuparametrisierung des Modells zu fördern.

Dave Fournier
quelle
1

Der Levenberg-Marquardt-Algorithmus kann helfen:

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)
Etienne Low-Décarie
quelle
1

Nach meiner Erfahrung ist die Verwendung eines evolutionären Algorithmus eine gute Möglichkeit, Startwerte für Parameter von NLR-Modellen zu finden. Wählen Sie aus einer anfänglichen Population (100) zufälliger Schätzungen (Eltern) in einem Suchraum die besten 20 (Nachkommen) aus und verwenden Sie diese, um eine Suche in einer nachfolgenden Population zu definieren. Wiederholen bis zur Konvergenz. Keine Notwendigkeit für Gradienten oder Hessen, nur SSE-Bewertungen. Wenn Sie nicht zu gierig sind, funktioniert dies sehr oft. Die Probleme, die Menschen häufig haben, bestehen darin, dass sie eine lokale Suche (Newton-Raphson) verwenden, um die Arbeit einer globalen Suche auszuführen. Wie immer geht es darum, das richtige Werkzeug für den jeweiligen Auftrag zu verwenden. Es ist sinnvoller, eine globale EA-Suche zu verwenden, um Startwerte für die lokale Newton-Suche zu finden, und diese dann auf das Minimum zu reduzieren. Aber wie bei allen Dingen steckt der Teufel im Detail.

Adrian O'Connor
quelle