Wie minimiere ich die Restquadratsumme einer Exponentialanpassung?

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.

r nonlinear-regression fitting nls Strohmi
quelle

In diesem Fall , wenn Sie ein Regressionsmodell betrachten

, wobei

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

{Emissions}_{i} = f ({Days}_{i}, a, b) + ϵ_{i}

$\text{Emissions}_i=f(\text{Days}_i,a,b)+\epsilon_i$

ϵ_{i} \sim N (0, σ)

$\epsilon_i\sim N(0,\sigma)$

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

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

@whuber - vielleicht solltest du das als Antwort posten?

Bogenschütze

Ein (negatives) Exponentialgesetz hat die Form . Wenn Sie Änderungen der Einheiten in den und Werten zulassen , sagen wir jedoch und , dann wird das Gesetz ausgedrückt als $y=-\exp(-x)$ $x$ $y$ $y = \alpha y' + \beta$ $x = \gamma x' + \delta$

α y^{'} + β = y = - \exp (- x) = - \exp (- γ x^{'} - δ),

$\alpha y' + \beta = y = -\exp(-x) = -\exp(-\gamma x' - \delta),$

das ist algebraisch äquivalent zu

y^{'} = \frac{- 1}{α} \exp (- γ x^{'} - δ) - β = a (1 - u \exp (- b x^{'}))

$y' = \frac{-1}{\alpha} \exp(-\gamma x' - \delta) - \beta = a\left(1 - u\exp(-b x')\right)$

unter Verwendung von drei Parametern , und . Wir können als Skalierungsparameter für , als Skalierungsparameter für $a = -\beta/\alpha$ $u = 1/(\beta\exp(\delta))$ $b = \gamma$ $a$ $y$ $b$ $x$ und wie aus einem Ableiten Position für den Parameter . $u$ $x$

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 . $a$ $2000$
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 ; relativ sind das ungefähr der Asymptote. $u$ $2000 - 937$ $0.55$
Da , sollte die Kurve , wenn dem dreifachen Wert von auf etwa oder ihrer Gesamtheit angestiegen sein. des Anstiegs von auf fast um ; Das Scannen des Diagramms ergab, dass dies bis Tage dauerte . Nennen wir es der Einfachheit halber, von wo aus $\exp(-3) \approx 0.05$ $x$ $1/b$ $1-0.05$ $95\%$ $95\%$ $937$ $2000$ $1950$ $20$ $25$ $24$ . (Diese Methode zur Ermittlung einer Exponentialskala ist in einigen Bereichen, in denen Exponentialkurven häufig verwendet werden, Standard.) $b \approx 3/24 = 0.125$ $95\%$

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

$2$ $1945$ $1995$

whuber
quelle

res <- residuals(fit); res %*% res

u

$u$

2724

$2724$

147073

$147073$

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

x

$x$

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

Wie minimiere ich die Restquadratsumme einer Exponentialanpassung?

Antworten: