Dopasuj sinusoidalny termin do danych

26

Chociaż czytam ten post, nadal nie mam pojęcia, jak zastosować to do moich danych i mam nadzieję, że ktoś może mi pomóc.

Mam następujące dane:

y <- c(11.622967, 12.006081, 11.760928, 12.246830, 12.052126, 12.346154, 12.039262, 12.362163, 12.009269, 11.260743, 10.950483, 10.522091,  9.346292,  7.014578,  6.981853,  7.197708,  7.035624,  6.785289, 7.134426,  8.338514,  8.723832, 10.276473, 10.602792, 11.031908, 11.364901, 11.687638, 11.947783, 12.228909, 11.918379, 12.343574, 12.046851, 12.316508, 12.147746, 12.136446, 11.744371,  8.317413, 8.790837, 10.139807,  7.019035,  7.541484,  7.199672,  9.090377,  7.532161,  8.156842,  9.329572, 9.991522, 10.036448, 10.797905)
t <- 18:65

A teraz chcę po prostu dopasować falę sinusoidalną

y(t)=Asin(ωt+ϕ)+C.

z czterech niewiadomych , , i do niego.ω ϕ C.AωϕC

Reszta mojego kodu wygląda następująco

res <- nls(y ~ A*sin(omega*t+phi)+C, data=data.frame(t,y), start=list(A=1,omega=1,phi=1,C=1))
co <- coef(res)

fit <- function(x, a, b, c, d) {a*sin(b*x+c)+d}

# Plot result
plot(x=t, y=y)
curve(fit(x, a=co["A"], b=co["omega"], c=co["phi"], d=co["C"]), add=TRUE ,lwd=2, col="steelblue")

Ale wynik jest naprawdę słaby.

Dopasowanie sinusoidalne

Byłbym bardzo wdzięczny za każdą pomoc.

Twoje zdrowie.

Pascal
źródło
Próbujesz dopasować falę sinusoidalną do danych, czy próbujesz dopasować jakiś model harmoniczny do sinusu i komponentu cosinus? W pakiecie TSA w R znajduje się funkcja harmonicznych, którą możesz chcieć sprawdzić. Dopasuj swój model za pomocą tego i zobacz, jakie wyniki uzyskasz.
Eric Peterson
5
Czy próbowałeś różnych wartości początkowych? Twoja funkcja utraty nie jest wypukła, więc różne wartości początkowe mogą prowadzić do różnych rozwiązań.
Stefan Wager
1
Powiedz nam więcej o danych. Zazwyczaj znana jest okresowość, więc nie trzeba jej szacować na podstawie danych. Czy to szereg czasowy czy coś innego? Jest to o wiele łatwiejsze, jeśli można dopasować oddzielne warunki sinus i cosinus za pomocą modelu liniowego.
Nick Cox
2
Nieznany okres sprawia, że ​​Twój model jest nieliniowy (o takim zdarzeniu wspomina się w wybranej odpowiedzi w łączonym poście). Biorąc to pod uwagę, pozostałe parametry są warunkowo liniowe; w przypadku niektórych nieliniowych procedur LS ta informacja jest ważna i może poprawić zachowanie. Jedną z opcji może być użycie metod spektralnych w celu uzyskania okresu i warunku; innym byłoby zaktualizowanie okresu i innych parametrów odpowiednio poprzez optymalizację nieliniową i liniową w sposób iteracyjny.
Glen_b
(Właśnie edytowałem tam odpowiedź, aby konkretny przypadek nieznanego okresu był wyraźnym przykładem tego, co może uczynić go nieliniowym.)
Glen_b

Odpowiedzi:

18

Jeśli chcesz tylko dobrze oszacować i nie przejmujesz się jego standardowym błędem:ω

ssp <- spectrum(y)  
per <- 1/ssp$freq[ssp$spec==max(ssp$spec)]
reslm <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t))
summary(reslm)

rg <- diff(range(y))
plot(y~t,ylim=c(min(y)-0.1*rg,max(y)+0.1*rg))
lines(fitted(reslm)~t,col=4,lty=2)   # dashed blue line is sin fit

# including 2nd harmonic really improves the fit
reslm2 <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t)+sin(4*pi/per*t)+cos(4*pi/per*t))
summary(reslm2)
lines(fitted(reslm2)~t,col=3)    # solid green line is periodic with second harmonic

fabuła sinusoidy

(Nadal lepsze dopasowanie może w jakiś sposób tłumaczyć wartości odstające w tej serii, zmniejszając ich wpływ).

---

Jeśli chcesz mieć pojęcie o niepewności w , możesz użyć prawdopodobieństwa profilu ( pdf1 , pdf2 - odniesienia do uzyskania przybliżonych CI lub SE z prawdopodobieństwa profilu lub jego wariantów nie są trudne do zlokalizowania)ω

(Alternatywnie, możesz podać te szacunki do nls ... i rozpocząć już zbieżne).

Glen_b - Przywróć Monikę
źródło
(+1) ładna odpowiedź. Próbowałem dopasować model liniowy, lm(y~sin(2*pi*t)+cos(2*pi*t)ale to nie zadziałało ( costermin zawsze wynosił 1). Z ciekawości: co robią pierwsze dwie linie (wiem, że spectrumszacuje gęstość widmową)?
COOLSerdash,
1
t2*pi*t
1
@COOLSerdash (ctd) - 2. linia znajduje częstotliwość powiązaną z największym pikiem w widmie i odwraca się w celu identyfikacji okresu. Przynajmniej w tym przypadku (ale podejrzewam, że szerzej), wartości domyślne w nim zasadniczo określają okres, który maksymalizuje prawdopodobieństwo tak dokładnie, że usunąłem kroki, które miałem, aby zmaksymalizować prawdopodobieństwo profilu w regionie w tym okresie. Funkcja specw TSA może być lepsza (wydaje się, że ma więcej opcji, z których jedna może być czasem ważna), ale w tym przypadku główny szczyt znajdował się dokładnie w tym samym miejscu co, spectrumwięc nie zawracałem sobie głowy.
Glen_b
@Glen_b ta metoda działa cuda w moim przypadku użycia. Ja też trzeba dopasować cos (x) krzywą, ale to nie działa jak dobrze ... Zmieniłem reslmsię reslm <- lm(y ~ cos(2*pi/per*t)+tan(2*pi/per*t)), ale to nie wygląda dobrze. jakieś wskazówki?
Amit Kohli
Dlaczego masz tam opaleniznę?
Glen_b
15

2π/20

Kiedy kładę, że w nls„s startlisty, mam krzywą, która była znacznie bardziej rozsądne, chociaż nadal ma pewne systematyczne uprzedzeń.

W zależności od tego, jaki jest cel tego zestawu danych, możesz spróbować poprawić dopasowanie, dodając dodatkowe warunki lub stosując podejście nieparametryczne, takie jak proces Gaussa z okresowym jądrem.

Dopasowanie sinusoidalne

Automatyczny wybór wartości początkowej

Jeśli chcesz wybrać dominującą częstotliwość, możesz użyć szybkiej transformaty Fouriera (FFT). To jest poza moim obszarem specjalizacji, więc pozwolę innym osobom wypełnić szczegóły, jeśli chcą (szczególnie o krokach 2 i 3), ale Rponiższy kod powinien działać.

# Step 1: do the FFT
raw.fft = fft(y)

# Step 2: drop anything past the N/2 - 1th element.
# This has something to do with the Nyquist-shannon limit, I believe
# (https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem)
truncated.fft = raw.fft[seq(1, length(y)/2 - 1)]

# Step 3: drop the first element. It doesn't contain frequency information.
truncated.fft[1] = 0

# Step 4: the importance of each frequency corresponds to the absolute value of the FFT.
# The 2, pi, and length(y) ensure that omega is on the correct scale relative to t.
# Here, I set omega based on the largest value using which.max().
omega = which.max(abs(truncated.fft)) * 2 * pi / length(y)

Możesz również wykreślić, abs(truncated.fft)czy istnieją inne ważne częstotliwości, ale będziesz musiał nieco pogrzebać przy skalowaniu osi x.

Uważam też, że @Glen_b ma rację, że problem jest wypukły, gdy znasz omegę (a może też musisz znać phi? Nie jestem pewien). W każdym razie znajomość wartości początkowych dla innych parametrów nie powinna być tak ważna, jak w przypadku omegi, jeśli są one na właściwym miejscu. Prawdopodobnie można uzyskać przyzwoite oszacowania innych parametrów z FFT, ale nie jestem pewien, jak to by działało.

David J. Harris
źródło
1
Dzięki za podpowiedź. Wystarczy trochę wyjaśnić: dane są częścią mikromacierzy, w której okresowość genów była mierzona w czasie, tj. Pokazane dane są danymi dotyczącymi ekspresji jednego genu. Problem polega na tym, że chcę zastosować tę metodę do około 40 000 genów, z których wszystkie mają różne okresowości i amplitudy. Dlatego bardzo ważne jest, aby znaleźć dobre dopasowanie niezależnie od warunków początkowych.
Pascal
1
@Pascal Zobacz moje aktualizacje powyżej, aby uzyskać zalecenia dotyczące automatycznego wyboru wartości początkowej dla omega.
David J. Harris
2
ϕab
Zastanawiam się, gdzie wchodzą tutaj wartości x. Jasne, że ma to znaczenie dla omegi, niezależnie od tego, czy podane wartości y są oddzielone 1, czy 5 x krokami, prawda?
knub
1
Wskazówka programistyczna niezwiązana z pytaniem: ostrożność przy nazywaniu obiektów R jako foo.bar. Wynika to ze sposobu, w jaki R określa metody dla klas .
Firebug
10

Jako alternatywę do tego, co już powiedziano, warto zauważyć, że model AR (2) z klasy modeli ARIMA może być wykorzystywany do generowania prognoz z wzorem fali sinusoidalnej.

yt=C+ϕ1yt1+ϕ2yt2+at
Cϕ1ϕ2at

ϕ12+4ϕ2<0.

Panratz (1991) mówi nam o cyklach stochastycznych:

Stochastyczny wzór cyklu można pomyśleć o zniekształconym wzorze fali sinusoidalnej we wzorze prognozy: Jest to fala sinusoidalna z okresem stochastycznym (probabilistycznym), amplitudą i kątem fazowym.

Aby sprawdzić, czy taki model można dopasować do danych, skorzystałem z auto.arima()funkcji z pakietu prognozy, aby dowiedzieć się, czy sugerowałby model AR (2). Okazuje się, że auto.arima()funkcja sugeruje model ARMA (2,2); nie jest to czysty model AR (2), ale to jest OK. Jest OK, ponieważ model ARMA (2,2) zawiera element AR (2), więc obowiązuje ta sama zasada (o cyklach stochastycznych). Oznacza to, że nadal możemy sprawdzić wyżej wspomniany warunek, aby sprawdzić, czy zostaną wygenerowane prognozy fali sinusoidalnej.

Wyniki auto.arima(y)pokazano poniżej.

Series: y 
ARIMA(2,0,2) with non-zero mean 

Coefficients:
         ar1      ar2      ma1     ma2  intercept
      1.7347  -0.8324  -1.2474  0.6918    10.2727
s.e.  0.1078   0.0981   0.1167  0.1911     0.5324

sigma^2 estimated as 0.6756:  log likelihood=-60.14
AIC=132.27   AICc=134.32   BIC=143.5

ϕ12+4ϕ2<01.73472+4(0.8324)<00.3202914<0

Poniższy wykres przedstawia oryginalną serię y, dopasowanie modelu ARMA (2,2) i 14 prognoz poza próbą. Jak można zauważyć, prognozy poza próbą są zgodne z wzorem fali sinusoidalnej.

wprowadź opis zdjęcia tutaj

Pamiętaj o dwóch rzeczach. 1) To tylko bardzo szybka analiza (przy użyciu zautomatyzowanego narzędzia), a właściwe leczenie wymagałoby zastosowania metodologii Boxa-Jenkinsa. 2) Prognozy ARIMA są dobre w prognozowaniu krótkoterminowym, więc może okazać się, że prognozy długoterminowe z modeli w odpowiedziach @Davida J. Harrisa i @Glen_b są bardziej wiarygodne.

Wreszcie, mam nadzieję, że jest to miły dodatek do niektórych już bardzo pouczających odpowiedzi.

Odniesienie : Prognozowanie za pomocą modeli regresji dynamicznej: Alan Pankratz, 1991, (John Wiley and Sons, New York), ISBN 0-471-61528-5

Graeme Walsh
źródło
1

Obecne metody dopasowania krzywej grzechu do danego zestawu danych wymagają pierwszego odgadnięcia parametrów, a następnie procesu interakcyjnego. Jest to problem regresji nieliniowej. Inna metoda polega na przekształceniu regresji nieliniowej w regresję liniową dzięki wygodnemu równaniu całkowemu. Wtedy nie ma potrzeby wstępnego odgadywania i nie ma potrzeby iteracyjnego procesu: dopasowanie jest uzyskiwane bezpośrednio. W przypadku funkcji y = a + r * sin (w * x + phi) lub y = a + b * sin (w * x) + c * cos (w * x), patrz strony 35-36 artykułu „Régression sinusoidale” opublikowane na Scribd: http://www.scribd.com/JJacquelin/documents W przypadku funkcji y = a + p * x + r * sin (w * x + phi): strony 49–51 rozdziału „Mieszane regresje liniowe i sinusoidalne”. W przypadku bardziej skomplikowanych funkcji ogólny proces wyjaśniono w rozdziale „Uogólniona regresja sinusoidalna” na stronach 54–61, a następnie w przykładzie numerycznym y = r * sin (w * x + phi) + (b / x) + c * ln (x), strony 62-63

JJacquelin
źródło
0

Jeśli znasz najniższy i najwyższy punkt danych wyglądających na cosinus, możesz użyć tej prostej funkcji do obliczenia wszystkich współczynników cosinus:

getMyCosine <- function(lowest_point=c(pi,-1), highest_point=c(0,1)){
  cosine <- list(
    T = pi / abs(highest_point[1] - lowest_point[1]),
    b = - highest_point[1],
    k = (highest_point[2] + lowest_point[2]) / 2,
    A = (highest_point[2] - lowest_point[2]) / 2
  )
  return(cosine)
}

Poniżej służy do symulacji zmian temperatury w ciągu dnia za pomocą funkcji cosinus, wprowadzając godziny i wartości temperatury dla najniższej i najcieplejszej godziny:

c <- getMyCosine(c(4,10),c(17,25)) 
# lowest temprature at 4:00 (10 degrees), highest at 17:00 (25 degrees)

x = seq(0,23,by=1);  y = c$A*cos(c$T*(x +c$b))+c$k ; 
library(ggplot2);   qplot(x,y,geom="step")

Dane wyjściowe są poniżej: Cosinus obliczany od najniższych i najwyższych punktów

IVIM
źródło
3
Podejście to wydaje się szczególnie wrażliwe na wszelkie przypadkowo wyglądające odstępstwa od czysto sinusoidalnego zachowania, co uniemożliwiłoby zastosowanie go do prawie wszystkich zbiorów danych, takich jak ten przedstawiony w pytaniu. Można sobie wyobrazić, że można go użyć do podania wartości początkowych dla niektórych innych podejść iteracyjnych sugerowanych w tym wątku.
whuber
zgadzam się, jest najprostszy, przydałoby się pewne przybliżenie przy pewnych założeniach
IVIM
0

Inną opcją jest użycie funkcji ogólnej optim lub nls. Próbowałem obu, ale żaden z nich nie jest całkowicie solidny

Poniższe funkcje pobierają dane w y i obliczają parametry.

calc.period <- function(y,t)
{     
   fs <- 1/(t[2]-t[1])
   ssp <- spectrum(y,plot=FALSE )  
   fN <- ssp$freq[which.max(ssp$spec)]
   per <- 1/(fN*fs)
   return(per)
 }

fit.sine<- function(y, t)
{ 
  data <- data.frame(x = as.vector(t), y=as.vector(y))
  min.RSS <- function (data, par){
    with(data, sum((par[1]*sin(2*pi*par[2]*x + par[3])+par[4]-y )^2))
  }  
  amp = sd(data$y)*2.**0.5
  offset = mean(data$y)
  fest <- 1/calc.period(y,t)
  guess = c( amp, fest,  0,   offset)
  #res <- optim(par=guess, fn = min.RSS, data=data ) 
  r<-nls(y~offset+A*sin(2*pi*f*t+phi), 
     start=list(A=amp, f=fest, phi=0, offset=offset))
  res <- list(par=as.vector(r$m$getPars()))
  return(res)
}

 genSine <- function(t, params)
     return( params[1]*sin(2*pi*params[2]*t+ params[3])+params[4])

zastosowanie jest następujące:

t <- seq(0, 10, by = 0.01)
A <- 2 
f <- 1.5
phase <- 0.2432
offset <- -2

y <- A*sin(2*pi*f*t +phase)+offset + rnorm(length(t), mean=0, sd=0.2)

reslm1 <- fit.sine(y = y, t= t)

Poniższy kod porównuje dane

ysin <- genSine(as.vector(t), params=reslm1$par)
ysin.cor <- genSine(as.vector(t), params=c(A, f, phase, offset))

plot(t, y)
lines(t, ysin, col=2)
lines(t, ysin.cor, col=3)
NMech
źródło