Ajustar un término sinusoidal a los datos

Question

Aunque he leído este post, aún no tengo idea de cómo aplicar esto a mis propios datos y espero que alguien pueda ayudarme.

Tengo los siguientes datos:

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

Y ahora simplemente quiero ajustar una onda sinusoidal

$$ y(t)=A\cdot sin(\omega t+\phi) +C. $$

con las cuatro incógnitas $A$ , $\omega$ , $\phi$ y $C$ a ella.

El resto de mi código es el siguiente

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

Pero el resultado es realmente pobre.

Agradecería mucho cualquier ayuda.

Salud.

Answer 1

Si sólo quiere una buena estimación de $\omega$ y no se preocupan mucho por su error estándar:

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

(Un ajuste mejor aún podría tener en cuenta los valores atípicos de esa serie de alguna manera, reduciendo su influencia).

---

Si quiere tener una idea de la incertidumbre en $\omega$ , podría utilizar la probabilidad del perfil ( pdf1 , pdf2 - no es difícil encontrar referencias sobre cómo obtener ICs o SEs aproximados a partir de la probabilidad del perfil o sus variantes)

(Como alternativa, se podrían introducir estas estimaciones en nls... y empezar ya convergiendo).

Answer 2

Como sugirió @Stefan, los diferentes valores de partida parecen mejorar el ajuste de forma espectacular. He calculado los datos para sugerir que omega debería ser de aproximadamente $2 \pi / 20$ ya que los picos parecían estar separados por unas 20 unidades.

Cuando pongo eso en nls 's start lista, obtuve una curva mucho más razonable, aunque sigue teniendo algunos sesgos sistemáticos.

Dependiendo de cuál sea su objetivo con este conjunto de datos, podría intentar mejorar el ajuste añadiendo términos adicionales o utilizando un enfoque no paramétrico como un proceso gaussiano con un núcleo periódico.

Elegir un valor inicial automáticamente

Si quieres elegir la frecuencia dominante, puedes utilizar una transformada rápida de Fourier (FFT). Esto está fuera de mi área de experiencia, así que dejaré que otras personas completen los detalles si lo desean (especialmente sobre los pasos 2 y 3), pero la R el código de abajo debería funcionar.

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

También puede trazar abs(truncated.fft) para ver si hay otras frecuencias importantes, pero tendrás que juguetear un poco con la escala del eje x.

Además, creo que @Glen_b está en lo cierto al afirmar que el problema es convexo una vez que se conoce omega (¿o tal vez también es necesario conocer phi? No estoy seguro). En cualquier caso, conocer los valores iniciales de los otros parámetros no debería ser tan importante como los de omega si están en el rango correcto. Probablemente podrías obtener estimaciones decentes de los otros parámetros a partir de la FFT, pero no estoy seguro de cómo funcionaría.

Answer 3

Como alternativa a lo que ya se ha dicho, cabe señalar que un modelo AR(2) de la clase de los modelos ARIMA puede utilizarse para generar previsiones con un patrón sinusoidal.

Un modelo AR(2) puede escribirse como sigue: \begin{equation} y_{t} = C + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + a_{t} \end{equation} donde $C$ es una constante, $\phi_{1}$ , $\phi_{2}$ son parámetros a estimar y $a_{t}$ es un término de choque aleatorio.

Ahora bien, no todos los modelos AR(2) producen patrones sinusoidales (también conocidos como ciclos estocásticos ) en sus previsiones, pero sí ocurre cuando se cumple la siguiente condición: \begin{equation} \phi_{1}^{2} + 4 \phi_{2} < 0. \end{equation}

Panratz(1991) nos dice lo siguiente sobre los ciclos estocásticos:

Un patrón de ciclo estocástico puede pensarse como una onda sinusoidal distorsionada distorsionada en el patrón de previsión: Es una onda sinusoidal con un período (probabilística) con un período, una amplitud y un ángulo de fase estocásticos.

Para ver si dicho modelo podía ajustarse a los datos, utilicé el auto.arima() del paquete de previsión para saber si sugiere un modelo AR(2). Resulta que la función auto.arima() sugiere un modelo ARMA(2,2); no un modelo AR(2) puro, pero esto está bien. Está bien porque un modelo ARMA(2,2) contiene un componente AR(2), por lo que se aplica la misma regla (sobre los ciclos estocásticos). Es decir, podemos seguir comprobando la condición mencionada para ver si se producen previsiones sinusoidales.

Los resultados de auto.arima(y) se muestran a continuación.

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

Ahora vamos a comprobar la condición: \begin{equation} \phi_{1}^{2} + 4 \phi_{2} < 0\\ 1.7347^{2} + 4 (-0.8324) < 0 \\ -0.3202914 < 0 \end{equation} y comprobamos que, efectivamente, la condición se cumple.

El gráfico siguiente muestra la serie original, y, el ajuste del modelo ARMA(2,2) y 14 previsiones fuera de muestra. Como puede verse, las previsiones fuera de muestra siguen un patrón de onda sinusoidal.

Tenga en cuenta dos cosas. 1) Este es sólo un análisis muy rápido (con una herramienta automatizada) y un tratamiento adecuado implicaría seguir la metodología Box-Jenkins. 2) Las previsiones ARIMA son buenas para las previsiones a corto plazo, por lo que es posible que las previsiones a largo plazo de los modelos de las respuestas de @David J. Harris y @Glen_b sean más fiables.

Por último, espero que esto sea un buen complemento a unas respuestas ya muy informativas.

Referencia : Previsión con modelos de regresión dinámica: Alan Pankratz, 1991, (John Wiley and Sons, Nueva York), ISBN 0-471-61528-5

Answer 4

Los métodos actuales para ajustar una curva de pecado a un conjunto de datos dado requieren una primera estimación de los parámetros, seguida de un proceso interativo. Se trata de un problema de regresión no lineal. Un método diferente consiste en transformar la regresión no lineal en una regresión lineal gracias a una ecuación integral conveniente. Entonces, no hay necesidad de conjeturas iniciales ni de procesos iterativos: el ajuste se obtiene directamente. En el caso de la función y = a + r*sin(w*x+phi) o y=a+b*sin(w*x)+c*cos(w*x), véanse las páginas 35-36 del documento "Régression sinusoidale" publicado en Scribd : http://www.scribd.com/JJacquelin/documents En el caso de la función y = a + p*x + r*sin(w*x+phi) : páginas 49-51 del capítulo "Regresiones lineales y sinusoidales mixtas". En el caso de funciones más complicadas, el proceso general se explica en el capítulo "Regresión sinusoidal generalizada" páginas 54-61, seguido de un ejemplo numérico y = r*sin(w*x+phi)+(b/x)+c*ln(x), páginas 62-63

Answer 5

Si conoce el punto más bajo y el más alto de sus datos con aspecto de coseno, puede utilizar esta sencilla función para calcular todos los coeficientes del coseno:

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

A continuación se utiliza para simular la variación de la temperatura a lo largo del día con una función coseno, introduciendo las horas y los valores de temperatura de la hora más baja y la más cálida:

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

El resultado es el siguiente: