Ajuste de la curva sinusoidal y ajuste del modelo mediante nls en R

Question

He adquirido datos (adaptación del motor =y en función de los retrasos =t ) que espero que se parezcan a una onda sinusoidal. Estoy intentando (1) ajustar una curva sinusoidal a mis datos y (2) estimar el mejor modelo/parámetros. He leído varios posts aquí, pero sigo teniendo problemas.

Luego traté de usar nls, con el fin de afinar mi modelo mediante un bucle a través de varios parámetros. (También usé nls2 + optimización que me dio resultados similares a este Correo electrónico: )

CÓDIGO

t<-c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
y<-c(0.310 ,0.630 ,0.430 ,0.245, 0.650 ,0.085 ,0.370, 0.560 ,0.250, 0.520)
A <- (max(y)-min(y)/2)
C<-((max(y)+min(y))/2)
pp <- expand.grid(omega=(c(2.094395, 1.570796,  1.256637)), phi=(-1:1), A=A, C=C) # omega = 2*pi/3, pi/2 , 2*pi/5
#View(pp)

fit_AIC<- vector()
fit_BIC<- vector()

coef_A<- vector()
coef_ome<- vector()
coef_phi<- vector()
coef_C<- vector()
RSS<-vector()

for (ii in 1:nrow(pp))
{
  res<- nls(y ~ A*sin(omega*t+phi)+C, data=data.frame(t,y), start=list(A=pp$A[ii],omega=pp$omega[ii],phi=pp$phi[ii],C=pp$C[ii]), trace = TRUE)
  fit_AIC[ii]<-AIC(res)
  fit_BIC[ii]<-BIC(res)

  coef_A[ii]<- coef(res)[1]
  coef_ome[ii]<- coef(res)[2]
  coef_phi[ii]<- coef(res)[3]
  coef_C[ii]<- coef(res)[4]
  RSS<-sum(resid(res)^2)

}

results<-data.frame(RSS, fit_AIC,  fit_BIC,  coef_A,  coef_ome,  coef_phi,  coef_C)
View(results)

SALIDA

Me da un error...

1.405742 :   0.607500  2.094395 -1.000000  0.367500
0.1448148 :   0.1563179  2.1441802 -0.9937729  0.4172079

...

0.05018035 :   0.2195573  2.2852482 -1.1577097  0.4114573
2.085664 :  0.607500 1.570796 1.000000 0.367500
0.3104012 :  0.01321257 1.60518024 0.83201816 0.40437498
0.3098916 :   0.0180852  3.0888764 -5.9933691  0.4060743
Error in nls(y ~ A * sin(omega * t + phi) + C, data = data.frame(t, y),  : 
  le pas 0.000488281 est devenu inférieur à 'minFactor' de 0.000976562

        RSS    fit_AIC    fit_BIC     coef_A  coef_ome    coef_phi    coef_C
1 0.05018035 -14.568398 -13.055473 0.21955754 2.2852455  -1.1576955 0.4114573
2 0.05018035   2.753153   4.266079 0.07487110 0.8575642   0.2299909 0.3916769
3 0.05018035   2.753153   4.266079 0.07487109 0.8575736   0.2299951 0.3916763
4 0.05018035 -14.568398 -13.055473 0.21955763 2.2852443  -1.1576894 0.4114573
5 0.05018035 -14.568398 -13.055473 0.21955729 2.2852490 -32.5736406 0.4114573
6 0.05018035   2.753153   4.266079 0.07487105 0.8575619   0.2300021 0.3916770
7 0.05018035 -14.568398 -13.055473 0.21955735 2.2852482  -1.1577097 0.4114573

PREGUNTA

-Así que este error parece deberse a que mis parámetros iniciales son incorrectos. ¿Es esto correcto? Pero, ¿cómo puedo estimar los mejores parámetros si la mayoría de los parámetros no funcionan?

-Además, no entiendo por qué el RSS es siempre el mismo a pesar de los diferentes parámetros

-¿Y por qué observo sólo 2 AIC y 2 BIC diferentes mientras los modelos son diferentes?

Cualquier tipo de ayuda sería muy apreciada, gracias

Answer 1

Usted está ajustando un mismo modelo utilizando diferentes parámetros de partida. El óptimo global de un ajuste (suponiendo que exista) es independiente de los parámetros de partida. Si obtienes diferentes parámetros, o bien estos parámetros combinados dan como resultado exactamente el mismo ajuste o algunos de estos ajustes no convergen al óptimo global.

Esto se podría ver mucho mejor, si se fija la línea que almacena la suma residual de los cuadrados en RSS[ii] <- sum(resid(res)^2) que conduce a esta salida:

#         RSS    fit_AIC    fit_BIC     coef_A  coef_ome    coef_phi    coef_C
#1 0.05018035 -14.568398 -13.055473 0.21955754 2.2852455  -1.1576955 0.4114573
#2 0.28366065   2.753153   4.266079 0.07487110 0.8575642   0.2299909 0.3916769
#3 0.28366065   2.753153   4.266079 0.07487109 0.8575736   0.2299951 0.3916763
#4 0.05018035 -14.568398 -13.055473 0.21955763 2.2852443  -1.1576894 0.4114573
#5 0.05018035 -14.568398 -13.055473 0.21955729 2.2852490 -32.5736406 0.4114573
#6 0.28366065   2.753153   4.266079 0.07487105 0.8575619   0.2300021 0.3916770
#7 0.05018035 -14.568398 -13.055473 0.21955735 2.2852482  -1.1577097 0.4114573

Este post muestra un buen enfoque para derivar un buen valor inicial para omega: https://stats.stackexchange.com/a/60997/11849

plot(y ~ t)

raw.fft = fft(y)
truncated.fft = raw.fft[seq(1, length(y)/2 - 1)]
truncated.fft[1] = 0
omega = which.max(abs(truncated.fft)) * 2 * pi / length(y)

res<- nls(y ~ A * sin(omega * t + phi) + C, data = data.frame(t, y), 
          start=list(A = A, omega = omega, phi=1, C=0.4), trace = TRUE)
summary(res)
curve(predict(res, newdata = data.frame(t = x)), add = TRUE)

El ajuste de las funciones trigonométricas siempre es un reto y lo mejor es utilizar enfoques del procesamiento de señales. No estoy convencido de que el modelo que intentas ajustar sea un buen modelo para estos datos. Sería bueno que tuvieras mediciones más frecuentes a lo largo del eje t.