Ajustar un término sinusoidal a los datos

26

Aunque leí esta publicación, todavía 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)=UNAsyonorte(ωt+ϕ)+do.

con las cuatro incógnitas , \ omega , \ phi y C a la misma.UNAωϕdo

El resto del aspecto 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.

Ajuste sinusoidal

Agradecería mucho cualquier ayuda.

Aclamaciones.

Pascal
fuente
¿Está tratando de ajustar una onda sinusoidal a los datos o está tratando de ajustar algún tipo de modelo armónico con un componente seno y coseno? Hay una función armónica en el paquete TSA en R que es posible que desee consultar. Ajuste su modelo usando eso y vea qué tipo de resultados obtiene.
Eric Peterson
55
¿Has probado diferentes valores iniciales? Su función de pérdida no es convexa, por lo que diferentes valores iniciales pueden conducir a diferentes soluciones.
Stefan Wager
1
Cuéntanos más sobre los datos. Por lo general, hay una periodicidad conocida, por lo que no es necesario estimarla a partir de los datos. ¿Es esta una serie temporal o algo más? Es mucho más fácil si puede ajustar términos seno y coseno separados por un modelo lineal.
Nick Cox
2
Tener un período desconocido hace que su modelo no sea lineal (se alude a dicho evento en la respuesta seleccionada en la publicación vinculada). Dado que, los otros parámetros son condicionalmente lineales; para algunas rutinas LS no lineales, esa información es importante y puede mejorar el comportamiento. Una opción podría ser utilizar métodos espectrales para obtener el período y la condición de eso; otro sería actualizar el período y los otros parámetros mediante una optimización no lineal y lineal, respectivamente, de manera iterativa.
Glen_b -Reinstale a Monica el
(Acabo de editar la respuesta allí para hacer que el caso particular de un período desconocido sea un ejemplo explícito de lo que puede hacerlo no lineal.)
Glen_b -Reinstate Monica el

Respuestas:

18

Si solo desea una buena estimación de y no le importa mucho 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

trama sinusoidal

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

---

Si desea tener una idea de la incertidumbre en , puede utilizar la probabilidad de perfil ( pdf1 , pdf2 : referencias sobre cómo obtener CI o SE aproximados de la probabilidad de perfil o sus variantes no son difíciles de localizar)ω

(Alternativamente, puede alimentar estas estimaciones en nls ... y comenzar ya convergido).

Glen_b -Reinstate a Monica
fuente
(+1) buena respuesta. Traté de ajustar el modelo lineal lm(y~sin(2*pi*t)+cos(2*pi*t)pero esto no funcionó (el costérmino siempre fue 1). Solo por curiosidad: ¿qué hacen las dos primeras líneas (sé que spectrumestima la densidad espectral)?
COOLSerdash
1
@COOLSerdash Sí, debes tener las unidades de como el período (como estaba en la pregunta vinculada) para que funcione. Debería volver y enfatizar eso en la otra respuesta. (ctd)t2*pi*t
Glen_b -Reinstate Monica
1
@COOLSerdash (ctd): la segunda línea encuentra la frecuencia asociada con el pico más grande del espectro y se invierte para identificar el período. Al menos en este caso (pero sospecho más ampliamente), los valores predeterminados esencialmente identifican el período que maximiza la probabilidad tan de cerca que eliminé los pasos que tenía para maximizar la probabilidad de perfil en la región alrededor de ese período. La función specen TSA puede ser mejor (parece tener más opciones, una de las cuales puede ser importante a veces), pero en este caso el pico principal estaba exactamente en el mismo lugar que con, spectrumasí que no me molesté.
Glen_b -Reinstate a Monica el
@Glen_b este método funciona de maravilla para mi caso de uso. También necesito para ajustar una curva cos (x), pero no funciona así ... He cambiado la reslma reslm <- lm(y ~ cos(2*pi/per*t)+tan(2*pi/per*t)), pero que no se ven bien. alguna pista?
Amit Kohli
¿Por qué tienes un término bronceado allí?
Glen_b -Reinstala a Monica
15

2π/ /20

Cuando puse eso en nlsla startlista, obtuve una curva que era mucho más razonable, aunque todavía tiene algunos sesgos sistemáticos.

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

Ajuste sinusoidal

Elegir un valor inicial automáticamente

Si desea elegir la frecuencia dominante, puede usar una transformación rápida de Fourier (FFT). Esto está fuera de mi área de especialización, por lo que dejaré que otras personas completen los detalles si lo desean (especialmente sobre los pasos 2 y 3), pero el Rcódigo a continuación 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á que jugar un poco con la escala del eje x.

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

David J. Harris
fuente
1
Gracias por esa pista. Solo para aclarar un poco: los datos son parte de una micromatriz en la que se midió la periodicidad de los genes a lo largo del tiempo, es decir, los datos mostrados son los datos de expresión de un gen. El problema ahora es que quiero aplicar este método a unos 40k genes, todos con diferentes periodicidades y amplitudes. Por lo tanto, es crucial que se encuentre un buen ajuste independientemente de las condiciones iniciales.
Pascal
1
@Pascal Consulte mis actualizaciones anteriores para obtener una recomendación para elegir automáticamente el valor inicial de omega.
David J. Harris
2
ϕunasi
Me pregunto dónde entran en juego los valores de x aquí. Claro que hace una diferencia para omega, si los valores de y dados están separados por 1 o por 5 x pasos, ¿no?
knub
1
Consejo de programación no relacionado con la pregunta: precaución al nombrar objetos R como foo.bar. Esto se debe a cómo R especifica los métodos para las clases .
Firebug
10

Como alternativa a lo que ya se ha dicho, vale la pena señalar que un modelo AR (2) de la clase de modelos ARIMA se puede utilizar para generar pronósticos con un patrón de onda sinusoidal.

yt=do+ϕ1yt-1+ϕ2yt-2+unat
doϕ1ϕ2unat

ϕ12+4 4ϕ2<0.

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

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

Para ver si tal modelo podría ajustarse a los datos, utilicé la auto.arima()función del paquete de pronóstico para averiguar si sugeriría un modelo AR (2). Resulta que la auto.arima()función sugiere un modelo ARMA (2,2); no es 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 ciclos estocásticos). Es decir, todavía podemos verificar la condición antes mencionada para ver si se producirán pronósticos de onda sinusoidal.

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

ϕ12+4 4ϕ2<0 01.73472+4 4(-0.8324)<0 0-0.3202914<0 0

La siguiente gráfica muestra la serie original, y, el ajuste del modelo ARMA (2,2) y 14 pronósticos fuera de muestra. Como se puede ver, los pronósticos fuera de muestra siguen un patrón de onda sinusoidal.

ingrese la descripción de la imagen aquí

Ten en cuenta dos cosas. 1) Este es solo un análisis muy rápido (usando una herramienta automatizada) y un tratamiento adecuado implicaría seguir la metodología de Box-Jenkins. 2) Los pronósticos ARIMA son buenos para el pronóstico a corto plazo, por lo que puede encontrar que los pronósticos a largo plazo de los modelos en las respuestas de @David J. Harris y @Glen_b son más confiables.

Por último, espero que esta sea una buena adición a algunas respuestas ya muy informativas.

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

Graeme Walsh
fuente
1

Los métodos actuales para ajustar una curva sin a un conjunto de datos determinado requieren una primera aproximación de los parámetros, seguido de un proceso interactivo. Este es 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), consulte 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

Jaquelin
fuente
0

Si conoce el punto más bajo y más alto de sus datos de aspecto coseno, puede usar esta función simple para calcular todos los coeficientes 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 temperatura a lo largo del día con una función coseno, ingresando las horas y los valores de temperatura para la hora más baja y 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")

La salida está abajo: Coseno calculado desde los puntos más bajos y más altos

IVIM
fuente
3
Este enfoque parece ser particularmente sensible a cualquier desviación de aspecto aleatorio del comportamiento sinusoidal puro, lo que lo haría inaplicable a casi cualquier conjunto de datos como el ilustrado en la pregunta. Posiblemente, podría usarse para proporcionar valores iniciales para algunos de los otros enfoques iterativos sugeridos en este hilo.
whuber
están de acuerdo, es el más simple, sería bueno para aproximación simple bajo ciertas suposiciones
IVIM
0

Otra opción es utilizar la función genérica optim o nls. He intentado ambos ninguno de ellos es completamente robusto

Las siguientes funciones toman los datos en y y calculan los parámetros.

calc.period <- function(y,t)
{     
   fs <- 1/(t[2]-t[1])
   ssp <- spectrum(y,plot=FALSE )  
   fN <- ssp$freq[which.max(ssp$spec)]
   per <- 1/(fN*fs)
   return(per)
 }

fit.sine<- function(y, t)
{ 
  data <- data.frame(x = as.vector(t), y=as.vector(y))
  min.RSS <- function (data, par){
    with(data, sum((par[1]*sin(2*pi*par[2]*x + par[3])+par[4]-y )^2))
  }  
  amp = sd(data$y)*2.**0.5
  offset = mean(data$y)
  fest <- 1/calc.period(y,t)
  guess = c( amp, fest,  0,   offset)
  #res <- optim(par=guess, fn = min.RSS, data=data ) 
  r<-nls(y~offset+A*sin(2*pi*f*t+phi), 
     start=list(A=amp, f=fest, phi=0, offset=offset))
  res <- list(par=as.vector(r$m$getPars()))
  return(res)
}

 genSine <- function(t, params)
     return( params[1]*sin(2*pi*params[2]*t+ params[3])+params[4])

El uso es el siguiente:

t <- seq(0, 10, by = 0.01)
A <- 2 
f <- 1.5
phase <- 0.2432
offset <- -2

y <- A*sin(2*pi*f*t +phase)+offset + rnorm(length(t), mean=0, sd=0.2)

reslm1 <- fit.sine(y = y, t= t)

El siguiente código compara los datos.

ysin <- genSine(as.vector(t), params=reslm1$par)
ysin.cor <- genSine(as.vector(t), params=c(A, f, phase, offset))

plot(t, y)
lines(t, ysin, col=2)
lines(t, ysin.cor, col=3)
NMech
fuente