Non-linear regression problem

Question

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.

Answer 1

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.

Answer 2

Here's how I would approach the problem:

1. 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.
2. 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.
3. $\theta_0=\theta_{\max}-\theta_{\text{mean}}.$
4. 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:

Answer 3

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.