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

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 = \log_{10}$$y = \log_{10}$

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?$r^2$

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.

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?

$f$$L$

• $f$
• $L$
• $A$$A>0$$A$
• $B$
• $C$$C$$-C$
• $D$
• $E$$L$$E$$L$$f=0$

$E$$E$$f=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.

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$

$y$$x$

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

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:

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)

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$

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

$B$$C$$B$$B$$C$

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]

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

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.