1
$\begingroup$

I acquired data (motor adaptation =y in function of delays =t ) which I expect to look like a sine wave. I am trying (1) to fit a sine curve in my data and (2)to estimate the best model/parameters. I read several posts here here here but I am sill struggling.

Then I tried to use nls, in order to tune my model by looping through several parameters. (I aslo used nls2 + optimization which gave me similar results to this post )

CODE

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)

OUTPUT

I get an 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

QUESTION

-So this error seems to be because my initial parameters are wrong. Is this correct? But how can I estimate the best parameters if the majority of parameters do not work?

-Also, I do not understand why the RSS is always the same despite different parameters

-And why do I observe only 2 different AIC and 2 different BIC while the models are different?

Any kind of help would much appreciated, thanks

$\endgroup$
1
  • $\begingroup$ Some sine wave will fit your data perfectly. Read about aliasing for some intuition. So, unless you place an upper limit on the frequency, this is a useless exercise. $\endgroup$
    – whuber
    Commented Oct 14, 2020 at 16:14

1 Answer 1

1
$\begingroup$

You are fitting one and the same model using different starting parameters. The global optimum of a fit (assuming it exists) is independent of starting parameters. If you get different parameters, either these parameters combined result in the exact same fit or some of these fits did not converge to the global optimum.

You could see this much better, if you fixed the line that stores the residual sum of squares to RSS[ii] <- sum(resid(res)^2), which leads to this output:

#         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

This post does show a nice approach for deriving a good starting value for 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)

resulting plot depicting the fitted sine wave

Fitting trigonometric functions is always challenging and it's best to use approaches from signal processing. I'm not convinced that the model you are trying to fit is a good model for this data. It would be good if you had more frequent measurements along the t axis.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.