Non-linear regression problem A physical phenomenon I'm studying obeys the simple law (harmonic oscillator):
$$\theta=\theta_0\cos(\omega t),$$
where $\theta$ is an angle, $\theta_0$ the amplitude of the oscillation, $\omega$ its angular frequency and $t$ is time.
I have a large data set $(\theta, t)$ with which to determine $\theta_0$, $\omega$ and from the latter the period $T$ (from $\omega=\frac{2\pi}{T}$).
My first thought was to linearize with:
$$\arccos\left(\frac{\theta}{\theta_0}\right)=\omega t,$$
but in reality I don't know the value of $\theta_0$, even though I could get a good estimate from the data plot.
Wikipedia has a nice entry about non-linear regression but I'm unfamiliar with the listed methods.
Any help would be appreciated.
 A: Defining the adjusting error as
$$
E(\theta_0,\omega) = \sum_{k=1}^n\left(\theta_k-\theta_0\cos(\omega t_k)\right)^2
$$
the error stationary points are solutions for
$$
\nabla E = \begin{cases}
E_{\theta_0} = \sum_{k=1}^n\left(\theta_k-\theta_0\cos(\omega t_k)\right)\cos(\omega t_k) = 0\\
E_{\omega} = \sum_{k=1}^n\left(\theta_k-\theta_0\cos(\omega t_k)\right)\sin(\omega t_k) = 0
\end{cases}
$$
Those two nonlinear equations should be solved regarding $\theta_0,\omega$ using an iterative procedure like Newton's. This procedure can be described as
$$
(\theta_0,\omega)_{k+1} = (\theta_0,\omega)_k-J_k^{-1}\cdot (E_{\theta_0},E_{\omega})_k
$$
with
$$
J(\theta_0,\omega) = \nabla^2 E(\theta_0,\omega)
$$
Follows a basic python script which implements the procedure

import numpy as np
import random as rd
import matplotlib as mp


n     = 100
sigma = 0.5
theta = 1
omega = 3
tkmax = 10.0
niter = 50
error = 0.0001


def noise(sigma):
    return (rd.random()-0.5)*sigma


thetak = []
tk = list(np.linspace(0.0,tkmax,n))


for k in range(n):
    thetak.append(theta*np.cos(omega*tk[k])+noise(sigma))


mp.pyplot.plot(tk,thetak,'o', color = 'Red')


"""
Data is contained in tk and thetak. This signal is shown in red
"""


def dE(tk,thetak,theta,omega):
    dEtheta = 0.0
    dEomega = 0.0
    for k in range(n):
        wk       = omega*tk[k]
        dEtheta += thetak[k]*np.cos(wk)-0.5*theta*(np.cos(2*wk)+1)
        dEomega += thetak[k]*np.sin(wk)-0.5*theta*np.sin(2*wk)
    return np.array([dEtheta,dEomega])


def d2E(tk,thetak,theta,omega):
    dEthetaomega = 0
    dEtheta2     = 0
    dEomega2     = 0
    for k in range(n):
        wk            = omega*tk[k]
        dEthetaomega += (thetak[k]*np.sin(wk)-theta*np.sin(2*wk))
        dEtheta2     += (np.cos(2*wk)+1)
        dEomega2     += (thetak[k]*np.cos(wk)-theta*np.cos(2*wk))
    a11 = -0.5*dEtheta2
    a12 = omega*dEthetaomega
    a21 = a12
    a22 = omega*dEomega2
    return np.array([[a11,a12],[a21,a22]])


"""
Main iterative procedure
"""


eka = np.array([theta+noise(sigma),omega+noise(sigma)])
print("initial values theta0 = ",eka[0], " omega = ", eka[1])
for j in range(niter):
    invd2E = np.linalg.inv(d2E(tk,thetak,eka[0],eka[1]))
    dEk    = dE(tk,thetak,eka[0],eka[1])
    dek    = np.dot(invd2E,dEk)
    if max(abs(dek)) < error:
        break
    ek  = eka - dek
    eka = ek


print("final values   theta0 = ",eka[0], " omega = ", eka[1])
print("               in ", j," iterations within an error ", error)


theta0 = []
for j in range(n):
    theta0.append(ek[0]*np.cos(ek[1]*tk[j]))


"""
The adjusted function y = eka[0]*cos(eka[1]*tk) is shown in green
"""


mp.pyplot.plot(tk,theta0,'o-',color='Black')
Follow a plot showing in red the data points and in black the adjusting solution.

`
NOTE
Depending on the noise added, for large sigma, strange solutions may appear due to the aliasing phenomena. This problem can be handled by introducing suitable terms into the error definition. 
A: Here's how I would approach the problem:


*

*Find a suitable maximum value of $\theta$ (you might need to filter the data in advance, perhaps using a median filter). You can use a peak detector for this.

*Find the mean of $\theta$. This is best done with a complete number of cycles, but the more cycles you have, the less important this is.

*$\theta_0=\theta_{\max}-\theta_{\text{mean}}.$

*To find $\omega,$ you need to find where your signal $\theta$ crosses $\theta_{\text{mean}}$ on a rising edge. Alternatively, you could reuse the peaks you found before. There are usually library functions in most programming languages to do this. Sort the $t$ values where these crossings or peaks occur, then take the deltas. The median delta is going to be your $T,$ and you can compute $\omega=\frac{2\pi}{T},$ if you really need it.


Here's some LabVIEW code and a front panel that accomplishes this (the code is a VI snippet, so you can download and use the code yourself if you have LabVIEW 2014 SP1 or higher):
Code:

Results:

A: Use a statistics package. For example, R is an environment for statistical computing. In R you can load your data into a 'data frame' (call it mydat) with columns named (say) t and y. The nonlinear model
$$
y = A\cos(\phi + \omega t)\tag1
$$
with amplitude $A$, angular frequency $\omega$ and phase $\phi$ is represented in R as the formula
y ~ A * cos(phi + omega * t)

where A, omega, and phi are the parameters to be fitted and t and y are the columns in your data. To perform a nonlinear regression you invoke the function nls() with the formula, the data, and starting guesses for the parameters. A good guess for $A$ would be half the range of the $y$ values, so you might do:
A0 <- (max(mydat$y) - min(mydat$y))/2
nlfit <- nls(y ~ A * cos(phi + omega * t), data=mydat, start=list(A=A0, omega=1, phi=0))

The estimated parameters are found by inspecting the output object nlfit:
print(summary(nlfit))
ahat <- coef(nlfit)['A']
omegahat <- coef(nlfit)['omega']
phihat <- coef(nlfit)['phi']

You can also use R to plot your data along with the fitted model, and much more.
