Warum gibt mir nls () die Fehlermeldung "singuläre Gradientenmatrix bei anfänglichen Parameterschätzungen"?

21

Ich habe einige grundlegende Daten zur Emissionsreduzierung und zu den Kosten pro Auto:

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

Ich weiß, dass dies eine Exponentialfunktion ist, daher erwarte ich, ein Modell finden zu können, das zu Folgendem passt:

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

aber ich erhalte eine Fehlermeldung:

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

Ich habe eine Menge Fragen zu dem Fehler durchgelesen, den ich sehe, und ich stelle fest , dass das Problem wahrscheinlich darin besteht, dass ich bessere / andere startWerte brauche (das initial parameter estimatesmacht ein bisschen mehr Sinn), aber ich bin mir angesichts der Tatsache nicht sicher, ob das so ist Daten, die ich habe, wie ich vorgehen würde, um bessere Parameter abzuschätzen.

Amanda
quelle
Ich würde vorschlagen, dass Sie mit der Entschlüsselung beginnen, indem Sie auf unserer Website nach der Fehlermeldung suchen .
Whuber
3
Tatsächlich habe ich das getan und meine Suche nach dem vollständigen Fehler ergab eine halbherzige Frage mit drei Datenpunkten und keiner Antwort. Ihre genauere Suche liefert jedoch einige Ergebnisse. Möglicherweise, weil Sie hier mehr Erfahrung haben und wissen, welche Begriffe relevant sind.
Amanda
Eine Sache, die ich bei Softwarefehlern festgestellt habe, ist, dass eine Suche nach der spezifischen Fehlermeldung (normalerweise in Anführungszeichen) der sicherste Weg ist, um herauszufinden, ob sie zuvor besprochen wurde. (Dies gilt für das gesamte Internet, nicht nur für SE-Sites.) Wie in unserer Wartemeldung angegeben, sollten Sie sich bei der Behebung Ihres Problems zurückmelden, um uns etwas zurückzudrängen: Diese Frage ist an die Schnittstelle von Statistik und Informatik und könnte hier einige Themen von großem Interesse aufdecken.
whuber
1
Die Anpassung an Ihre Startwerte ist sehr weit von den Daten entfernt. vergleiche exp(50)und exp(95)mit den y-Werten bei x = 50 und x = 95. Wenn Sie c=0y setzen und protokollieren (eine lineare Beziehung erstellen ), können Sie mithilfe der Regression anfängliche Schätzungen für log ( ) und b abrufen, die für Ihre Daten ausreichen (oder wenn Sie eine Linie durch den Ursprung ziehen, können Sie diese verlassen a bei 1 und verwende einfach die Schätzung für b ; das reicht auch für deine Daten). Wenn b außerhalb eines relativ engen Intervalls um diese beiden Werte liegt, treten einige Probleme auf. [Alternativ versuchen Sie es mit einem anderen Algorithmus]einbeinbb
Glen_b -Reinstate Monica
1
Danke @Glen_b. Ich hatte gehofft, ich könnte R anstelle eines Grafikrechners verwenden, um ein Statistik-Intro-Lehrbuch durchzuarbeiten (und den Kurs selbst zu überspringen), also beginne ich nur mit der bloßsten statistischen Einsicht, aber mit viel Erfahrung beim anderen Schneiden und Würfeln in R .
Amanda

Antworten:

37

Gute Startwerte für ein nichtlineares Modell automatisch zu finden, ist eine Kunst. (Für einmalige Datasets ist es relativ einfach, wenn Sie die Daten einfach grafisch darstellen und einige gute Schätzungen vornehmen können.) Ein Ansatz besteht darin, das Modell zu linearisieren und Schätzungen der kleinsten Quadrate zu verwenden.

In diesem Fall hat das Modell die Form

E(Y.)=einexp(bx)+c

für unbekannte Parameter . Das Vorhandensein des Exponentials ermutigt uns, Logarithmen zu verwenden - aber das Hinzufügen von c macht es schwierig, dies zu tun. Beachten Sie jedoch, dass wenn a positiv ist, c kleiner als der kleinste erwartete Wert von Y ist und daher möglicherweise etwas kleiner als der kleinste beobachtete Wert von Y ist . (Wenn a negativ sein könnte, müssen Sie auch einen Wert von c berücksichtigen, der etwas größer als der größte beobachtete Wert von Y ist .)ein,b,cceincY.Y.eincY.

Kümmern wir uns also um indem wir als Anfangsschätzung c 0 etwa die Hälfte des Minimums der Beobachtungen y i verwenden . Das Modell kann jetzt ohne diesen heiklen Zusatzbegriff als umgeschrieben werdencc0yich

E(Y.)-c0einexp(bx).

Dass wir das Protokoll nehmen können von:

Log(E(Y.)-c0)Log(ein)+bx.

Log(ein)b

Hier ist der überarbeitete Code:

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)

Seine Ausgabe (für die Beispieldaten) ist

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

Die Konvergenz sieht gut aus. Lass es uns plotten:

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

Zahl

Es hat gut funktioniert!

yein<0


Eine andere Methode zur Schätzung von Anfangswerten beruht auf dem Verständnis ihrer Bedeutung, die auf Erfahrung, physikalischer Theorie usw. beruhen kann . Ein erweitertes Beispiel für eine (mäßig schwierige) nichtlineare Anpassung, deren Anfangswerte auf diese Weise bestimmt werden können, ist in meiner Antwort beschrieben unter /stats//a/15769 .

Die visuelle Analyse eines Streudiagramms (zur Bestimmung der anfänglichen Parameterschätzungen) wird unter /stats//a/32832 beschrieben und veranschaulicht .

Unter bestimmten Umständen wird eine Folge nichtlinearer Anpassungen erstellt, bei denen Sie erwarten können, dass sich die Lösungen langsam ändern. In diesem Fall ist es oft bequem (und schnell), die vorherigen Lösungen als erste Schätzungen für die nächsten zu verwenden . Ich erinnere mich, diese Technik (ohne Kommentar) unter /stats//a/63169 verwendet zu haben .

whuber
quelle
2

Diese Bibliothek konnte mein Problem mit nls lösen singular gradient: http://www.r-bloggers.com/a-better-nls/ Ein Beispiel:

library(minpack.lm)
nlsLM(function, start=list(variable=2,variable2=12))
Joshring
quelle
Diese Funktion scheint jetzt aufgerufen zu werden nls.lm.
Matt
-1

Also ... ich glaube, ich habe dies als Exponentialfunktion falsch gelesen. Ich brauchte nurpoly()

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

Oder mit 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
quelle
2
Dies beantwortet die von Ihnen gestellte Frage nicht - es ändert sich zu etwas anderem (und weniger interessant, IMHO).
Whuber
6
Obwohl dies das Problem der Anpassung einer Funktion zur Darstellung der Daten lösen kann, sind Ihre akzeptierten Antworten nicht die erwarteten für Ihre Frage. Mr. @whuber hat Ihnen eine hervorragende Erklärung gegeben und die akzeptierte Antwort verdient.
Lourenco