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

37

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))
GaRyu
źródło

Odpowiedzi:

44

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

sincos

sincos

AφAsin(x+φ)asinx+bcosxabA=a2+b2sinφ=ba2+b2

asin(x)+bcos(x)=a2+b2(aa2+b2sin(x)+ba2+b2cos(x))=A[sin(x)cos(φ)+cos(x)sin(φ)]=Asin(x+φ).

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)

grzech pasuje

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πttω(2π/ω)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)

grzech pasuje 2

... 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 - 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).

grzech pasuje 3 i 4

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).grzechcossin

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 ).

Glen_b
źródło
Niesamowite! Dziękuję, naprawdę powinienem spróbować dowiedzieć się więcej na temat metod radzenia sobie z częstotliwościami. Nie do końca rozumiem, dlaczego potrzebna jest część Cos, ale znajomość zasady ułatwia jej wdrożenie.
GaRyu
@COOLSerdash - właściwie żałuję, że nie usunąłeś swojej odpowiedzi (rzeczywiście ją głosowałem); ma tę zaletę, że działa w znacznie szerszym zakresie okoliczności; Popraw kilka rzeczy na temat problemu, a możesz stracić liniowość - a wtedy moje podejście jest bezużyteczne, ale twoje nadal działa. Myślę, że wiele można powiedzieć o tym, że można to zrobić w ten sposób.
Glen_b
@Glen_b Ach, przepraszam, myślałem, że twój post zwolnił mój, ponieważ nie zastosowałem standardowego sposobu rozwiązania problemu. Cofnąłem to.
COOLSerdash
cos
1
To nie byłem ja… Mówisz przesunięcie fazowe, jakby to nazywało to, co się dzieje, i robi to matematycznie. Ale dla ciebie kluczową kwestią jest to, że 31 grudnia / 1 stycznia jest arbitralnym początkiem pory roku, biorąc pod uwagę opóźnienia w reakcji temperatury na zmiany odbioru promieniowania. Więc przesunięcie fazowe jest tutaj nazwą czegoś klimatycznego, czasu minimalnej i maksymalnej temperatury w stosunku do twojego systemu zapisu. (To drobny szczegół, ale wolę określać porę roku na 12 miesięcy jako 1/24, 3/24, ..., 23/24.)
Nick Cox
10

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")

Dopasowanie NLS

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")

NLS pasuje 2

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")

ARIMA

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

COOLSerdash
źródło
4
Nie ma powodu, dla którego nieliniowe najmniejsze kwadraty są tutaj, nie dlatego, że nie będą działać dobrze. Oblicz z góry sin (2 * pi * ToY), cos (2 * pi * ToY) z wyprzedzeniem i karm je lm()tak, jak wszystkie inne predyktory. Innymi słowy, lm()nie trzeba wcale widzieć trygonometrii. Jednak może być potrzebny inny model, aby dobrze uchwycić zaznaczoną asymetrię. Nie jestem zwykłym użytkownikiem R, ale często korzystałem z tego podejścia gdzie indziej (patrz stata-journal.com/sjpdf.html?articlenum=st0116 ).
Nick Cox
@NickCox Dzięki Nick, to bardzo pomocna rada. Za chwilę zaktualizuję swoją odpowiedź.
COOLSerdash
Glen był szybszy :)
COOLSerdash
1
@COOLserdash Nie widziałem tam nawet komentarza Nicka Coxa; przyszło, kiedy generowałem odpowiedź. (Takie podejście jest dość oczywiste, jeśli widziałeś jakąś serię Fouriera.)
Glen_b
2
Jak sugeruje @Glen_b, jest to standardowe podejście, które po prostu nie jest powszechnie znane.
Nick Cox