Ich versuche, eine exponentielle Abklingfunktion an y-Werte anzupassen, die bei hohen x-Werten negativ werden, kann meine nls
Funktion jedoch nicht richtig konfigurieren .
Ziel
Ich interessiere mich für die Steigung der Abklingfunktion (nach einigen Quellen ). Wie ich diese Steigung erhalte, ist nicht wichtig, aber das Modell sollte so gut wie möglich zu meinen Daten passen (dh eine Linearisierung des Problems ist akzeptabel , wenn die Anpassung gut ist; siehe "Linearisierung"). In früheren Arbeiten zu diesem Thema wurde jedoch eine folgende exponentielle Abklingfunktion verwendet ( Artikel mit geschlossenem Zugriff von Stedmon et al., Gleichung 3 ):
Wo S
ist die Steigung, die mich interessiert, K
der Korrekturfaktor, um negative Werte zuzulassen, und a
der Anfangswert für x
(dh Achsenabschnitt).
Ich muss dies in R tun, da ich eine Funktion schreibe, die Rohmessungen von chromophoren gelösten organischen Stoffen (CDOM) in Werte umwandelt, an denen Forscher interessiert sind.
Beispieldaten
Aufgrund der Art der Daten musste ich PasteBin verwenden. Die Beispieldaten finden Sie hier .
Schreiben Sie dt <-
den Code von PasteBin und kopieren Sie ihn in Ihre R-Konsole. Dh
dt <- structure(list(x = ...
Die Daten sehen folgendermaßen aus:
library(ggplot2)
ggplot(dt, aes(x = x, y = y)) + geom_point()
Negative y-Werte treten auf, wenn .
Ich versuche eine Lösung zu finden mit nls
Der erste Versuch, nls
eine zu verwenden, erzeugt eine Singularität, was keine Überraschung sein sollte, da ich nur die Startwerte für Parameter betrachtet habe:
nls(y ~ a * exp(-S * x) + K, data = dt, start = list(a = 0.5, S = 0.1, K = -0.1))
# Error in nlsModel(formula, mf, start, wts) :
# singular gradient matrix at initial parameter estimates
Nach dieser Antwort kann ich versuchen, die Startparameter besser anzupassen, um die nls
Funktion zu unterstützen:
K0 <- min(dt$y)/2
mod0 <- lm(log(y - K0) ~ x, data = dt) # produces NaNs due to the negative values
start <- list(a = exp(coef(mod0)[1]), S = coef(mod0)[2], K = K0)
nls(y ~ a * exp(-S * x) + K, data = dt, start = start)
# Error in nls(y ~ a * exp(-S * x) + K, data = dt, start = start) :
# number of iterations exceeded maximum of 50
Die Funktion scheint nicht in der Lage zu sein, eine Lösung mit der Standardanzahl von Iterationen zu finden. Erhöhen wir die Anzahl der Iterationen:
nls(y ~ a * exp(-S * x) + K, data = dt, start = start, nls.control(maxiter = 1000))
# Error in nls(y ~ a * exp(-S * x) + K, data = dt, start = start, nls.control(maxiter = 1000)) :
# step factor 0.000488281 reduced below 'minFactor' of 0.000976562
Weitere Fehler. Chuck es! Erzwingen wir einfach die Funktion, um uns eine Lösung zu geben:
mod <- nls(y ~ a * exp(-S * x) + K, data = dt, start = start, nls.control(maxiter = 1000, warnOnly = TRUE))
mod.dat <- data.frame(x = dt$x, y = predict(mod, list(wavelength = dt$x)))
ggplot(dt, aes(x = x, y = y)) + geom_point() +
geom_line(data = mod.dat, aes(x = x, y = y), color = "red")
Nun, das war definitiv keine gute Lösung ...
Das Problem linearisieren
Viele Menschen haben ihre exponentiellen Zerfallsfunktionen mit Erfolg linearisiert (Quellen: 1 , 2 , 3 ). In diesem Fall müssen wir sicherstellen, dass kein y-Wert negativ oder 0 ist. Lassen Sie uns den minimalen y-Wert innerhalb der Gleitkomma-Grenzen von Computern so nahe wie möglich an 0 bringen :
K <- abs(min(dt$y))
dt$y <- dt$y + K*(1+10^-15)
fit <- lm(log(y) ~ x, data=dt)
ggplot(dt, aes(x = x, y = y)) + geom_point() +
geom_line(aes(x=x, y=exp(fit$fitted.values)), color = "red")
Viel besser, aber das Modell verfolgt y-Werte bei niedrigen x-Werten nicht perfekt.
Beachten Sie, dass die nls
Funktion immer noch nicht zum exponentiellen Zerfall passt:
K0 <- min(dt$y)/2
mod0 <- lm(log(y - K0) ~ x, data = dt) # produces NaNs due to the negative values
start <- list(a = exp(coef(mod0)[1]), S = coef(mod0)[2], K = K0)
nls(y ~ a * exp(-S * x) + K, data = dt, start = start)
# Error in nlsModel(formula, mf, start, wts) :
# singular gradient matrix at initial parameter estimates
Sind die negativen Werte wichtig?
Die negativen Werte sind offensichtlich ein Messfehler, da die Absorptionskoeffizienten nicht negativ sein können. Was ist, wenn ich die y-Werte großzügig positiv mache? Es ist die Steigung, die mich interessiert. Wenn die Addition die Steigung nicht beeinflusst, sollte ich abgerechnet werden:
dt$y <- dt$y + 0.1
fit <- lm(log(y) ~ x, data=dt)
ggplot(dt, aes(x = x, y = y)) + geom_point() + geom_line(aes(x=x, y=exp(fit$fitted.values)), color = "red")
Nun, das lief nicht so gut ... Hohe x-Werte sollten offensichtlich so nahe wie möglich bei Null liegen.
Die Frage
Ich mache hier offensichtlich etwas falsch. Was ist der genaueste Weg, um die Steigung für eine exponentielle Abklingfunktion zu schätzen, die an Daten mit negativen y-Werten unter Verwendung von R angepasst ist?
quelle
nls
konvergierte für mich mit den Startwertennls(y~SSasymp(x, Asym, r0, lrc), data = dt)
. Das konvergiert auch für mich.Antworten:
Verwenden Sie eine Selbststartfunktion:
Ich würde jedoch ernsthaft darüber nachdenken, ob Ihr Domain-Wissen es nicht rechtfertigt, die Asymptote auf Null zu setzen. Ich glaube, dass dies der Fall ist und das obige Modell nicht anderer Meinung ist (siehe Standardfehler / p-Wert des Koeffizienten).
quelle
SSasymp
Funktion. Vielen Dank! Ich glaube, die Forscher möchten sich auf den Artikel beziehen, den ich in der Frage zitiert habe, und denK
Begriff verwenden, aber ich werde ihnen vorschlagen, ihre Gleichung zu ändern. Ich denke, sie wollen das behaltenK
, weil negative Werte bedeuten, dass sich das Instrument nicht wie erwartet verhalten hat, aber sie sind an der Steigung interessiert. Das Entfernen der negativen Asymptote kann in einigen Fällen die Steigung beeinflussen.Diese Frage hat Beziehungen zu mehreren anderen Fragen
Ich habe drei zusätzliche Bemerkungen zu einigen Punkten in dieser Frage.
1: Warum das linearisierte Modell nicht gut zu den großen Werten von passty
Die linearisierte Anpassung minimiert nicht die gleichen Residuen. Auf der logarithmischen Skala sind die Residuen für kleinere Werte größer. Das Bild unten zeigt den Vergleich durch Auftragen der y-Achse auf einer logarithmischen Skala im rechten Bild:
Bei Bedarf können Sie der Verlustfunktion für kleinste Quadrate Gewichte hinzufügen.
2: Verwenden der linearisierten Anpassung als Startwerte
Nachdem Sie Schätzungen mit Ihrer linearisierten Anpassung erhalten haben, können Sie diese als Ausgangspunkt für die nichtlineare Anpassung verwenden. *
3: Verwenden einer allgemeineren Methode, um den Ausgangspunkt zu erhalten
Wenn Sie genügend Punkte haben, können Sie auch die Steigung erhalten, ohne sich um asymptotische Werte und negative Werte kümmern zu müssen (keine Berechnung eines Logarithmus erforderlich).
Sie können dies tun, indem Sie die Datenpunkte integrieren. Dann mity=aesx+k und Y=asesx+kx+Const Sie können ein lineares Modell verwenden, um den Wert von zu erhalten s durch beschreiben y als lineare Kombination der Vektoren Y , x und ein Abschnitt:
Der Vorteil dieser Methode (siehe Tittelbach und Helmrich 1993 "Eine Integrationsmethode zur Analyse multiexponentieller transienter Signale" ) besteht darin, dass Sie sie auf mehr als eine einzelne exponentiell abfallende Komponente erweitern können (Hinzufügen weiterer Integrale).
Fußnote: * Diese Verwendung der Steigung im linearisierten Problem ist genau das, was die
SSasymp
Selbststartfunktion bewirkt. Es schätzt zuerst die Asymptoteund dann die Steigung durch (Subtrahieren des Asymptotenwerts und Nehmen der logarithmischen Werte)
Beachten Sie die Zeile
start = list(lrc = log(-coef(lm(log(abs(y - rt)) ~ x, data = xy))[[2L]]))
Nebenbemerkung: Im Sonderfall dasK=0 Sie können verwenden
welches den beobachteten Parameter modellierty wie
quelle