Skip to content

Lecture notes. Written from the material of Berkeley Stat 153, licensed CC BY 4.0. These are notes, not a transcript: the material has been reorganised and rewritten. This adaptation carries the same licence, and the original is linked above.

7. Sinusoidal Nonlinear Regression

What this covers

This chapter asks what changes when a regression model is nonlinear in its parameters, using the running example of a periodic signal of unknown frequency. It works out how a sinusoid's amplitude, frequency and phase relate to the coefficients a computer actually fits, why sampling a continuous signal at discrete times caps which frequencies can be recovered at all, and how the frequency is estimated by least squares despite entering the model nonlinearly. It assumes multiple linear regression — the normal equations, the design matrix, residual sum of squares — from the previous lecture.

From linear to nonlinear regression

Multiple linear regression models \(y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} + w_i\): a relationship that is linear in the parameters \(\beta\). Many real series are not. The lecture's example is the sunspot record: smoothed monthly sunspot counts do not settle near a fixed level but oscillate with a very stable period, and no linear combination of a few predictors captures that on its own.

The same shape of problem — recovering an unknown frequency from noisy data — recurs across fields:

  • Climate data, e.g. El Niño / the Southern Oscillation Index.
  • Astronomy: variable stars, whose brightness oscillates with a period \(1/f\) that is related to the star's intrinsic luminosity, and hence to its distance. Edwin Hubble used exactly this relationship in the 1920s to estimate the distance to the Andromeda galaxy, showing it lay far outside the Milky Way — evidence that the Milky Way is only one galaxy among many.
  • Signal processing: estimating the pitch and formants of a person's voice, part of automatic speech recognition.
  • Economics, "but generally less stable."

The natural model for a periodic signal observed with noise is

\[y_t = \beta_0 + \beta_1 \cos(2\pi f t) + \beta_2 \sin(2\pi f t) + \epsilon_t, \qquad \epsilon_t \overset{\text{iid}}{\sim} N(0,\sigma^2).\]

If the frequency \(f\) were known in advance, this would already be multiple linear regression. Writing

\[y = \begin{pmatrix} y_1 \\ \vdots \\ y_n \end{pmatrix}, \qquad X_f = \begin{pmatrix} 1 & \cos(2\pi f(1)) & \sin(2\pi f(1)) \\ \vdots & \vdots & \vdots \\ 1 & \cos(2\pi f(n)) & \sin(2\pi f(n)) \end{pmatrix}, \qquad \beta = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{pmatrix},\]

the model reads \(y = X_f \beta + \epsilon\) and ordinary least squares applies exactly as before. What makes this a genuinely nonlinear regression problem is that \(f\) is usually not known, and the model is not linear in \(f\): \(\cos(2\pi f t)\) does not change linearly as \(f\) varies. This particular model is called a sinusoid, and it has properties worth pinning down before estimating anything.

What a sinusoid's parameters mean

Write the same curve in amplitude–phase form:

\[s(t) := \beta_0 + R\cos(2\pi f t + \phi).\]
  • \(R\) — the amplitude: the height of the oscillation above (and below) the center line \(\beta_0\).
  • \(f\) — the frequency: cycles per unit time. If \(t\) is in seconds, \(f\) is in Hertz.
  • \(1/f\) — the period: the time to complete one full oscillation.
  • \(\phi\) — the phase: shifts the whole oscillation left or right in time. At \(\phi = 0\) the curve is at its maximum when \(t = 0\).
  • \(\omega := 2\pi f\) — the angular frequency: the rate of change of the angle inside the cosine.
β₀ phase φ R period = 1/f t
The sinusoid $s(t) = \beta_0 + R\cos(2\pi f t + \phi)$: $\beta_0$ sets the center line, $R$ the height of the swings above it, $1/f$ the time between successive peaks, and $\phi$ how far the first peak is shifted away from $t=0$.

This doesn't look like \(y_t = \beta_0 + \beta_1\cos(2\pi f t) + \beta_2 \sin(2\pi f t) + \epsilon_t\) above, but it's the same model. Expand with the identity \(\cos(a+b) = \cos a \cos b - \sin a \sin b\), taking \(a = 2\pi f t\) and \(b = \phi\):

\[s(t) = \beta_0 + R\big(\cos(2\pi f t)\cos\phi - \sin(2\pi f t)\sin\phi\big) = \beta_0 + R\cos\phi \cos(2\pi f t) - R\sin\phi \sin(2\pi f t).\]

Setting \(\beta_1 = R\cos\phi\) and \(\beta_2 = -R\sin\phi\) turns this into exactly \(s(t) = \beta_0 + \beta_1 \cos(2\pi f t) + \beta_2 \sin(2\pi f t)\). The map is invertible, so a fitted \((\hat\beta_1, \hat\beta_2)\) always gives back the amplitude and phase:

\[R = \sqrt{\beta_1^2 + \beta_2^2}, \qquad \phi = \arctan\!\left(-\frac{\beta_2}{\beta_1}\right).\]

So \((\beta_1,\beta_2)\) and \((R,\phi)\) are two parameterizations of the same two degrees of freedom. Regression estimates the first pair directly; the second is read off afterward.

Sampling a continuous signal: the Nyquist limit

The sinusoid above is written in continuous time, but data are collected at discrete sampling times. If the sampling rate is \(f_s\) (samples per second), observations occur at \(t_n = n/f_s\) for \(n = 0,1,2,\dots\), so what is actually fit is

\[y_n = \beta_0 + R\cos\!\left(\frac{2\pi f n}{f_s} + \phi\right) + \epsilon_n.\]

The sampling rate fundamentally constrains which frequencies can be recovered. Frequencies up to half the sampling rate,

\[f_{\max} = \frac{f_s}{2},\]

called the Nyquist frequency, can be reliably estimated. A true component at a higher frequency \(f > f_s/2\) does not vanish from the data — it produces aliasing: the frequency that shows up in the fit is \(|f - n f_s|\) for whichever integer \(n\) brings that value down into range. Two different continuous sinusoids can agree exactly at every sampling time while disagreeing everywhere between them.

samples common to both curves true signal: f = 0.9 aliased estimate: |0.9 − 1| = 0.1 0 1 2 3 4 5 t, in samples (f_s = 1)
With sampling rate $f_s=1$, a true frequency $f=0.9$ (solid) and the aliased frequency $|0.9-1|=0.1$ (dashed) give identical values at every integer sampling time, even though the two curves disagree everywhere between the samples.

Because of this, searching for \(f\) by least squares never needs to test frequencies above the Nyquist limit: it suffices to search \(f \in [0, 1/2]\) once frequency is expressed in units of the sampling rate.

Estimating the frequency by least squares

All five parameters \(\beta_0, \beta_1, \beta_2, f, \sigma\) are fit by minimizing

\[S(\beta_0,\beta_1,\beta_2,f,\sigma) := \sum_{t=1}^n \Big(y_t - \beta_0 - \beta_1\cos(2\pi f t) - \beta_2 \sin(2\pi f t)\Big)^2.\]

Minimizing over all five at once is a genuine nonlinear least squares problem. But for any fixed value of \(f\), this is ordinary linear regression on the design matrix \(X_f\) from the second section — \(f\) enters only through the columns of \(X_f\), not through \(\beta_0,\beta_1,\beta_2\). That lets the frequency be searched for separately from the linear coefficients:

  1. Take a grid of candidate values of \(f\) in \([0, 1/2]\).
  2. For each \(f\) on the grid, build \(X_f\), regress \(y\) on \(X_f\), and record the residual sum of squares \(\mathrm{RSS}(f)\).
  3. Take \(\hat f\) to be the grid value that minimizes \(\mathrm{RSS}(f)\).
  4. With \(f\) fixed at \(\hat f\), take \(\hat\beta\) and \(\hat\sigma\) to be the usual least-squares estimates from regressing \(y\) on \(X_{\hat f}\).

Sources

← 6. Simple and Multiple Linear Regression · Contents · 8. Improving upon the linear model →