Cómo encontrar un buen ajuste para la semi-sinusoidal modelo en R?

Question

Quiero suponer que la temperatura de la superficie marina del Mar Báltico es la misma año tras año, y, a continuación, describir con una función / modelo lineal. La idea que yo tenía era sólo la entrada de año como un número decimal (o num_months/12) y obtener lo que la temperatura debe ser de alrededor de ese tiempo. Lo arrojaba lm() en función de R, que no reconoce sinusoidal de datos, por lo que sólo se produce una línea recta. Así que me puse el pecado() función dentro de un I() soporte y trató de un par de valores para ajustar de forma manual la función, y que se acerca a lo que quiero. Pero el mar se está calentando más rápido en el verano y luego enfriamiento más lento en el otoño... Así que el modelo está mal el primer año, y luego se pone más correcto después de un par de años, y, a continuación, en el futuro, supongo que se vuelve más y más mal de nuevo.

¿Cómo puedo obtener R para estimar el modelo para mí, así que no tienes que adivinar los números de mí mismo? La clave aquí es que lo quiero para producir los mismos valores, año tras año, no sólo de ser la correcta para un año. Si supiera más acerca de las matemáticas, tal vez yo podría conjeturar que es algo como una Poisson o de Gauss en lugar de pecado(), pero no sé cómo hacer eso. Cualquier ayuda para acercarnos a una buena respuesta sería muy apreciada.

Estos son los datos que yo tengo, y el código para mostrar los resultados hasta el momento:

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

Answer 1

Se puede hacer con la regresión lineal -

Usted sólo necesita una $\sin$ $\cos$ plazo en cada frecuencia.

---

La razón por la que usted puede utilizar un $\sin$ $\cos$ plazo en una regresión lineal para controlar la estacionalidad con alguna amplitud y la fase es debido a la siguiente identidad trigonométrica:

Un 'general' de la onda sinusoidal con amplitud $A$ y la fase $\varphi$, $A \sin (x + \varphi)$, puede ser escrito como la combinación lineal $a\sin x+b\cos x$ donde $a$ $b$ son tales que $A=\sqrt{a^2+b^2}$$\sin\varphi = \frac{b}{\sqrt{a^2+b^2}}$. Vamos a ver que las dos son equivalentes:

\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}

---

Aquí está el 'básicas' del modelo:

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

[snip]

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)

---

Edit: nota Importante - la $2\pi\,t$ término las obras debido a que el periodo de la función ha sido configurada de tal forma que un período = 1 unidad de $t$. Si el período es diferente de 1, dicen que el período es $\omega$, entonces usted necesita $(2\pi/\omega)\, t$ lugar.

---

Aquí está el modelo con el segundo armónico:

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

[snip]

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)

... y así sucesivamente, con 6*pi*ToY etc. Si había un poco de ruido en los datos probablemente me deje con este segundo modelo, aunque.

Con suficiente términos que usted puede ajustarse exactamente asimétrica e incluso irregulares periódico de las secuencias, pero el resultado se ajusta puede 'meneo'. Aquí está una función asimétrica (es un diente de sierra ) que se agrega a una versión a escala de su función periódica), con terceros (rojo) y la cuarta (verde) armónicos. El verde de ajuste es, en promedio, un poco más cerca, pero "ondulante" (incluso cuando el ajuste pasa a través de cada punto, el ajuste puede ser muy floja entre los puntos).

La periodicidad con que se significa aquí sólo hay 12 d.f. disponible para una temporada en el modelo de datos. Con el intercepto en el modelo, usted sólo tiene suficientes grados de libertad para 11 adicional de temporada parámetros. Puesto que usted está agregando dos términos con cada uno de los armónicos, la última armónico que puede caber sólo permite uno de ellos para el último término, el sexto armónico (y que uno tiene que ser $\cos$; el $\sin$ plazo será todo de cero, mientras que el cos alterna entre 1 y -1).

Si desea cabe que son más suaves que este enfoque produce en la no-suave de la serie, puede que desee ver en el periódico de la spline se ajusta.

Sin embargo, otro enfoque es el uso de dummies estacionales, pero la sin/cos enfoque a menudo es mejor si es un buen periódico de la función.

Este tipo de enfoque a la estacionalidad también puede adaptarse a situaciones donde la estacionalidad es cambiante, como el uso trigonométricas o ficticio de estacionalidad con modelos de espacio de estado.

---

Mientras que el modelo lineal enfoque discutido aquí es simple de usar, una de las ventajas de @COOLSerdash de regresión no lineal de enfoque es que se puede tratar con una gama mucho más amplia de situaciones - usted no tiene que cambiar mucho antes de que usted esté en una situación en la regresión lineal ya no es el más adecuado, pero no lineal de mínimos cuadrados puede ser utilizado (teniendo un período desconocido sería uno de esos casos).

Answer 2

La temperatura que usted proporcione en su pregunta se repite exactamente cada año. Sospecho que esto no es realmente mide temperaturas de más de cuatro años. En tu ejemplo, no se necesita un modelo, debido a que la temperatura de repetir exactamente. Pero de lo contrario se podría utilizar la nls de la función de ajuste de una curva sinusoidal:

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

Pero el ajuste no es muy bueno, sobre todo al principio. Parece que sus datos no pueden ser adecuadamente modelados por un simple curva sinusoidal. Tal vez una más compleja de la función trigonométrica de hacer el truco?

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

La curva roja se ajusta a los datos mejor. Con la nls función, usted puede poner en el modelo que usted piensa que es apropiado.

O tal vez usted podría hacer uso de la forecast paquete. En el ejemplo de abajo, he asumido que el tiempo de la serie comenzó en enero del 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")

Porque los datos es determinista, no hay confianza bandas se muestran.