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


But the result is really poor.

I would very much appreciate any help.

• 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. Jun 5, 2013 at 18:31
• Have you tried different starting values? Your loss function is non-convex, so different starting values can lead to different solutions. Jun 5, 2013 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. Jun 5, 2013 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. Jun 5, 2013 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.) Jun 5, 2013 at 23:19

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)? Jun 6, 2013 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) Jun 6, 2013 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. Jun 6, 2013 at 10:07 • @Glen_b this method works wonders for my use case. I also need to fit a cos(x) curve, but it doesn't work as well... I changed the reslm to reslm <- lm(y ~ cos(2*pi/per*t)+tan(2*pi/per*t)) but that don't look right. any hints? Feb 12, 2019 at 10:05 • Why do you have a tan term in there? Feb 12, 2019 at 12:03 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. Jun 5, 2013 at 21:08 • @Pascal See my updates above for a recommendation for automatically choosing the starting value for omega. Jun 6, 2013 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. Jun 6, 2013 at 3:51 • I wonder where the x values come into play here. Sure it makes a difference for omega, whether the given y values are separated by 1 or by 5 x steps, doesn't it? – knub Jan 7, 2014 at 15:47 • Programming tip not related to the question: caution when naming R objects as foo.bar. This is due to how R specifies methods for classes. Aug 8, 2019 at 12:43 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. Lastly, hopefully this is a nice addition to some already very informative answers. Reference: Forecasting with dynamic regression models: Alan Pankratz, 1991, (John Wiley and Sons, New York), ISBN 0-471-61528-5 This isn't so different than some of the other solutions but we can simplify the use of nls. y and t and the parameters are as defined in the question. We use: • the plinear algorithm of nls to avoid estimating A and C as plinear only requires initial values for parameters that do not enter linearly. In this case the RHS of the formula should be specified as a matrix with each column multiplying one of the parameters that enter linearly. The parameters themselves are implicit but we can name the columns so that the output is easier to read. • spec.ar to derive an initial value for omega as shown in the code below • 0 for the initial value of phi. This only uses base R -- no packages. spec <- spec.ar(y, plot = FALSE, n.freq = 1000) f <- spec$$freq[which.max(spec$$spec)] fm <- nls(y ~ cbind(C = 1, A = cos(omega*t + phi)), start = list(omega = 2*pi*f, phi = 0), algorithm = "plinear") fm ## Nonlinear regression model ## model: y ~ cbind(C = 1, A = cos(omega * t + phi)) ## data: parent.frame() ## omega phi .lin.C .lin.A ## 0.2636 0.3819 10.1501 2.5838 ## residual sum-of-squares: 22.99 ## ## Number of iterations to convergence: 5 ## Achieved convergence tolerance: 2.346e-06 plot(y ~ t) lines(fitted(fm) ~ t)  The current methods to fit a sin curve to a given data set require a first guess of the parameters, followed by an interative process. This is a non-linear regression problem. A different method consists in transforming the non-linear regression to a linear regression thanks to a convenient integral equation. Then, there is no need for initial guess and no need for iterative process : the fitting is directly obtained. In case of the function y = a + r*sin(w*x+phi) or y=a+b*sin(w*x)+c*cos(w*x), see pages 35-36 of the paper "Régression sinusoidale" published on Scribd : http://www.scribd.com/JJacquelin/documents In case of the function y = a + p*x + r*sin(w*x+phi) : pages 49-51 of the chapter "Mixed linear and sinusoidal regressions". In case of more complicated functions, the general process is explained in the chapter "Generalized sinusoidal regression" pages 54-61, followed by a numerical example y = r*sin(w*x+phi)+(b/x)+c*ln(x), pages 62-63 Another option would be the functions sinusoid and mvrm from package BNSP. data <- data.frame(y, t) model <- y ~ sinusoid(t, harmonics = 2, amplitude = 1, period = 24) m1 <- mvrm(formula = model, data = data, sweeps = 10000, burn = 5000, thin = 2, seed = 1, StorageDir = getwd()) plotOptionsM <- list(geom_point(data = data, aes(x = t, y = y))) plot(x = m1, term = 1, plotOptions = plotOptionsM, intercept = TRUE, quantiles = c(0.025, 0.975), grid = 100)  The argument harmonics specifies the number of sins and cosines to be included. In this case, the model includes 2 pairs of sins and cosines: $$\sin(2 \pi t /per)$$, $$\cos(2\pi t/per)$$, and $$\sin(4 \pi t per)$$, $$\cos(4 \pi t/per)$$. The argument amplitude defines the model for the amplitude. Here amplitude = 1 signifies that the amplitude is constant over t (which I think is appropriate for this case). If you know the lowest and highest point of your cosine-looking data, you can use this simple function to compute all cosine coefficients: 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) }  Below it is used to simulate the variation of temperature throughout the day with a cosine function, by entering the hours and temperature values for the lowest and warmest hour: 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")  The output is below: • This approach would seem to be particularly sensitive to any random-looking departures from pure sinusoidal behavior, which would render it inapplicable to almost any datasets like the one illustrated in the question. Conceivably, it could be used to provide starting values for some of the other iterative approaches suggested in this thread. – whuber Jul 4, 2017 at 20:58 • agree, it is the simplest, would be good for simple approximation under certain assumptions – IVIM Jul 5, 2017 at 21:41 Another option is using the generic function optim or nls. I've tried both none of them is completely robust The following functions takes the data in y and calculates the parameters. 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])  the use is the following: 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)  The following code compares the data ysin <- genSine(as.vector(t), params=reslm1$par)