Strategie zur Anpassung einer stark nichtlinearen Funktion

Zur Analyse von Daten aus einem biophysikalischen Experiment versuche ich derzeit, eine Kurvenanpassung mit einem stark nichtlinearen Modell durchzuführen. Die Modellfunktion sieht grundsätzlich so aus:

$y = ax + bx^{-1/2}$

Hier ist insbesondere der Wert von $b$ von großem Interesse.

Ein Plot für diese Funktion:

Funktionsdiagramm

(Beachten Sie, dass die Modellfunktion auf einer gründlichen mathematischen Beschreibung des Systems basiert und sehr gut zu funktionieren scheint. Es ist nur so, dass automatisierte Anpassungen schwierig sind.)

Natürlich ist die Modellfunktion problematisch: Anpassungsstrategien, die ich bisher ausprobiert habe, scheitern an der scharfen Asymptote bei $x=0$ , insbesondere bei verrauschten Daten.

Mein Verständnis des Problems hier ist, dass die einfache Anpassung der kleinsten Quadrate (ich habe in MATLAB sowohl mit linearer als auch mit nichtlinearer Regression gespielt; meistens Levenberg-Marquardt) für die vertikale Asymptote sehr empfindlich ist, da kleine Fehler in x enorm verstärkt werden .

Könnte mich jemand auf eine passende Strategie hinweisen, die das umgehen könnte?

Ich habe einige Grundkenntnisse in Statistik, aber das ist immer noch ziemlich begrenzt. Ich würde gerne lernen, wenn ich nur wüsste, wo ich anfangen soll :)

Vielen Dank für Ihren Rat!

Bearbeiten Bitte um Verzeihung, dass Sie vergessen haben, die Fehler zu erwähnen. Das einzige signifikante Rauschen ist in $x$ und es ist additiv.

Bearbeiten 2 Einige zusätzliche Informationen zum Hintergrund dieser Frage. Die obige Grafik modelliert das Streckverhalten eines Polymers. Wie @whuber in den Kommentaren ausführte, benötigen Sie $b \approx -200 a$ , um ein Diagramm wie oben zu erhalten.

Wie die Leute diese Kurve bis zu diesem Punkt angepasst haben: Es scheint, dass die Leute im Allgemeinen die vertikale Asymptote abschneiden, bis sie eine gute Passform finden. Die Auswahl der Abschaltung ist jedoch immer noch willkürlich, wodurch das Anpassungsverfahren unzuverlässig und nicht reproduzierbar wird.

3 & 4 Fixed Graph bearbeiten .

curve-fitting nonlinear onnodb
quelle

Treten die Fehler in

oder in auf ?

x

$x$

oder in beiden auf? In welcher Form soll das Rauschen eintreten (multiplikativ, additiv usw.)?

y

$y$

Wahrscheinlichkeit

@onnodb: Mein Anliegen ist, könnte dies nicht grundsätzlich in Frage stellen, wie robust Ihr Modell selbst ist? Egal, welche Anpassstrategie Sie verwenden, bleibt

nicht hochsensibel? Können Sie jemals ein hohes Vertrauen in eine solche Schätzung für

b

$b$

b

$b$ ?

curious_cat

Das geht leider immer noch nicht. Es gibt einfach keine mögliche Kombination von

und

, die den von Ihnen gezeichneten Graphen qualitativ reproduziert. (Offensichtlich ist

negativ.

muss kleiner sein als die geringste Steigung im Graphen, aber positiv, was es in ein enges Intervall bringt. Aber wenn

in diesem Intervall ist, ist es einfach nicht groß genug, um die große negative Spitze bei zu überwinden der Ursprung durch die eingeführten

term.) Was haben Sie gezogen? Daten? Irgendeine andere Funktion?

a

$a$

b

$b$

b

$b$

a

$a$

a

$a$

b x^{1 / 2}

$b x^{1/2}$

whuber

Danke, aber es ist immer noch falsch. Verlängerung der Tangente zu diesem Graphen rückwärts von jedem Punkt

(x, a x + b x^{1 / 2})

$(x,ax+bx^{1/2})$ , wo

, werden Sie die y-Achse bei abfängt

x > 0

$x\gt 0$

. Weil die Abwärtsspitze bei

zeigt

(0, 3 b / (2 x^{1 / 2}))

$(0,3b/(2x^{1/2}))$

0

$0$

b

$b$ negativ ist, muss auch dieser y-Achsenabschnitt negativ sein. In Ihrer Abbildung ist jedoch sehr deutlich, dass die meisten dieser Abschnitte positiv sind und sich auf bis zu

. Somit ist es mathematisch unmöglich , daß eine Gleichung wie kann die programmierte Kurve beschreiben , auch nicht annähernd. Zumindest müssen Sie so etwas wie passen

15.5

$15.5$ $y=ax + bx^{1/2}$

y = a x + b x^{1 / 2} + c

$y=ax+bx^{1/2}+c$

whuber

Bevor ich daran gearbeitet habe, wollte ich die Aussage der Frage überprüfen: Deshalb ist es wichtig, dass die Funktion korrekt ist. Ich habe jetzt keine Zeit, eine vollständige Antwort zu geben, möchte aber bemerken, dass "andere Leute" vielleicht falsch liegen - aber es kommt auf noch mehr Details an, leider. Wenn Ihr

Fehler wirklich additiv ist, muss er meines Erachtens immer noch stark heteroskedastisch sein, da sonst seine Varianz bei kleinen Werten von

wirklich klein wäre. Was können Sie uns quantitativ über diesen Fehler sagen?

x

$x$

x

$x$

Whuber

Antworten:

Die Methoden, die wir verwenden würden, um dies manuell anzupassen (d. H. Von Exploratory Data Analysis), können mit solchen Daten bemerkenswert gut funktionieren.

Ich möchte das Modell leicht umparametrieren , um die Parameter positiv zu machen:

y = a x - b / \sqrt{x} .

$y = a x - b / \sqrt{x}.$

wir für ein gegebenes , dass es ein eindeutiges reales , das diese Gleichung erfüllt. nenne dies oder der Kürze halber wenn $y$ $x$ $f(y; a,b)$ $f(y)$ verstanden wird. $(a,b)$

Wir beobachten eine Ansammlung geordneter Paare denen von abweicht $(x_i, y_i)$ $x_i$ durch unabhängige Zufallsvariablen mit dem Mittelwert Null abweichen. In dieser Diskussion gehe ich davon aus, dass sie alle eine gemeinsame Varianz haben, aber eine Erweiterung dieser Ergebnisse (unter Verwendung gewichteter kleinster Quadrate) ist möglich, offensichtlich und einfach zu implementieren. Hier ist ein simuliertes Beispiel einer solchen Sammlung von Werten mit , und einer gemeinsamen Varianz von $f(y_i; a,b)$ $100$ $a=0.0001$ $b=0.1$ . $\sigma^2=4$

Datenplot

Dies ist ein (absichtlich) schwieriges Beispiel, wie die nichtphysikalischen (negativen) Werte und ihre außergewöhnliche Streuung (die normalerweise horizontale Einheiten beträgt , aber bis zu betragen kann) erkennen lassen $x$ $\pm 2$ $5$ oder auf der Achse ) erkennen lassen. Wenn wir eine vernünftige Übereinstimmung mit diesen Daten erzielen können, die der Schätzung von , und nahekommt , dann sind wir in der Tat erfolgreich. $6$ $x$ $a$ $b$ $\sigma^2$

Eine explorative Anpassung ist iterativ. Jede Stufe besteht aus zwei Schritten: Schätzung (basierend auf den Daten und früheren Schätzungen und von und , von der vorherigen vorhergesagten Werten kann für die erhalten wird , ) , und dann schätzen . Da die Fehler in x sind , schätzen die Anpassungen aus und nicht umgekehrt. Um zuerst die Fehler in , wenn $a$ $\hat{a}$ $\hat{b}$ $a$ $b$ $\hat{x}_i$ $x_i$ $b$ $x_i$ $(y_i)$ $x$ $x$ ausreichend groß ist,

x_{i} \approx \frac{1}{a} (y_{i} + \frac{\hat{b}}{\sqrt{{\hat{x}}_{i}}}) .

$x_i \approx \frac{1}{a}\left(y_i + \frac{\hat{b}}{\sqrt{\hat{x}_i}}\right).$

Daher können wir aktualisieren , durch den Einbau dieses Modell mit der kleinsten Quadrate (Anmerkung es nur einen Parameter aufweist - eine Steigung, --Und kein intercept) und der Kehrwert des Koeffizienten als die aktualisierte Schätzung des Nehmens . $\hat{a}$ $a$ $a$

Als nächstes dominiert , wenn ausreichend klein ist, der umgekehrte quadratische Term und wir finden (wieder in erster Ordnung in den Fehlern), dass $x$

x_{i} \approx b^{2} \frac{1 - 2 \hat{a} \hat{b} {\hat{x}}^{3 / 2}}{y_{i}^{2}} .

$x_i \approx b^2\frac{1 - 2 \hat{a} \hat{b} \hat{x}^{3/2}}{y_i^2}.$

Erneut unter Verwendung der kleinsten Quadrate (mit nur einer Steigung Begriff ) haben wir eine aktualisierte Schätzung erhalten über die Quadratwurzel der angepassten Steigung. $b$ $\hat{b}$

Um zu sehen, warum dies funktioniert, kann eine grobe explorative Annäherung an diese Anpassung erhalten werden, indem gegen für das kleinere . Noch besser wäre es, weil die mit Fehler gemessen werden , und die monoton ändern sich mit der , sollten wir uns auf die Daten konzentrieren mit den größeren Werten von . Hier ist ein Beispiel aus unserem simulierten Datensatz, der die größte Hälfte von $x_i$ $1/y_i^2$ $x_i$ $x_i$ $y_i$ $x_i$ $1/y_i^2$ $y_i$ in rot, die kleinste hälfte in blau und eine linie durch den ursprung passen zu den roten punkten.

Zahl

Die Punkte richten sich ungefähr aus, obwohl es bei den kleinen Werten von und eine leichte Krümmung gibt . (Beachten Sie die Auswahl der Achsen: weil $x$ $y$ $x$ das Maß ist, ist es üblich, es auf der vertikalen Achse zu zeichnen .) Durch Fokussieren der Anpassung auf die roten Punkte, bei denen die Krümmung minimal sein sollte, sollten wir eine vernünftige Schätzung von . Der im Titel angezeigte Wert von ist die Quadratwurzel der Steigung dieser Linie: Es sind nur % weniger als der wahre Wert! $b$ $0.096$ $4$

Zu diesem Zeitpunkt können die vorhergesagten Werte über aktualisiert werden

{\hat{x}}_{i} = f (y_{i}; \hat{a}, \hat{b}) .

$\hat{x}_i = f(y_i; \hat{a}, \hat{b}).$

Iterieren Sie, bis sich die Schätzungen entweder stabilisieren (was nicht garantiert ist) oder durch kleine Wertebereiche laufen (was immer noch nicht garantiert werden kann).

Es stellt sich heraus, dass nur schwer abzuschätzen ist, wenn wir eine gute Menge sehr großer Werte von , aber dass - das die vertikale Asymptote in der ursprünglichen Darstellung (in der Frage) bestimmt und im Mittelpunkt der Frage steht - kann ziemlich genau festgehalten werden, vorausgesetzt, es gibt einige Daten innerhalb der vertikalen Asymptote. In unserem laufenden Beispiel tun die Iterationen konvergieren zu (die fast zweimal der richtige Wert von ist ) und $a$ $x$ $b$ $\hat{a} = 0.000196$ $0.0001$ (was nahe dem korrekten Wert von ist $\hat{b} = 0.1073$ $0.1$ ). Dieses Diagramm zeigt noch einmal die Daten, die überlagert sind (a) die wahre Kurve in grau (gestrichelt) und (b) die geschätzte Kurve in rot (durchgehend):

Passt

Diese Anpassung ist so gut, dass es schwierig ist, die wahre Kurve von der angepassten Kurve zu unterscheiden: Sie überlappen sich fast überall. Im Übrigen liegt die geschätzte Fehlervarianz von sehr nahe am wahren Wert von . $3.73$ $4$

Bei diesem Ansatz gibt es einige Probleme:

Die Schätzungen sind voreingenommen. Die Abweichung wird deutlich, wenn der Datensatz klein ist und relativ wenige Werte nahe an der x-Achse liegen. Die Passform ist systematisch etwas niedrig.
Das Schätzverfahren erfordert ein Verfahren, um "große" von "kleinen" Werten von . Ich könnte explorative Wege vorschlagen, um optimale Definitionen zu identifizieren, aber aus praktischen Gründen können Sie diese als "Tuning" -Konstanten belassen und sie ändern, um die Empfindlichkeit der Ergebnisse zu überprüfen. Ich habe sie willkürlich festgelegt, indem ich die Daten gemäß dem Wert von in drei gleiche Gruppen $y_i$ $y_i$ und die beiden äußeren Gruppen verwendete.
Die Prozedur funktioniert nicht für alle möglichen Kombinationen von und oder alle möglichen Datenbereiche. Es sollte jedoch immer dann gut funktionieren, wenn im Datensatz genügend Kurvenmaterial vorhanden ist, um beide Asymptoten wiederzugeben: die vertikale an einem Ende und die geneigte am anderen Ende. $a$ $b$

Code

Das Folgende ist in Mathematica geschrieben .

estimate[{a_, b_, xHat_}, {x_, y_}] := 
  Module[{n = Length[x], k0, k1, yLarge, xLarge, xHatLarge, ySmall, 
    xSmall, xHatSmall, a1, b1, xHat1, u, fr},
   fr[y_, {a_, b_}] := Root[-b^2 + y^2 #1 - 2 a y #1^2 + a^2 #1^3 &, 1];
   k0 = Floor[1 n/3]; k1 = Ceiling[2 n/3];(* The tuning constants *)
   yLarge = y[[k1 + 1 ;;]]; xLarge = x[[k1 + 1 ;;]]; xHatLarge = xHat[[k1 + 1 ;;]];
   ySmall = y[[;; k0]]; xSmall = x[[;; k0]]; xHatSmall = xHat[[;; k0]];
   a1 = 1/
     Last[LinearModelFit[{yLarge + b/Sqrt[xHatLarge], 
          xLarge}\[Transpose], u, u]["BestFitParameters"]];
   b1 = Sqrt[
     Last[LinearModelFit[{(1 - 2 a1 b  xHatSmall^(3/2)) / ySmall^2, 
          xSmall}\[Transpose], u, u]["BestFitParameters"]]];
   xHat1 = fr[#, {a1, b1}] & /@ y;
   {a1, b1, xHat1}
   ];

Wenden Sie dies auf Daten an (gegeben durch parallele Vektoren xund ygebildet in eine zweispaltige Matrix data = {x,y}) bis zur Konvergenz, beginnend mit Schätzungen von : $a=b=0$

{a, b, xHat} = NestWhile[estimate[##, data] &, {0, 0, data[[1]]}, 
                Norm[Most[#1] - Most[#2]] >= 0.001 &,  2, 100]

whuber
quelle

Dies ist eine erstaunliche Antwort. Ich bin sehr verpflichtet! Ich habe damit gespielt und die Ergebnisse sehen sehr vielversprechend aus. Ich brauche jedoch etwas mehr Zeit, um die Argumentation vollständig zu verstehen :) Außerdem: Könnte ich Sie über Ihre Website wegen einer weiteren (privaten) Frage zu Danksagungen kontaktieren?

Onnodb

Lesen Sie die wichtigen Fragen unter @probabilityislogic

Wenn Sie nur Fehler in y haben und diese additiv sind und Sie eine konstante Varianz haben (dh Ihre Annahmen stimmen mit dem überein, wie Sie sich angehört haben), dann lassen Sie $y^* = y\sqrt{x}$ $y^*$ $x^* = x^{3/2}$ $1/x$

$b$

$x$

Bearbeiten Sie, um die zusätzlichen Informationen zu berücksichtigen:

$y^* = b + a x^*$

Wir haben jetzt, dass die Fehler in x und additiv sind. Wir wissen immer noch nicht, ob die Varianz auf dieser Skala konstant ist.

$x^* = y^*/a - b/a = m y^* + c$

Let $x_o^* = x^* + \eta$ , where this error term may be heteroskedastic (if the original $x$ has constant spread, it will be heteroskedastic, but of known form)

(where the $o$ in $x_o^*$ stands for 'observed')

Then $x^*_o = c + m y^* + \epsilon$ where $\epsilon = -\zeta$ looks nice but now has correlated errors in the $x$ and $y$ variables; so it's a linear errors-in-variables model, with heteroskedasticity and known form of dependence in the errors.

I am not sure that improves things! I believe there are methods for that kind of thing, but it's not really my area at all.

I mentioned in the comments that you might like to look at inverse regression, but the particular form of your function may preclude getting far with that.

You might even be stuck with trying fairly robust-to-errors-in-x methods in that linear form.

Now a huge question: if the errors are in x, how the heck were you fitting the nonlinear model? Were you just blindly minimizing the sum of squared errors in $y$ ? That might well be your problem.

I suppose one could try to rewrite the original thing as a model with errors in the $x$ and try to optimize the fit but I am not sure I see how to set that up right.

Glen_b -Reinstate Monica
quelle

Thanks! It's an interesting transformation, hadn't thought of that --- even though the errors are in

x

$x$ , I'll play around with it anyway!

onnodb

"even though the errors are in x" -- yikes, that's kind of important. You may want to check up on inverse regression.

Glen_b -Reinstate Monica

...or you could directly fit the model

x = \frac{1}{3} (\frac{2 y}{a} + \frac{2^{1 / 3} y^{2}}{{(27 a^{4} b^{2} - 2 a^{3} y^{3} + 3 \sqrt{3} \sqrt{27 a^{8} b^{4} - 4 a^{7} b^{2} y^{3}})}^{1 / 3}} + \frac{{(27 a^{4} b^{2} - 2 a^{3} y^{3} + 3 \sqrt{3} \sqrt{27 a^{8} b^{4} - 4 a^{7} b^{2} y^{3}})}^{1 / 3}}{2^{1 / 3} a^{2}})

$x = \frac{1}{3} \left(\frac{2 y}{a}+\frac{2^{1/3} y^2}{\left(27 a^4 b^2-2 a^3 y^3+3 \sqrt{3} \sqrt{27 a^8 b^4-4 a^7 b^2 y^3}\right)^{1/3}}+\frac{\left(27 a^4 b^2-2 a^3 y^3+3 \sqrt{3} \sqrt{27 a^8 b^4-4 a^7 b^2 y^3}\right)^{1/3}}{2^{1/3} a^2}\right)$ :-).

whuber

@whuber Hmm. Solving for the cubic, clever. If we write the original in terms of

x_{o}

$x_o$ where

x_{o}

$x_o$ is

x + ζ

$x + \zeta$ , this would leave us with

x = (t h a t m o n s t e r) + ϵ

$x = (\rm{that\,\, monster}) + \epsilon\,$ , (again with

ϵ = - ζ

$\epsilon = -\zeta$ ) which at least notionally can be done with nonlinear least squares. So that looks like it takes care of the error propagation properly. It might actually work if the OP was to use the linear form I was playing with (using some robust-to-errors-in-the-IV-and-hetero estimation) to get good starting values for the parameters and then try to use this nonlinear LS form to polish it.

Glen_b -Reinstate Monica

I believe linearizing the function

x (y)

$x(y)$ and (ironically) applying nonlinear (weighted) least squares would work, especially if the data were restricted to relatively small values of

y

$y$ where the curve is primarily determined by

b

$b$ .

whuber

After some more weeks of experimenting, a different technique seems to work the best in this particular case: Total Least Squares fitting. It's a variant of the usual (nonlinear) Least Squares fitting, but instead of measuring fit errors along just one of the axes (which causes problems in highly nonlinear cases such as this one), it takes both axes into account.

There's a plethora of articles, tutorials and books avaiable on the subject, although the nonlinear case is more elusive. There's even some MATLAB code available.

onnodb
quelle

Thanks for sharing this. I accept that it it might produce good-looking results in your case, but I have two concerns. The first you mention: how exactly does one apply total least squares/errors-in-variables regression/orthogonal regression/Deming regression to nonlinear fits? The second is that this approach does not seem appropriate for your data, in which

y

$y$ is measured essentially without error. When that's the case, you should not be allowing for residuals in the

y

$y$ variable and doing so ought to produce unreliable, biased results.

whuber

@whuber Thanks for expressing your concerns! Right now, I'm still working on running simulations to probe the reliability of TLS fitting for this problem. What I've seen thus far, though, is that TLS' consideration of both variables helps greatly in overcoming the high non-linearity of the model. Fits of simulated data are reliable and converge very well. More work needs to be done though, and I'll definitely have to stack your method up to this one, once we have more actual data available --- and look in detail into your concerns.

onnodb

OK--don't forget I have comparable concerns about the method I proposed!

whuber