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

Otra opción es utilizar la función genérica optim o nls. He probado 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)