Wie minimiere ich die Restquadratsumme einer Exponentialanpassung?

14

Ich habe folgende Daten und möchte ein negatives exponentielles Wachstumsmodell hinzufügen:

Days <- c( 1,5,12,16,22,27,36,43)
Emissions <- c( 936.76, 1458.68, 1787.23, 1840.04, 1928.97, 1963.63, 1965.37, 1985.71)
plot(Days, Emissions)
fit <- nls(Emissions ~ a* (1-exp(-b*Days)), start = list(a = 2000, b = 0.55))
curve((y = 1882 * (1 - exp(-0.5108*x))), from = 0, to =45, add = T, col = "green", lwd = 4)

Der Code funktioniert und eine Anpassungslinie wird gezeichnet. Die Passform ist jedoch optisch nicht optimal und die verbleibende Quadratsumme scheint ziemlich groß zu sein (147073).

Wie können wir unsere Passform verbessern? Ermöglichen die Daten überhaupt eine bessere Anpassung?

Wir konnten im Internet keine Lösung für diese Herausforderung finden. Jede direkte Hilfe oder Verknüpfung zu anderen Websites / Posts wird sehr geschätzt.

Strohmi
quelle
1
In diesem Fall , wenn Sie ein Regressionsmodell betrachten , wobei ε i ~ N ( 0 , σ ) , dann erhalten Sie Ähnliche Schätzer. Durch Auftragen der Vertrauensbereiche kann beobachtet werden, wie diese Werte in den Vertrauensbereichen enthalten sind. Sie können keine perfekte Anpassung erwarten, wenn Sie die Punkte nicht interpolieren oder ein flexibleres nichtlineares Modell verwenden. Emissionsi=f(Tageich,ein,b)+ϵichϵichN(0,σ)
Ich habe den Titel geändert, weil "negatives Exponentialmodell" etwas anderes bedeutet als in der Frage beschrieben.
Whuber
Vielen Dank für die Klarstellung der Frage (@whuber) und für Ihre Antwort (@Procrastinator). Wie kann ich die Vertrauensbereiche berechnen und zeichnen? Und was wäre ein flexibleres nichtlineares Modell?
Strohmi
4
Sie benötigen einen zusätzlichen Parameter. Sehen Sie, was passiert mit fit <- nls(Emissions ~ a* (1- u*exp(-b*Days)), start = list(a = 2000, b = 0.1, u=.5)); beta <- coefficients(fit); curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T).
Whuber
1
@whuber - vielleicht solltest du das als Antwort posten?
Bogenschütze

Antworten:

16

Ein (negatives) Exponentialgesetz hat die Form . Wenn Sie Änderungen der Einheiten in den x- und y- Werten zulassen , sagen wir jedoch y = α y ' + β und x = γ x ' + δ , dann wird das Gesetz ausgedrückt alsy=exp(x)xyy=αy+βx=γx+δ

αy+β=y=exp(x)=exp(γxδ),

das ist algebraisch äquivalent zu

y=1αexp(γxδ)β=a(1uexp(bx))

unter Verwendung von drei Parametern , u = 1 / ( β exp ( δ ) ) und b = γ . Wir können a als Skalierungsparameter für y , b als Skalierungsparameter für x erkennena=β/αu=1/(βexp(δ))b=γaybx und wie aus einem Ableiten Position für den Parameter x .ux

Als Faustregel können diese Parameter auf einen Blick aus der Grafik identifiziert werden :

  • Der Parameter ist der Wert der horizontalen Asymptote, etwas weniger als 2000 .a2000

  • Der Parameter ist der relative Betrag, um den die Kurve vom Ursprung zu ihrer horizontalen Asymptote ansteigt. Hier beträgt der Anstieg also etwas weniger als 2000 - 937 ; relativ sind das ungefähr 0,55 der Asymptote.u20009370.55

  • Da , sollte die Kurve , wenn x dem dreifachen Wert von 1 / b entspricht, auf etwa 1 - 0,05 oder 95 % ihrer Gesamtheit angestiegen sein. 95 % des Anstiegs von 937 auf fast 2000 liegen um 1950 ; Das Scannen des Diagramms ergab, dass dies 20 bis 25 Tage dauerte . Nennen wir es 24 der Einfachheit halber, von wo aus b 3 / 24exp(3)0.05x1/b10.0595%95%93720001950202524 . (Diese 95- prozentige Methode zur Ermittlung einer Exponentialskala ist in einigen Bereichen, in denen Exponentialkurven häufig verwendet werden, Standard.)b3/24=0.12595%

Mal sehen, wie das aussieht:

plot(Days, Emissions)
curve((y = 2000 * (1 - 0.56 * exp(-0.125*x))), add = T)

Eyeball fit

Nicht schlecht für den Anfang! (Trotz des Tippens 0.56anstelle von 0.55, was sowieso eine grobe Annäherung war.) Wir können es polieren mit nls:

fit <- nls(Emissions ~ a * (1- u * exp(-b*Days)), start=list(a=2000, b=1/8, u=0.55))
beta <- coefficients(fit)
plot(Days, Emissions)
curve((y = beta["a"] * (1 - beta["u"] * exp(-beta["b"]*x))), add = T, col="Green", lwd=2)

NLS fit

Die Ausgabe von nlsenthält umfangreiche Informationen zur Parameterunsicherheit. Zum Beispiel kann eine einfache summaryliefert Standardfehler der Schätzungen:

> summary(fit)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a 1.969e+03  1.317e+01  149.51 2.54e-10 ***
b 1.603e-01  1.022e-02   15.69 1.91e-05 ***
u 6.091e-01  1.613e-02   37.75 2.46e-07 ***

Wir können die gesamte Kovarianzmatrix der Schätzungen lesen und damit arbeiten. Dies ist nützlich für die Schätzung simultaner Konfidenzintervalle (zumindest für große Datensätze):

> vcov(fit)
             a             b             u
a 173.38613624 -8.720531e-02 -2.602935e-02
b  -0.08720531  1.044004e-04  9.442374e-05
u  -0.02602935  9.442374e-05  2.603217e-04

nls unterstützt Profildiagramme für die Parameter und gibt detailliertere Informationen zu deren Unsicherheit:

> plot(profile(fit))

Hier ist einer der drei Ausgabediagramme, die die Variation von a :

Profile plot

219451995

whuber
quelle
res <- residuals(fit); res %*% resu2724147073
Alles schön und gut whuber. Aber vielleicht hatte das OP einen Grund, sich für das Exponentialmodell zu entscheiden (oder vielleicht nur, weil es bekannt ist). Ich denke, zuerst sollten die Residuen für das Exponentialmodell betrachtet werden. Plotten Sie sie gegen potenzielle Kovariaten, um zu sehen, ob dort Struktur vorhanden ist und nicht nur großes zufälliges Rauschen. Bevor Sie sich mit anspruchsvolleren Modellen beschäftigen, versuchen Sie, herauszufinden, ob ein schickeres Modell möglicherweise Abhilfe schaffen könnte.
Michael R. Chernick
3
x
2
Ich habe deine Antwort nicht kritisiert! Ich habe keine Restparzellen gesehen. Alles, was ich vorgeschlagen habe, ist, dass Diagramme von Residuen gegen potenzielle Kovariaten der erste Schritt sein sollten, um ein besseres Modell zu finden. Wenn ich geglaubt hätte, eine Antwort zu haben, hätte ich eine Antwort gegeben, anstatt meinen Punkt als Konstante zu erwähnen. Ich dachte, du hast eine großartige Antwort gegeben und ich war einer derjenigen, die dir +1 gegeben haben.
Michael R. Chernick