Como encontrar um bom ajuste para o modelo semi-sinusoidal em R?

Quero assumir que a temperatura da superfície do mar no Mar Báltico é a mesma ano após ano e depois descrevê-la com um modelo de função / linear. A idéia que tive foi apenas inserir o ano como um número decimal (ou num_months / 12) e descobrir qual deveria ser a temperatura naquele momento. Jogando-o na função lm () em R, ele não reconhece dados sinusoidais, portanto apenas produz uma linha reta. Então, coloquei a função sin () dentro de um colchete I () e tentei alguns valores para ajustá-la manualmente, e isso se aproxima do que eu quero. Mas o mar está esquentando mais rápido no verão e depois esfriando mais devagar no outono ... Então, o modelo está errado no primeiro ano, fica mais correto depois de alguns anos e, no futuro, acho que fica mais e mais errado de novo.

Como posso obter R para estimar o modelo para mim, para que eu não precise adivinhar números sozinho? A chave aqui é que eu quero que produza os mesmos valores ano após ano, e não apenas esteja correto por um ano. Se eu soubesse mais sobre matemática, talvez pudesse calculá-la como algo como um Poisson ou Gaussiano em vez de sin (), mas também não sei como fazer isso. Qualquer ajuda para se aproximar de uma boa resposta seria muito apreciada.

Aqui estão os dados que eu uso e o código para mostrar os resultados até agora:

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

r regression time-series lm GaRyu
fonte

Respostas:

Isso pode ser feito com regressão linear -

Você só precisa de um e um $\sin$ $\cos$ termo em cada frequência.

A razão pela qual você pode usar um termo e em uma regressão linear para lidar com a sazonalidade com qualquer amplitude e fase é devido ao seguinte $\sin$ $\cos$ identidade trigonométrica :

Uma onda senoidal 'geral' com amplitude e fase , , pode ser escrita como a combinação linear onde e $A$ $\varphi$ $A \sin (x + \varphi)$ $a\sin x+b\cos x$ $a$ são tais que $b$ e $A=\sqrt{a^2+b^2}$ $\sin\varphi = \frac{b}{\sqrt{a^2+b^2}}$ . Vamos ver que os dois são equivalentes:

\begin{array}{rcl} uma pecado (x) + b \cos (x) & = & \sqrt{{uma}^{2} + b^{2}} (\frac{uma}{\sqrt{{uma}^{2} + b^{2}}} pecado (x) + \frac{b}{\sqrt{{uma}^{2} + b^{2}}} \cos (x)) \\ = & UMA [pecado (x) \cos (φ) + \cos (x) pecado (φ)] \\ = & UMA pecado (x + φ) . \end{array}

$\begin{eqnarray} a \sin(x) + b \cos(x) &=& \sqrt{a^2+b^2} \left(\frac{a}{\sqrt{a^2+b^2}} \sin(x) + \frac{b}{\sqrt{a^2+b^2}} \cos(x)\right)\\ &=& A\left[\sin(x)\cos(\varphi) + \cos(x)\sin(\varphi)\right]\\ &=& A\sin(x+\varphi)\,\text{.} \end{eqnarray}$

Aqui está o modelo 'básico':

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

[recorte]

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)

pecado ajuste

Edit: Nota importante - os $2\pi\,t$ $t$ $\omega$ $(2\pi/\omega)\, t$

Aqui está o modelo com o segundo harmônico:

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

[recorte]

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)

pecado ajuste 2

... e assim por diante, com 6*pi*ToYetc. Se houvesse um pouquinho de ruído nos dados, eu provavelmente pararia com este segundo modelo.

Com termos suficientes, você pode ajustar exatamente sequências periódicas assimétricas e até irregulares, mas os ajustes resultantes podem 'mexer'. Aqui está uma função assimétrica (é uma serra - dente de serra ) adicionada a uma versão em escala de sua função periódica), com terceiros (vermelho) e quarto (verde) harmônicos. O ajuste verde é, em média, um pouco mais próximo, mas "ondulado" (mesmo quando o ajuste passa por todos os pontos, o ajuste pode ser muito ondulado entre os pontos).

pecado ajuste 3 e 4

$\cos$ $\sin$ termo será zero), enquanto o cos alterna entre 1 e -1).

Se você deseja ajustes mais suaves do que essa abordagem produz em séries não suaves, convém examinar splines periódicos ajustes .

Ainda outra abordagem é usar manequins sazonais, mas a abordagem sin / cos geralmente é melhor se for uma função periódica suave.

Esse tipo de abordagem da sazonalidade também pode se adaptar a situações em que a sazonalidade está mudando, como o uso de sazonalidade trigonométrica ou simulada com modelos de espaço de estados.

Embora a abordagem do modelo linear discutida aqui seja simples de usar, uma vantagem da abordagem de regressão não linear do @COOLSerdash é que ela pode lidar com uma variedade muito maior de situações - você não precisa mudar muito antes de estar em uma situação linear. a regressão não é mais adequada, mas os mínimos quadrados não lineares ainda podem ser usados (ter um período desconhecido seria um desses casos).

Glen_b -Reinstate Monica
fonte

Impressionante! Obrigado, eu realmente deveria tentar aprender mais sobre métodos para lidar com frequências. Não entendo por que a parte cos é necessária, mas conhecer o princípio facilita a implementação.

GaRyu 31/05

@COOLSerdash - na verdade, eu gostaria que você não tivesse excluído sua resposta (de fato, eu a votei); tem a vantagem de trabalhar em uma variedade muito maior de circunstâncias; ajustar algumas coisas sobre o problema e você pode perder a linearidade - e então minha abordagem é inútil, mas a sua ainda funciona. Eu acho que há muito a ser dito por ser capaz de fazer dessa maneira.

Glen_b -Reinstate Monica

@Glen_b Ah, desculpe, achei que sua postagem tornava a minha redundante porque eu não usava a maneira padrão de lidar com o problema. Eu a desfiz.

COOLSerdash

\cos

$\cos$

Aquele não fui eu ... Você diz que o deslocamento de fase é como o nome do que está acontecendo, e isso acontece matematicamente. Mas, para você, é mais provável que o ponto principal seja que 31 de dezembro / 1º de janeiro seja uma origem arbitrária para a época do ano, devido a atrasos na resposta da temperatura a variações no recebimento de radiação. Portanto, offset de fase também é um nome para algo climatológico, o tempo da temperatura mínima e máxima em relação ao seu sistema de gravação. (É um detalhe menor, mas eu prefiro quantificar época do ano para 12 meses, 1/24, 3/24, ..., 23/24.)

Nick Cox

A temperatura que você fornece na sua pergunta se repete exatamente a cada ano. Suspeito que não sejam realmente temperaturas medidas ao longo de quatro anos. No seu exemplo, você não precisaria de um modelo, porque as temperaturas se repetem exatamente. Mas, caso contrário, você poderia usar a nlsfunção para ajustar uma curva senoidal:

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

Ajuste NLS

Mas o ajuste não é muito bom, especialmente no começo. Parece que seus dados não podem ser modelados adequadamente por uma simples curva senoidal. Talvez uma função trigonométrica mais complexa funcione?

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 fit 2

A curva vermelha ajusta melhor os dados. Com a nlsfunção, você pode inserir o modelo que achar apropriado.

Ou talvez você possa fazer uso do forecastpacote. No exemplo abaixo, presumi que a série cronológica começou em janeiro de 2010:

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

Como os dados são determinísticos, nenhuma faixa de confiança é mostrada.

COOLSerdash
fonte

Não há razão para mínimos quadrados não lineares aqui, não que não funcione razoavelmente bem. Calcule o pecado (2 * pi * ToY), cos (2 * pi * ToY) com antecedência e alimente-os lm()como qualquer outro preditor. Em outras palavras, lm()não precisa ver nenhuma trigonometria. No entanto, você pode precisar de outro modelo para capturar bem a assimetria marcada. Não sou usuário regular de R, mas frequentemente usei essa abordagem em outros lugares (consulte stata-journal.com/sjpdf.html?articlenum=st0116 ).

Nick Cox

@ NickCox Obrigado Nick, esse é um conselho muito útil. Vou atualizar minha resposta daqui a pouco.

COOLSerdash

O Glen foi mais rápido :)

COOLSerdash

@COOLserdash Eu nem vi o comentário de Nick Cox lá; veio enquanto eu estava gerando minha resposta. (Essa abordagem é bastante óbvia se você já viu alguma série de Fourier.)

Glen_b -Reinstate Monica

Como o @Glen_b implica, essa é uma abordagem padrão, mas não universalmente conhecida.

Nick Cox