Strategia dopasowania wysoce nieliniowej funkcji

Do analizy danych z eksperymentu biofizyki próbuję obecnie dopasować krzywą za pomocą wysoce nieliniowego modelu. Funkcja modelu wygląda następująco:

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

Tutaj szczególnie duże znaczenie ma wartość . $b$

Wykres dla tej funkcji:

Wykres funkcji

(Zauważ, że funkcja modelu opiera się na dokładnym matematycznym opisie systemu i wydaje się działać bardzo dobrze --- po prostu zautomatyzowane dopasowania są trudne).

Oczywiście funkcja modelu jest problematyczna: strategie dopasowania, które wypróbowałem do tej pory, zawodzą z powodu ostrej asymptoty przy , szczególnie przy zaszumionych danych. $x=0$

Rozumiem ten problem tutaj: proste dopasowanie najmniejszych kwadratów (grałem z regresją zarówno liniową, jak i nieliniową w MATLAB; głównie Levenberg-Marquardt) jest bardzo wrażliwe na asymptotę pionową, ponieważ małe błędy w x są bardzo wzmocnione .

Czy ktoś mógłby wskazać mi odpowiednią strategię, która mogłaby obejść ten problem?

Mam podstawową wiedzę na temat statystyki, ale wciąż jest dość ograniczona. Chciałbym się uczyć, gdybym tylko wiedział, od czego zacząć :)

Bardzo dziękuję za radę!

Edytuj Błagam o wybaczenie za zapomnienie o wspominaniu błędów. Jedyny znaczący hałas występuje w i jest on addytywny. $x$

Edytuj 2 Dodatkowe informacje o tle tego pytania. Powyższy wykres modeluje zachowanie rozciągania polimeru. Jak zauważył @whuber w komentarzach, potrzebujesz aby uzyskać wykres podobny do powyższego. $b \approx -200 a$

Co do tego, jak ludzie dopasowywali tę krzywą do tego momentu: wydaje się, że ludzie zazwyczaj odcinają pionową asymptotę, dopóki nie znajdą dobrego dopasowania. Wybór odcięcia jest jednak nadal arbitralny, co powoduje, że procedura dopasowania jest niewiarygodna i odtwarzalna.

Edycja 3 i 4 Naprawiono wykres.

curve-fitting nonlinear

— onnodb
źródło

Czy błędy występują w w czy w obu? W jakiej formie spodziewany jest hałas (multiplikatywny, addytywny itp.)?

x

$x$

y

$y$

— probabilityislogic

@onnodb: Obawiam się, czy to nie może w zasadniczy sposób zakwestionować solidności samego modelu? Bez względu na to, jaką strategię dopasowania zastosujesz, nie pozostanie bardzo wrażliwy? Czy możesz kiedykolwiek mieć zaufanie do takiego oszacowania dla ?

b

$b$

b

$b$

— curious_cat 13.03.2013

Niestety to nadal nie zadziała. Tam po prostu nie jest możliwe połączenie i , które nawet jakościowo odtworzenia wykres został sporządzony. (Oczywiście jest ujemne. musi być mniejsze niż najmniejsze nachylenie na wykresie, ale dodatnie, co stawia go w wąskim przedziale. Ale kiedy jest w tym przedziale, po prostu nie jest wystarczająco duży, aby pokonać ogromny ujemny skok na pochodzenie wprowadzone terminem ). Co narysowałeś? Dane? Jakaś inna funkcja?

a

$a$

b

$b$

b

$b$

a

$a$

a

$a$

b x^{1 / 2}

$b x^{1/2}$

— whuber

Dzięki, ale nadal jest źle. Rozszerzając styczną do tego wykresu do tyłu z dowolnego punktu gdzie , przechwycisz oś y w . Ponieważ skok w dół przy pokazuje, że jest ujemny, to przecięcie również musi być ujemne. Ale na twojej figurze jest całkowicie jasne, że większość takich przechwyceń jest dodatnich, rozciągając się aż do . Jest zatem matematycznie niemożliwe, aby równanie takie jak mogło opisać twoją krzywą , nawet w przybliżeniu. Co najmniej musisz dopasować coś takiego jak .

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

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

x > 0

$x\gt 0$

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

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

0

$0$

b

$b$

15.5

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

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

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

— whuber

Zanim zacząłem nad tym pracować, chciałem się upewnić, że pytanie jest takie: ważne jest, aby funkcja była poprawna. Nie mam teraz czasu na udzielenie pełnej odpowiedzi, ale chciałbym zauważyć, że „inni ludzie” mogą się mylić - ale niestety zależy to od jeszcze więcej szczegółów. Jeśli twój błąd jest naprawdę addytywny , wydaje mi się, że nadal musi być silnie heteroscedastyczny, ponieważ w przeciwnym razie jego wariancja przy małych wartościach byłaby naprawdę niewielka. Co możesz nam powiedzieć - ilościowo - o tym błędzie?

x

$x$

x

$x$

— whuber

Odpowiedzi:

Metody, których użylibyśmy, aby dopasować to ręcznie (to znaczy eksploracyjną analizę danych), mogą zadziałać wyjątkowo dobrze z takimi danymi.

Chciałbym nieco sparametryzować model , aby jego parametry były dodatnie:

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

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

Dla danego , załóżmy, że istnieje unikalny rzeczywisty spełniający to równanie; nazwij to lub, dla zwięzłości, gdy zrozumiemy . $y$ $x$ $f(y; a,b)$ $f(y)$ $(a,b)$

Obserwujemy zbiór uporządkowanych par gdzie odbiegają od przez niezależne losowe zmienne ze średnimi zerowymi. W tej dyskusji założę, że wszystkie mają wspólną wariancję, ale rozszerzenie tych wyników (przy użyciu ważonych najmniejszych kwadratów) jest możliwe, oczywiste i łatwe do wdrożenia. Oto symulowany przykład takiej kolekcji wartości wartości , i wspólnej wariancji . $(x_i, y_i)$ $x_i$ $f(y_i; a,b)$ $100$ $a=0.0001$ $b=0.1$ $\sigma^2=4$

Wykres danych

Jest to (celowo) trudny przykład, co można docenić przez niefizyczne (ujemne) wartości i ich niezwykły rozkład (który zwykle wynosi jednostki poziome , ale może wynosić do lub na osi ). Jeśli uda nam się uzyskać rozsądne dopasowanie do tych danych, które są bliskie oszacowania zastosowanych wartości , i , naprawdę zrobimy to dobrze. $x$ $\pm 2$ $5$ $6$ $x$ $a$ $b$ $\sigma^2$

Wyposażenie eksploracyjne jest iteracyjne. Każdy etap składa się z dwóch etapów: oszacowanie (w oparciu o dane i wcześniejszych oszacowań i z i , z których poprzednie wartości przewidywanych można uzyskać za ), a następnie oszacuj . Ponieważ błędy są w x , dopasowania szacują na podstawie , a nie na odwrót. Aby najpierw uporządkować błędy w , gdy jest wystarczająco duże, $a$ $\hat{a}$ $\hat{b}$ $a$ $b$ $\hat{x}_i$ $x_i$ $b$ $x_i$ $(y_i)$ $x$ $x$

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

W związku z tym, możemy aktualizować montując ten model z najmniejszych kwadratów (zawiadomienie to ma tylko jeden parametr - zboczu, --and nie przechwycenia) oraz biorąc odwrotność współczynnika jako zaktualizowanego oszacowania . $\hat{a}$ $a$ $a$

Następnie, gdy jest wystarczająco małe, dominuje odwrotny kwadratowy wyraz i stwierdzamy (ponownie w celu pierwszego uporządkowania błędów), że $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}.$

Po ponownym użyciem metody najmniejszych kwadratów (z tylko zbocze termin ) otrzymujemy zaktualizowanej estymaty przez pierwiastek kwadratowy z zamontowanym zbocza. $b$ $\hat{b}$

Aby zrozumieć, dlaczego to działa, surowy rozpoznawcza przybliżenie tego ataku można uzyskać przez wykreślenie przed dla mniejszych . Jeszcze lepiej, ponieważ mierzy z błędem, a zmienia się monotonicznie z , powinniśmy skupić się na danych z większych wartości . Oto przykład z naszego symulowanego zbiorze pokazując największą połowę $x_i$ $1/y_i^2$ $x_i$ $x_i$ $y_i$ $x_i$ $1/y_i^2$ $y_i$ na czerwono, najmniejsza połowa na niebiesko, a linia przechodząca przez początek pasuje do czerwonych punktów.

Postać

Punkty w przybliżeniu się wyrównują, chociaż przy małych wartościach i występuje trochę krzywizny . (Zwróć uwagę na wybór osi: ponieważ jest pomiarem, zwykle wykreśla się go na osi pionowej .) Skupiając dopasowanie na czerwonych punktach, w których krzywizna powinna być minimalna, powinniśmy uzyskać rozsądne oszacowanie . Podana w tytule wartość jest pierwiastkiem kwadratowym nachylenia tej linii: to tylko % mniej niż prawdziwa wartość! $x$ $y$ $x$ $b$ $0.096$ $4$

W tym momencie prognozowane wartości można zaktualizować za pomocą

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

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

Iteruj, aż oszacowania ustabilizują się (co nie jest gwarantowane) lub przejdą przez małe przedziały wartości (których nadal nie można zagwarantować).

Okazuje się, że jest trudne do oszacowania, chyba że mamy dobry zestaw bardzo dużych wartości , ale to co determinuje pionową asymptotę na oryginalnym wykresie (w pytaniu) i jest w centrum pytania - można przypiąć dość dokładnie, pod warunkiem, że w pionowej asymptocie znajdują się pewne dane. W naszym przykładzie systemem iteracje są zbieżne dla (który jest prawie dwukrotnie prawidłowa wartość ) oraz (który jest zbliżony do właściwej wartości $a$ $x$ $b$ $\hat{a} = 0.000196$ $0.0001$ $\hat{b} = 0.1073$ $0.1$ ). Ten wykres pokazuje dane jeszcze raz, na które nakładają się (a) prawdziwa krzywa w kolorze szarym (przerywanym) i (b) oszacowana krzywa w kolorze czerwonym (ciągłym):

Pasuje

To dopasowanie jest tak dobre, że trudno jest odróżnić prawdziwą krzywą od dopasowanej krzywej: pokrywają się prawie wszędzie. Nawiasem mówiąc, szacowana wariancja błędu jest bardzo zbliżona do prawdziwej wartości . $3.73$ $4$

Z tym podejściem wiążą się pewne problemy:

Szacunki są tendencyjne. Odchylenie staje się widoczne, gdy zestaw danych jest mały i stosunkowo niewiele wartości znajduje się w pobliżu osi x. Dopasowanie jest systematycznie trochę niskie.
Procedura szacowania wymaga metody powiedzieć „duże” z „małych” wartości . Mógłbym zaproponować eksploracyjne sposoby identyfikacji optymalnych definicji, ale w praktyce można je pozostawić jako stałe „dostrajające” i zmieniać je, aby sprawdzić czułość wyników. Ustawiłem je arbitralnie, dzieląc dane na trzy równe grupy zgodnie z wartością używając dwóch zewnętrznych grup. $y_i$ $y_i$
Procedura nie będzie działać dla wszystkich możliwych kombinacji i lub wszystkich możliwych zakresów danych. Powinien on jednak działać dobrze, ilekroć w zestawie danych jest wystarczająca liczba krzywych, aby odzwierciedlić obie asymptoty: pionową na jednym końcu i ukośną na drugim końcu. $a$ $b$

Kod

W Mathematica napisano poniżej .

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}
   ];

Zastosuj to do danych (podanych przez wektory równoległe xi yuformowanych w macierz dwukolumnową data = {x,y}) aż do zbieżności, zaczynając od oszacowań : $a=b=0$

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

— Whuber
źródło

To niesamowita odpowiedź; Jestem bardzo zobowiązany! Bawiłem się tym, a wyniki wyglądają bardzo obiecująco. Potrzebuję jednak trochę więcej czasu, aby w pełni zrozumieć uzasadnienie :) Również: czy mogę skontaktować się z tobą za pośrednictwem Twojej witryny w celu uzyskania jednego dodatkowego (prywatnego) pytania dotyczącego potwierdzeń?

— onnodb

Zobacz ważne pytania opublikowane przez @probabilityislogic

$y^* = y\sqrt{x}$ $y^*$ $x^* = x^{3/2}$ $1/x$

$b$

$x$

Edytuj, aby uwzględnić dodatkowe informacje:

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

Mamy teraz, że błędy są w x i addytywnie. Nadal nie wiemy, czy wariancja jest stała w tej skali.

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

$x_o^* = x^* + \eta$ $x$

$o$ $x_o^*$

$x^*_o = c + m y^* + \epsilon$ $\epsilon = -\zeta$ $x$ $y$

Nie jestem pewien, czy to poprawi rzeczy! Sądzę, że istnieją metody tego typu, ale tak naprawdę to nie jest mój obszar.

Wspomniałem w komentarzach, że możesz chcieć spojrzeć na regresję odwrotną, ale szczególna forma twojej funkcji może uniemożliwić daleko do tego.

Możesz nawet utknąć w próbach dość solidnych metod na błędy w x w tej liniowej formie.

$y$

$x$

— Glen_b - Przywróć Monikę
źródło

x

$x$

„ nawet jeśli błędy są w x ” - tak, to trochę ważne. Możesz sprawdzić regresję odwrotną.

— Glen_b

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

x_{o}

$x_o$

x_{o}

$x_o$

x + ζ

$x + \zeta$

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

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

ϵ = - ζ

$\epsilon = -\zeta$

x (y)

$x(y)$

y

$y$

b

$b$

Po kilku tygodniach eksperymentów w tym konkretnym przypadku najlepsza wydaje się inna technika: dopasowanie Total Least Squares . Jest to wariant zwykłego (nieliniowego) dopasowania najmniejszych kwadratów, ale zamiast mierzyć błędy dopasowania wzdłuż jednej z osi (co powoduje problemy w wysoce nieliniowych przypadkach, takich jak ta), bierze pod uwagę obie osie.

Istnieje mnóstwo artykułów, samouczków i książek na ten temat, chociaż nieliniowy przypadek jest bardziej nieuchwytny. Jest nawet dostępny kod MATLAB .

— onnodb
źródło

y

$y$

y

$y$

@whuber Dziękujemy za wyrażenie swoich obaw! W tej chwili wciąż pracuję nad uruchomieniem symulacji, aby sprawdzić niezawodność dopasowania TLS do tego problemu. Do tej pory widziałem jednak, że uwzględnienie obu zmiennych przez TLS bardzo pomaga w przezwyciężeniu wysokiej nieliniowości modelu. Pasje symulowanych danych są niezawodne i bardzo dobrze zbiegają się. Trzeba jednak wykonać więcej pracy i na pewno będę musiał połączyć twoją metodę do tej, kiedy będziemy mieli więcej dostępnych danych - i przyjrzymy się szczegółowo twoim obawom.

— onnodb

OK - nie zapomnij, że mam podobne obawy dotyczące zaproponowanej przeze mnie metody!

— whuber