Jak wybrać początkowe wartości dopasowania nieliniowego najmniejszego kwadratu

13

Powyższe pytanie mówi wszystko. Zasadniczo moje pytanie dotyczy ogólnej funkcji dopasowania (może być arbitralnie skomplikowane), która będzie nieliniowa w parametrach, które próbuję oszacować, w jaki sposób wybrać wartości początkowe, aby zainicjować dopasowanie? Próbuję robić nieliniowe najmniejsze kwadraty. Czy jest jakaś strategia lub metoda? Czy zostało to zbadane? Jakieś referencje? Czy jest coś poza zgadywaniem ad hoc? Konkretnie, obecnie jedną z form dopasowania, z którymi pracuję, jest forma Gaussa plus forma liniowa z pięcioma parametrami, które próbuję oszacować, takie jak

y=Ae(xBC)2+Dx+E

gdzie (dane odciętych) iy = log 10 (dane rzędnych), co oznacza, że ​​w przestrzeni log-log moje dane wyglądają jak linia prosta plus wybrzuszenie, które aproksymuję gaussowskim. Nie mam teorii, nic, co by mnie poprowadziło o tym, jak zainicjować dopasowanie nieliniowe, z wyjątkiem być może wykresów i gałek ocznych, takich jak nachylenie linii i jaki jest środek / szerokość wypukłości. Ale mam do wyboru ponad sto takich dopasowań zamiast grafowania i zgadywania, wolałbym jakieś podejście, które można zautomatyzować.x=log10y=log10

Nie mogę znaleźć żadnych odniesień w bibliotece ani w Internecie. Jedyne, co mogę wymyślić, to po prostu losowe wybranie wartości początkowych. MATLAB oferuje losowe wybranie wartości z [0,1] równomiernie rozłożonych. Więc z każdym zestawem danych uruchamiam losowo zainicjowane dopasowanie tysiąc razy, a następnie wybieram ten z najwyższym ? Jakieś inne (lepsze) pomysły?r2


Dodatek nr 1

Po pierwsze, oto kilka wizualnych reprezentacji zestawów danych, aby pokazać wam, o jakich danych mówię. Publikuję zarówno dane w oryginalnej postaci bez jakiejkolwiek transformacji, a następnie ich wizualną reprezentację w przestrzeni dziennika, ponieważ wyjaśnia niektóre cechy danych, a inne zniekształcają. Zamieszczam próbkę zarówno dobrych, jak i złych danych.

dobre dane log-loguj dobre dane złe dane zaloguj złe dane

Każdy z sześciu paneli na każdej figurze pokazuje cztery zestawy danych naniesione razem na czerwono, zielono, niebiesko i cyjanowo, a każdy zestaw danych ma dokładnie 20 punktów danych. Staram się dopasować do każdego z nich linię prostą plus gaussa z powodu nierówności widocznych w danych.

Pierwsza cyfra to niektóre z dobrych danych. Druga cyfra to wykres log-log tych samych dobrych danych z rysunku 1. Trzecia liczba to niektóre złe dane. Czwarta cyfra jest logarytmicznym wykresem z figury trzeciej. Jest o wiele więcej danych, to tylko dwa podzbiory. Większość danych (około 3/4) jest dobra, podobnie jak dobre dane, które pokazałem tutaj.

Teraz kilka komentarzy, proszę o wyrozumiałość, ponieważ może to potrwać długo, ale myślę, że wszystkie te szczegóły są konieczne. Spróbuję być możliwie zwięzły.

Pierwotnie spodziewałem się prawa prostej mocy (co oznacza linię prostą w przestrzeni log-log). Kiedy narysowałem wszystko w przestrzeni dziennika, zobaczyłem nieoczekiwany wzrost przy około 4,8 MHz. Guz został dokładnie zbadany i odkryto go również w innych pracach, więc nie zawiedliśmy go. Jest tam fizycznie i wspominają o tym także inne opublikowane prace. Więc właśnie dodałem termin gaussowski do mojej formy liniowej. Zauważ, że to dopasowanie miało być wykonane w przestrzeni log-log (stąd moje dwa pytania łącznie z tym).

Teraz, po przeczytaniu odpowiedzi przez stumpy Joe Pete'a do innej kwestii kopalni (niezwiązane z tych danych w ogóle) i czytając to i to i odniesienia w nim (stuff Clauset), zdaję sobie sprawę, że nie powinien zmieścić się w logarytmiczny przestrzeń. Więc teraz chcę robić wszystko w przestrzeni przekształconej.

Pytanie 1: Patrząc na dobre dane, nadal uważam, że liniowy plus gaussowski w przestrzeni wstępnie przetworzonej jest nadal dobrą formą. Chciałbym usłyszeć od innych, którzy mają większe doświadczenie w zakresie danych, co myślą. Czy gaussowski + liniowy jest rozsądny? Czy powinienem robić tylko gaussa? A może zupełnie inna forma?

Pytania 2: Niezależnie od odpowiedzi na pytanie 1 nadal potrzebowałbym (najprawdopodobniej) nieliniowego dopasowania najmniejszych kwadratów, więc nadal potrzebuję pomocy przy inicjalizacji.

Dane, w których widzimy dwa zestawy, bardzo wolimy uchwycić pierwszy guz przy częstotliwości około 4-5 MHz. Więc nie chcę dodawać więcej terminów gaussowskich, a nasz termin gaussowski powinien być wyśrodkowany na pierwszym guzie, który prawie zawsze jest większy. Chcemy „większej dokładności” między 0,8 MHz a około 5 MHz. Nie zależy nam zbytnio na wyższych częstotliwościach, ale nie chcemy też ich całkowicie ignorować. Więc może jakieś ważenie? Czy B może być zawsze inicjowany w okolicach 4,8 MHz?

fL

L=Ae(fBC)2+Df+E.
  • f
  • L
  • AA>0A
  • B
  • CCC
  • D
  • ELELf=0

Ae(B/C)2+E.

EEf=0

L

Pytania 3: Co według was ekstrapolujesz w ten sposób w tym przypadku? Jakieś zalety / wady? Jakieś inne pomysły na ekstrapolację? Ponownie dbamy tylko o niższe częstotliwości, dlatego ekstrapolujemy między 0 a 1 MHz ... czasami bardzo małe częstotliwości, zbliżone do zera. Wiem, że ten post jest już zapakowany. Zadałem to pytanie tutaj, ponieważ odpowiedzi mogą być powiązane, ale jeśli wolicie, mogę oddzielić to pytanie i zadać kolejne pytanie później.

Na koniec oto dwa przykładowe zestawy danych na żądanie.

0.813010000000000   0.091178000000000   0.012728000000000
1.626000000000000   0.103120000000000   0.019204000000000
2.439000000000000   0.114060000000000   0.063494000000000
3.252000000000000   0.123130000000000   0.071107000000000
4.065000000000000   0.128540000000000   0.073293000000000
4.878000000000000   0.137040000000000   0.074329000000000
5.691100000000000   0.124660000000000   0.071992000000000
6.504099999999999   0.104480000000000   0.071463000000000
7.317100000000000   0.088040000000000   0.070336000000000
8.130099999999999   0.080532000000000   0.036453000000000
8.943100000000001   0.070902000000000   0.024649000000000
9.756100000000000   0.061444000000000   0.024397000000000
10.569000000000001   0.056583000000000   0.025222000000000
11.382000000000000   0.052836000000000   0.024576000000000
12.194999999999999   0.048727000000000   0.026598000000000
13.008000000000001   0.045870000000000   0.029321000000000
13.821000000000000   0.041454000000000   0.067300000000000
14.633999999999999   0.039596000000000   0.081800000000000
15.447000000000001   0.038365000000000   0.076443000000000
16.260000000000002   0.036425000000000   0.075912000000000

Pierwsza kolumna to częstotliwości w MHz, identyczne w każdym zestawie danych. Druga kolumna to dobry zestaw danych (dobre dane pierwszy i drugi, panel 5, czerwony znacznik), a trzecia kolumna to zły zestaw danych (złe dane trzeci i czwarty, panel 5, czerwony marker).

Mam nadzieję, że to wystarczy, aby zachęcić do bardziej oświeconej dyskusji. Dziękuję wszystkim.

Fixed Point
źródło
+1 za dodatkowe informacje, ale teraz wygląda to na nowe pytanie. Nawiasem mówiąc, jeśli chcesz teraz usunąć poprzednią, myślę, że to byłoby w porządku, wygląda na to, że masz teraz dostęp do dodatkowych informacji.
Glen_b
@Glen_b Dlaczego tak jest? Dlaczego to wygląda na nowe pytanie? Co do starego pytania, wszyscy dziwimy się na punkty; - D, a stary ma dwa głosy poparcia, jakikolwiek sposób, aby połączyć to z tym, abym mógł zachować te dwa głosy również?
Fixed Point
Cóż, na początek, teraz pytasz o to, co powinieneś pasować, zamiast określać, co pasować, jak poprzednio. Istnieje wiele innych różnic, z których niektóre uważam za dość znaczące. Spojrzę na zmianę mojej odpowiedzi, ale myślę, że ta może stać się oryginalnym pytaniem i odpowiedzią, a twoje nowe części, w których pytasz o inne rzeczy, mogą być nowe. Na razie pozostawię to twojemu osądowi.
Glen_b
@Glen_b W porządku, wybrałem dodatkowe pytania. Więc pytania są nadal: mam pewne dane, które chcę dopasować przy użyciu formy liniowej + gaussowskiej, czy mogę lepiej niż przypadkowa inicjalizacja?
Fixed Point
Myślę, że moja obecna odpowiedź pokazuje, że - przynajmniej w niektórych okolicznościach - możesz zrobić lepiej, a @whuber sugeruje coś jeszcze prostszego niż mój proces. Mógłbym cofnąć się i zobaczyć, jak radzę sobie z danymi, ale nawet w obecnej postaci daje pewien pomysł, jak ustawić takie punkty początkowe.
Glen_b

Odpowiedzi:

10

Gdyby istniała strategia zarówno dobra, jak i ogólna - taka, która zawsze działała - byłaby już wdrożona w każdym nieliniowym programie najmniejszych kwadratów, a wartości początkowe nie byłyby problemem.

W przypadku wielu konkretnych problemów lub rodzin problemów istnieją pewne całkiem dobre podejścia do wartości początkowych; niektóre pakiety są wyposażone w dobre obliczenia wartości początkowej dla określonych modeli nieliniowych lub w bardziej ogólne podejścia, które często działają, ale może być konieczne uzyskanie bardziej szczegółowych funkcji lub bezpośrednie wprowadzenie wartości początkowych.

Eksplorowanie przestrzeni jest konieczne w niektórych sytuacjach, ale myślę, że twoja sytuacja może być taka, że ​​bardziej szczegółowe strategie będą prawdopodobnie warte zachodu - ale zaprojektowanie dobrej wymaga właściwie dużej wiedzy w dziedzinie, której prawdopodobnie nie będziemy w stanie posiadać.

x

yx

A

Pomocne byłyby niektóre przykładowe dane - typowe przypadki i trudne, jeśli możesz.


Edycja: Oto przykład, w jaki sposób możesz zrobić całkiem dobrze, jeśli problem nie jest zbyt głośny:

Oto niektóre dane generowane z twojego modelu (wartości populacji to A = 1,9947, B = 10, C = 2,828, D = 0,09, E = 5):

dane nls

Wartości początkowe, które udało mi się oszacować, to
(As = 1,658, Bs = 10,001, Cs = 3,053, Ds = 0,0881, Es = 5,026)

Dopasowanie tego modelu początkowego wygląda następująco:

nlstart

Kroki były następujące:

  1. Dopasuj regresję Theila, aby uzyskać przybliżone oszacowanie D i E.
  2. Odejmij dopasowanie regresji Theila
  3. Użyj LOESS, aby dopasować gładką krzywą
  4. Znajdź pik, aby uzyskać przybliżone oszacowanie A, a wartość x odpowiadająca pikowi, aby uzyskać przybliżone oszacowanie B
  5. Przyjmij pasowania LOESS, których wartości y są> 60% oszacowania A, jako obserwacje i pasuj do kwadratu
  6. Użyj kwadratu, aby zaktualizować oszacowanie B i oszacować C
  7. Od oryginalnych danych odejmij oszacowanie Gaussa
  8. Dopasuj kolejną regresję Theila do skorygowanych danych, aby zaktualizować oszacowanie D i E.

W takim przypadku wartości są bardzo odpowiednie do rozpoczęcia dopasowania nieliniowego.

Napisałem to jako Rkod, ale to samo można zrobić w MATLAB.

Myślę, że możliwe są lepsze rzeczy niż to.

Jeśli dane są bardzo hałaśliwe, nie zadziała to wcale dobrze.


Edycja2: To jest kod, którego użyłem w R, jeśli ktoś jest zainteresowany:

gausslin.start <- function(x,y) {

  theilreg <- function(x,y){
    yy <- outer(y, y, "-")
    xx <- outer(x, x, "-")
    z  <- yy / xx
    slope     <- median(z[lower.tri(z)])
    intercept <- median(y - slope * x)
    cbind(intercept=intercept,slope=slope)
  }

  tr <- theilreg(x,y1)
  abline(tr,col=4)
  Ds = tr[2]
  Es = tr[1]
  yf  <- y1-Ds*x-Es
  yfl <- loess(yf~x,span=.5)

  # assumes there are enough points that the maximum there is 'close enough' to 
  #  the true maximum

  yflf   <- yfl$fitted    
  locmax <- yflf==max(yflf)
  Bs     <- x[locmax]
  As     <- yflf[locmax]

  qs     <- yflf>.6*As
  ys     <- yfl$fitted[qs]
  xs     <- x[qs]-Bs
  lf     <- lm(ys~xs+I(xs^2))
  bets   <- lf$coefficients
  Bso    <- Bs
  Bs     <-  Bso-bets[2]/bets[3]/2
  Cs     <- sqrt(-1/bets[3])
  ystart <- As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es

  y1a <- y1-As*exp(-((x-Bs)/Cs)^2)
  tr  <- theilreg(x,y1a)
  Ds  <- tr[2]
  Es  <- tr[1]
  res <- data.frame(As=As, Bs=Bs, Cs=Cs, Ds=Ds, Es=Es)
  res
}

.

# population parameters: A = 1.9947 , B = 10, C = 2.828, D = 0.09, E = 5
# generate some data
set.seed(seed=3424921)
x  <- runif(50,1,30)
y  <- dnorm(x,10,2)*10+rnorm(50,0,.2)
y1 <- y+5+x*.09 # This is the data
xo <- order(x)

starts <- gausslin.start(x,y1)
ystart <- with(starts, As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es)
plot(x,y1)
lines(x[xo],ystart[xo],col=2)
Glen_b - Przywróć Monikę
źródło
3
+1. Powtarzanie dopasowania tysiąc razy i wybranie najlepszego (jeśli dobrze to rozumiem) brzmi dziwnie: nieliniowe najmniejsze kwadraty powinny się zbiegać, jeśli model jest odpowiedni dla danych i istnieją dobre wartości początkowe. Oczywiście drugie pytanie dotyczy tego, o co pytasz. Ale sugerowanie, że konieczne może być wybranie różnych wartości początkowych dla każdego dopasowania, wydaje się pesymistyczne.
Nick Cox
1
@NickCox Sprowadza się to do szeregu napotkanych problemów - jeśli przypominam sobie z wcześniejszych postów, OP otrzymuje ogromną liczbę tych problemów, ale nie przypominałem sobie, aby widziałem wystarczająco dużo szczegółów, aby przedstawić dobre sugestie, chociaż zainwestowałem trochę czas zabawy z potencjalnymi podejściami (które nie przyniosły niczego definitywnego do opublikowania). OP prawdopodobnie ma wiedzę w dziedzinie, która może przynieść dobre wartości początkowe, które prawie zawsze rozwiązują jego problemy.
Glen_b
1
Właśnie. Przegapiłem wcześniejszy post na stats.stackexchange.com/questions/61724/…
Nick Cox
3
|A|BA>0CA1/4A>0A<0
2
BB
6

Istnieje ogólne podejście do dopasowania tego rodzaju modeli nieliniowych. Polega ona na ponownym parametryzowaniu parametrów liniowych za pomocą wartości zmiennej zależnej, powiedzmy na pierwszej, ostatniej wartości częstotliwości, i na dobrym punkcie pośrodku, powiedzmy na 6-tym punkcie. następnie możesz utrzymać te parametry na stałym poziomie i rozwiązać dla parametru nieliniowego w pierwszej fazie minimalizacji, a następnie zminimalizować ogólne 5 parametrów.

Schnute i ja wymyśliliśmy to około 1982 r., Montując modele wzrostu dla ryb.

http://www.nrcresearchpress.com/doi/abs/10.1139/f80-172

Jednak nie jest konieczne czytanie tego artykułu. Z uwagi na fakt, że parametry są liniowe, po prostu konieczne jest ustawienie i rozwiązanie liniowego układu równań 3x3, aby zastosować stabilną parametryzację modelu.

M

M=(exp(((x(1)B)/C)2)x(1)1exp(((x(6)B)/C)2)x(6)1exp(((x(n)B)/C)2)x(n)1)
n=20
DATA_SECTION
  init_int n
  int mid
 !! mid=6;
  init_matrix data(1,n,1,3)
  vector x(1,n)
  vector y(1,n)
 !! x=column(data,1);
 !! y=column(data,3);   //use column 3
PARAMETER_SECTION
  init_number L1(3)     //(3) means estimate in phase 3
  init_number Lmid(3)
  init_number Ln(3)

  vector L(1,3)
  init_number log_B       // estimate in phase 1
  init_number log_C(2)    // estimate in phase 2 
  matrix M(1,3,1,3);
  objective_function_value f
  sdreport_vector P(1,3)
  sdreport_number B
  sdreport_number C
  vector pred(1,n);
PROCEDURE_SECTION
  L(1)=L1;
  L(2)=Lmid;
  L(3)=Ln;
  B=exp(log_B);
  C=exp(log_C);
  M(1,1)=exp(-square((x(1)-B)/C));
  M(1,2)=x(1);
  M(1,3)=1;
  M(2,1)=exp(-square((x(mid)-B)/C));
  M(2,2)=x(mid);
  M(2,3)=1;
  M(3,1)=exp(-square((x(n)-B)/C));
  M(3,2)=x(n);
  M(3,3)=1;

  P=solve(M,L);  // solve for standard parameters 
                 // P is vector corresponding to A,D,E

  pred=P(1)*exp(-square((x-B)/C))+P(2)*x+P(3);
  if (current_phase()<4)
    f+=norm2(y-pred);
  else
    f+=0.5*n*log(norm2(y-pred))  //concentrated likelihood

BCBBC

wprowadź opis zdjęcia tutaj

W twoim przypadku ze złymi danymi pasuje dość łatwo, a (zwykle) szacunkowe parametry to:

         estimate    std dev
A      2.0053e-01 5.8723e-02
D      1.6537e-02 4.7684e-03
E     -1.8197e-01 7.3355e-02
B      3.0609e+00 5.0197e-01
C      5.6154e+00 9.4564e-01]
Dave Fournier
źródło
Dave, to ciekawe, ale rodzi pytania. Dokładnie co rozumiesz przez „tego rodzaju modele nieliniowe”? Pytanie zaczyna się od odwołania się do „ogólnej funkcji dopasowania”, ale twój opis odnosi się tylko do „ogólnych 5 parametrów”.
whuber
Mam na myśli modele takie jak vonbertalanffy, na przykład logistyczne lub podwójnie wykładnicze. We wszystkich przypadkach model jest liniowy w niektórych parametrach i nieliniowy w innych. Ludzie zazwyczaj próbują je przekształcić, aby uzyskać bardziej stabilną parametryzację, koncentrując się na parametrach nieliniowych. Jest to jednak niewłaściwe podejście. Należy zmodyfikować parametryzację liniową. Na przykład dla logistyki 4-parametrowej model jest liniowy w górnej i dolnej asymptocie, ale zamiast używać tych parametrów, należy użyć przewidywanych wartości dla najmniejszego i największego ind. var.
dave fournier
@davefournier Dziękujemy za odpowiedź i wskazanie swojego dokumentu. Twój artykuł wydaje się trochę trudny do zdobycia, ale technika brzmi interesująco, więc nie mogę się doczekać, aby go przeczytać.
Naprawiono punkt
2

Jeśli musisz to robić wiele razy, sugerowałbym, abyś używał ewolucyjnego algorytmu w funkcji SSE jako interfejsu, aby podać wartości początkowe.

Z drugiej strony można użyć GEOGEBRA do utworzenia funkcji za pomocą suwaków parametrów i zabawy z nimi w celu uzyskania wartości początkowych.

LUB wartości wyjściowe z danych można oszacować na podstawie obserwacji.

  1. D i E pochodzą ze zbocza i przechwytują dane (ignorując Gaussa)
  2. A jest pionową odległością maksimum Gaussa od oszacowania linii Dx + E.
  3. B jest wartością x maksimum Gaussa
  4. C to połowa pozornej szerokości Gaussa
Adrian O'Connor
źródło
1

Dla wartości początkowych można zastosować zwykłe dopasowanie najmniejszych kwadratów. Jego nachylenie i punkt przecięcia byłyby wartościami początkowymi dla D i E. Największa wartość resztkowa byłaby wartością początkową dla A. Pozycja największej wartości resztkowej byłaby wartością początkową dla B. Może ktoś inny może zasugerować wartość początkową dla sigma.

Jednak nieliniowe najmniejsze kwadraty bez wyprowadzania jakiegokolwiek mechanistycznego równania ze znajomości przedmiotu jest ryzykownym biznesem, a robienie wielu osobnych dopasowań czyni sprawy jeszcze bardziej wątpliwymi. Czy za proponowanym równaniem kryje się jakaś wiedza przedmiotowa? Czy istnieją inne niezależne zmienne, które odnoszą się do różnic między 100 osobnymi dopasowaniami? Pomocne może być uwzględnienie tych różnic w jednym równaniu, które będzie pasować do wszystkich danych jednocześnie.

Emil Friedman
źródło