# Fit a sinusoidal term to data

Although I read this post, I still have no idea how to apply this to my own data and hope that someone can help me out.

I have the following data:

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


And now I simply want to fit a sine wave

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

with the four unknowns $A$, $\omega$, $\phi$ and $C$ to it.

The rest of my code looks is the following

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


But the result is really poor.

I would very much appreciate any help.

Cheers.

-
You're trying to fit a sine wave to the data or are you trying to fit some kind of a harmonic model with a sine and a cosine component? There is a harmonic function in the TSA package in R that you might want to check out. Fit your model using that and see what kind of results you get. –  BackGammon Cube Jun 5 at 18:31
Have you tried different starting values? Your loss function is non-convex, so different starting values can lead to different solutions. –  Stefan Wager Jun 5 at 18:34
Tell us more about the data. Usually there is a known periodicity, so that need not be estimated from the data. Is this a time series or something else? It is much easier if you can fit separate sine and cosine terms by a linear model. –  Nick Cox Jun 5 at 18:46
Having an unknown period makes your model nonlinear (such an event is alluded to in the selected answer at the linked post). The given that, the other parameters are conditionally linear; for some nonlinear LS routines that information is important and can improve the behaviour. One option might be to use spectral methods to get the period and condition on that; another would be to update the period and the other parameters via a nonlinear and linear optimization respectively in an iterative fashion. –  Glen_b Jun 5 at 23:13
(I just edited the answer there to make the particular case of unknown period an explicit example of what can make it nonlinear.) –  Glen_b Jun 5 at 23:19

As @Stefan suggested, different starting values do seem to improve the fit dramatically. I eyeballed the data to suggest that omega should be about $2 \pi / 20$, since the peaks looked like they were about 20 units apart.

When I put that into nls's start list, I got a curve that was much more reasonable, although it still has some systematic biases.

Depending on what your goal is with this data set, you could try to improve the fit by adding additional terms or using a nonparametric approach like a Gaussian process with a periodic kernel.

## Choosing a starting value automatically

If you want to pick the dominant frequency, you can use a fast Fourier transform (FFT). This is way out of my area of expertise, so I'll let other folks fill in the details if they'd like (especially about steps 2 and 3), but the R code below should work.

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


You can also plot abs(truncated.fft) to see if there are other important frequencies, but you'll have to fiddle with the scaling of the x-axis a bit.

Also, I believe @Glen_b is correct that the problem is convex once you know omega (or maybe you need to know phi too? I'm not sure). In any case, knowing the starting values for the other parameters shouldn't be nearly as important as for omega if they're in the right ballpark. You could probably get decent estimates of the other parameters from the FFT, but I'm not certain how that would work.

-
Thanks for that hint. Just to clarify a little: the data is a part of a microarray in which the periodicity of genes was measured over time, i.e. the shown data is the expression data of one gene. The problem now is that I wanna apply this method to about 40k genes all having different periodicities and amplitudes. So, it's quite crucial that a good fit is found independent of the initial conditions. –  Pascal Jun 5 at 21:08
@Pascal See my updates above for a recommendation for automatically choosing the starting value for omega. –  David J. Harris Jun 6 at 0:25
@DavidJ.Harris You can estimate $\phi$ in a linear model as well (well, calculate it directly from $a$ and $b$ in the linear model), see the post the OP linked to. –  Glen_b Jun 6 at 3:51
@Glen_b cool, good to know. Thanks –  David J. Harris Jun 6 at 5:57

If you just want a good estimate of $\omega$ and don't care much about its standard error:

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  (A better fit still would perhaps account for the outliers in that series in some way, reducing their influence.) --- If you want some idea of the uncertainty in$\omega$, you could use profile likelihood (pdf1, pdf2 - references on getting approximate CIs or SEs from profile likelihood or its variants aren't hard to locate) (Alternatively, you could feed these estimates into nls ... and start it already converged.) - (+1) nice answer. I tried to fit the linear model with lm(y~sin(2*pi*t)+cos(2*pi*t) but this didn't work (cos term was always 1). Just out of curiosity: what do the first two lines do (I know that spectrum estimates the spectral density)? – COOLSerdash Jun 6 at 8:44 @COOLSerdash Yeah, you have to have the units of$t$being the period (as it was in the linked question) for 2*pi*t to work. I should go back and emphasize that in the other answer. (ctd) – Glen_b Jun 6 at 10:06 @COOLSerdash (ctd)- The 2nd line finds the frequency associated with the biggest peak in the spectrum and inverts to identify the period. At least in this case (but I suspect more widely), the defaults on it essentially identifies the period that maximizes the likelihood so closely that I deleted the steps I had in to maximize the profile likelihood in the region around that period. The function spec in TSA may be better (it seems to have more options, one of which may be important sometimes), but in this case the main peak was in exactly the same place as with spectrum so I didn't bother. – Glen_b Jun 6 at 10:07 Thanks a lot for your explanations. – COOLSerdash Jun 6 at 10:14 add comment As an alternative to what has already been said, it may be worth noting that an AR(2) model from the class of ARIMA models can be used to generate forecasts with a sine wave pattern. An AR(2) model can be written as follows: $$y_{t} = C + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + a_{t}$$ where$C$is a constant,$\phi_{1}$,$\phi_{2}$are parameters to be estimated and$a_{t}\$ is a random shock term.

Now, not all AR(2) models produce sine wave patterns (also known as stochastic cycles) in their forecasts, but it does happen when the following condition is satisfied: $$\phi_{1}^{2} + 4 \phi_{2} < 0.$$

Panratz(1991) tells us the following about stochastic cycles:

A stochastic cycle pattern can be thought of a distorted sine wave pattern in the forecast pattern: It is a sine wave with a stochastic (probabilistic) period, amplitude, and phase angle.

To see if such a model could be fitted to the data I used the auto.arima() function from the forecast package to find out if it would suggest an AR(2) model. It turns out that the auto.arima() function suggests an ARMA(2,2) model; not a pure AR(2) model, but this is OK. It's OK because an ARMA(2,2) model contains an AR(2) component, so the same rule (about stochastic cycles) applies. That is, we can still check the aforementioned condition to see if sine wave forecasts will be produced.

The results of auto.arima(y) are shown below.

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


Now let's check the condition: $$\phi_{1}^{2} + 4 \phi_{2} < 0\\ 1.7347^{2} + 4 (-0.8324) < 0 \\ -0.3202914 < 0$$ and we find that the condition is, indeed, satisfied.

The plot below shows the original series, y, the fit of the ARMA(2,2) model, and 14 out-of-sample forecasts. As can be seen, the out-of-sample forecasts follow a sine-wave pattern.

Keep in mind two things. 1) This is just a very quick analysis (using an automated tool) and a proper treatment would involve following the Box-Jenkins methodology. 2) ARIMA forecasts are good at short-term forecasting, so you may find that long term forecasts from the models in the answers by @David J. Harris and @Glen_b to be more reliable.