Jak znaleźć dobre dopasowanie do modelu półsinusoidalnego w R?

Chcę założyć, że temperatura powierzchni morza w Morzu Bałtyckim jest taka sama rok po roku, a następnie opisać to za pomocą modelu funkcyjnego / liniowego. Pomysł, jaki miałem, to po prostu wpisać rok jako liczbę dziesiętną (lub num_months / 12) i ustalić, jaka powinna być temperatura w tym czasie. Wrzucając go do funkcji lm () w R, nie rozpoznaje danych sinusoidalnych, więc po prostu tworzy linię prostą. Umieściłem więc funkcję sin () w nawiasie I () i wypróbowałem kilka wartości, aby ręcznie dopasować funkcję, i to zbliża się do tego, czego chcę. Ale morze rozgrzewa się szybciej latem, a potem stygnie wolniej jesienią ... Więc model jest błędny w pierwszym roku, potem poprawia się po kilku latach, a potem wydaje mi się, że w przyszłości będzie bardziej i jeszcze raz źle.

Jak sprawić, by R oszacował dla mnie model, więc nie muszę zgadywać liczb? Kluczem tutaj jest to, że chcę, aby produkowała te same wartości rok po roku, a nie tylko przez jeden rok. Gdybym wiedział więcej o matematyce, może mógłbym zgadnąć, że to coś w rodzaju Poissona lub Gaussa zamiast sin (), ale nie wiem też, jak to zrobić. Będziemy wdzięczni za każdą pomoc w zbliżeniu się do dobrej odpowiedzi.

Oto dane, których używam, i kod do pokazania dotychczasowych wyników:

# SST from Bradtke et al 2010
ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)
SSTlm <- lm(SST$Degrees ~ I(sin(pi*2.07*SST$ToY)))
summary(SSTlm)
plot(SST,xlim=c(0,4),ylim=c(0,17))
par(new=T)
plot(data.frame(ToY=SST$ToY,Degrees=8.4418-6.9431*sin(2.07*pi*SST$ToY)),type="l",xlim=c(0,4),ylim=c(0,17))

Można to zrobić za pomocą regresji liniowej -

$\sin$$\cos$

$\sin$$\cos$

$A$$\varphi$$A \sin (x + \varphi)$$a\sin x+b\cos x$$a$$b$$A=\sqrt{a^2+b^2}$$\sin\varphi = \frac{b}{\sqrt{a^2+b^2}}$

Oto „podstawowy” model:

 SSTlm <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY),data=SST)
 summary(SSTlm)

[fantastyczna okazja]

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)              8.292      0.135   61.41   <2e-16 *** 
sin(2 * pi * ToY)       -5.916      0.191  -30.98   <2e-16 ***  
cos(2 * pi * ToY)       -4.046      0.191  -21.19   <2e-16 *** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9355 on 45 degrees of freedom
Multiple R-squared: 0.969,      Adjusted R-squared: 0.9677 
F-statistic: 704.3 on 2 and 45 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylim=c(1.5,16.5),data=SST)
 lines(SST$ToY,SSTlm$fitted,col=2)

Edycja: Ważna uwaga - termin działa, ponieważ okres funkcji został ustawiony tak, że jeden okres = 1 jednostka . Jeśli kropka różni się od 1, powiedzmy, że kropka to , wtedy potrzebujesz zamiast.t ω ( 2 π / ω $2\pi\,t$$t$$\omega$$(2\pi/\omega)\, t$

Oto model z drugą harmoniczną:

 SSTlm2 <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY)
                        +sin(4*pi*ToY)+cos(4*pi*ToY),data=SST)
 summary(SSTlm2)

[fantastyczna okazja]

Coefficients:
                  Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        8.29167    0.02637  314.450  < 2e-16 ***  
sin(2 * pi * ToY) -5.91562    0.03729 -158.634  < 2e-16 ***  
cos(2 * pi * ToY) -4.04632    0.03729 -108.506  < 2e-16 ***  
sin(4 * pi * ToY)  1.21244    0.03729   32.513  < 2e-16 ***  
cos(4 * pi * ToY)  0.33333    0.03729    8.939 2.32e-11 ***  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1827 on 43 degrees of freedom
Multiple R-squared: 0.9989,     Adjusted R-squared: 0.9988 
F-statistic:  9519 on 4 and 43 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylab="Degrees",xlab="ToY",ylim=c(1.5,16.5),data=SST)
 lines(SSTlm2$fitted~ToY,col=2,data=SST)

... i tak dalej, z 6*pi*ToYitd. Gdyby w danych był niewielki szum, prawdopodobnie zatrzymałbym się przy tym drugim modelu.

Przy wystarczającej liczbie terminów można dokładnie dopasować asymetryczne, a nawet postrzępione sekwencje okresowe, ale powstałe dopasowania mogą „drgać”. Oto funkcja asymetryczna (to piłokształtny - ) dodana do skalowanej wersji funkcji okresowej), z trzecimi (czerwonymi) i czwartymi (zielonymi) harmonicznymi. Zielony krój jest średnio nieco bliżej, ale „kręci się” (nawet gdy dopasowanie przechodzi przez każdy punkt, dopasowanie może być bardzo kruche między punktami).

Okresowość oznacza tutaj, że w danych dostępnych jest tylko 12 df dla modelu sezonowego. Z przecięciem w modelu masz tylko wystarczającą liczbę stopni swobody dla 11 dodatkowych parametrów sezonowych. Ponieważ dodajesz dwa terminy do każdej harmonicznej, ostatnia harmoniczna, którą możesz dopasować, pozwoli ci tylko na jedną z nich na ostatnią kadencję, szóstą harmoniczną (i ta musi być ; termin będzie all- zero, podczas gdy cos zmienia się między 1 a -1).$\cos$$\sin$

Jeśli chcesz, aby dopasowania były gładsze niż w przypadku tego podejścia w przypadku serii nieładnych, możesz zajrzeć do okresowych dopasowań splajnów .

Jeszcze innym podejściem jest stosowanie sezonowych manekinów, ale podejście sin / cos jest często lepsze, jeśli jest to płynna funkcja okresowa.

Tego rodzaju podejście do sezonowości można również dostosować do sytuacji, w których sezonowość się zmienia, takich jak stosowanie trygonometrycznej lub pozornej sezonowości w modelach przestrzeni stanów.

Chociaż omówione tutaj podejście do modelu liniowego jest proste w użyciu, jedną z zalet nieliniowej metody regresji @ COOLSerdash jest to, że może on obsługiwać znacznie szerszy zakres sytuacji - nie musisz dużo zmieniać, zanim znajdziesz się w sytuacji, w której liniowy regresja nie jest już odpowiednia, ale nadal można stosować nieliniowe najmniejsze kwadraty ( jednym z takich przypadków byłby nieznany okres ).

Temperatura podana w pytaniu jest powtarzana dokładnie co roku. Podejrzewam, że tak naprawdę nie są to zmierzone temperatury przez cztery lata. W twoim przykładzie nie potrzebujesz modelu, ponieważ temperatury tylko powtarzają się dokładnie. Ale w przeciwnym razie możesz użyć nlsfunkcji, aby dopasować krzywą sinusoidalną:

ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)

par(cex=1.5, bg="white")
plot(Degrees~ToY,xlim=c(0,4),ylim=c(0,17), pch=16, las=1)

nls.mod <-nls(Degrees ~ a + b*sin(2*pi*c*ToY), start=list(a = 1, b = 1, c=1))

co <- coef(nls.mod) 
f <- function(x, a, b, c) {a + b*sin(2*pi*c*x) }

curve(f(x, a=co["a"], b=co["b"], c=co["c"]), add=TRUE ,lwd=2, col="steelblue")

Ale dopasowanie nie jest zbyt dobre, szczególnie na początku. Wygląda na to, że twoich danych nie można odpowiednio modelować za pomocą prostej krzywej sinusoidalnej. Może załatwi to bardziej złożona funkcja trygonometryczna?

nls.mod2 <-nls(Degrees ~ a + b*sin(2*pi*c*ToY)+d*cos(2*pi*e*ToY), start=list(a = 1, b = 1, c=1, d=1, e=1))

co2 <- coef(nls.mod2) 
f <- function(x, a, b, c, d, e) {a + b*sin(2*pi*c*x)+d*cos(2*pi*e*x) }

curve(f(x, a=co2["a"], b=co2["b"], c=co2["c"], d=co2["d"], e=co2["e"]), add=TRUE ,lwd=2, col="red")

Czerwona krzywa lepiej pasuje do danych. Dzięki tej nlsfunkcji możesz wprowadzić model, który Twoim zdaniem jest odpowiedni.

A może możesz skorzystać z forecastpakietu. W poniższym przykładzie założyłem, że szereg czasowy rozpoczął się w styczniu 2010 r .:

library(forecast)

Degrees.ts <- ts(Degrees, start=c(2010,1), frequency=12)

Degree.trend <- auto.arima(Degrees.ts)

degrees.forecast <- forecast(Degree.trend, h=12, level=c(80,95), fan=F)

plot(degrees.forecast, las=1, main="", xlab="Time", ylab="Degrees")

Ponieważ dane są deterministyczne, nie pokazano pasm ufności.