Dopasuj sinusoidalny termin do danych

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ą

z czterech niewiadomych , , i do niego.ω ϕ $A$$\omega$$\phi$$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.

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

Twoje zdrowie.

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

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

(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)$\omega$

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

$2 \pi / 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.

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.

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.

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

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.

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

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

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:

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)