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.

34. Inference in Sinusoid Models¶

What this covers¶

This chapter works through Code Lab 3 of the course: data are generated by a single sinusoid of unknown frequency buried in noise, and the question is how to estimate that frequency, attach an uncertainty statement to the estimate, and understand why some frequencies are indistinguishable from others once you sample at equally spaced times. It assumes ordinary least squares and its MLE interpretation, the basic mechanics of Bayesian inference (prior times likelihood gives a posterior, up to a normalizing constant), and the $\chi^2$ and multivariate normal distributions used when sampling from a Bayesian linear model.

The sinusoid regression model¶

The model is $$y_t = \beta_0 + \beta_1\cos(2\pi f t) + \beta_2 \sin(2\pi f t) + \epsilon_t, \qquad t=1,\dots,n,$$ with $\epsilon_t$ iid noise and $\beta_0,\beta_1,\beta_2,f,\sigma$ all unknown. The lab works with a simulated instance: $n=400$, true $f=0.2035$, $\beta=(0,3,5)$, $\sigma=10$.

The structural fact the rest of the lab exploits: for a fixed value of $f$, the design matrix $X_f = [\mathbf 1,\ \cos(2\pi f t),\ \sin(2\pi f t)]$ is completely determined, and the model is an ordinary linear regression in $\beta$. The only genuinely nonlinear parameter is $f$, since it sits inside the cosine and sine. That turns a nonlinear least-squares problem in four parameters into a one-dimensional search over $f$, with $\beta$ profiled out by OLS at each candidate value.

Estimating the frequency by profile least squares¶

Define the profile residual sum of squares $$RSS(f) = \min_{\beta_0,\beta_1,\beta_2} \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,$$ computed by fitting an OLS regression for each candidate $f$ and summing squared residuals:

def rss(f):
    xcos = np.cos(2*np.pi*f*t)
    xsin = np.sin(2*np.pi*f*t)
    X = np.column_stack([np.ones(n), xcos, xsin])
    md = sm.OLS(y, X).fit()
    return np.sum(md.resid**2)

$RSS(f)$ is not convex in $f$ — it oscillates, with local minima at spurious frequencies that happen to fit part of the sinusoid — so the lab does not use gradient-based optimization. Instead it evaluates $RSS(f)$ on a dense grid of $100{,}000$ points over $f\in[0,\tfrac12]$ (the restriction to $[0,\tfrac12]$ is explained below) and takes $\hat f = \arg\min_f RSS(f)$. On the simulated data this gives $\hat f = 0.20355$, matching the generating value $f = 0.2035$ to four digits.

Estimating the remaining parameters¶

Once $f$ is fixed at $\hat f$, the model is a genuine linear regression, and the remaining parameters are estimated exactly as in ordinary least squares:

$\hat\beta$ is the OLS coefficient vector of $y$ on $X_{\hat f}$ — already computed inside rss;
the MLE of $\sigma^2$ is $\widehat{\sigma^2}_{\mathrm{MLE}} = RSS(\hat f)/n$ (biased downward, as always for OLS);
the unbiased estimator divides by the residual degrees of freedom instead, $\widehat{\sigma^2} = RSS(\hat f)/(n-p)$ with $p=3$ parameters in $\beta$.

On the simulated data: $\hat\beta = (-0.053,\ 3.004,\ 5.639)$ against the true $(0,3,5)$, and $\hat\sigma_{\mathrm{MLE}}=9.50$, $\hat\sigma_{\mathrm{unbiased}}=9.54$ against the true $\sigma=10$ — recovered reasonably well from a single noisy realization.

A Bayesian posterior for the frequency¶

Point estimates say nothing about how confident to be in $\hat f$. The lab quantifies that with a posterior for $f$: $$p(f\mid \text{data}) \ \propto\ \mathbb 1\{0\le f\le \tfrac12\}\ \big|X_f^\top X_f\big|^{-1/2} \left(\frac{1}{RSS(f)}\right)^{(n-p)/2}.$$ This is the marginal posterior of $f$ once $\beta$ and $\sigma$ have been integrated out under a flat/reference prior; that derivation belongs to the main lecture on Bayesian linear models and is not repeated in the lab — here the formula is the starting point for computation.

Working on the log scale¶

With $n=400$, the exponent $(n-p)/2\approx 198$, so $RSS(f)^{-(n-p)/2}$ overflows or underflows in floating point. The lab instead computes $$\log p(f\mid\text{data}) = \frac{p-n}{2}\log RSS(f) \;-\; \tfrac12 \log\big|X_f^\top X_f\big| + \text{const},$$ using np.linalg.slogdet for a numerically safe log-determinant. There is a second hazard at the boundary: as $f\to 0$ or $f\to\tfrac12$, the columns of $X_f$ become collinear (the sine column vanishes as $f\to 0$, for instance), so $X_f^\top X_f$ becomes near-singular and $|X_f^\top X_f|^{-1/2}$ blows up. The grid used for the posterior avoids the immediate neighborhood of $0$ and $\tfrac12$ — the lab zooms in on $f\in[0.05,0.35]$, informed by the earlier RSS-versus-frequency curve (the lab calls this "the periodogram") showing the maximizing frequency near $0.2$.

The log posterior traces a curve similar in shape to that periodogram — both peak at the same frequencies — but once exponentiated the posterior has essentially one dominant peak: the exponent $(n-p)/2$ that separates the true frequency's tiny $RSS$ from every competing frequency's larger $RSS$ is large enough, at $n=400$, to swamp every secondary peak.

Credible intervals from a grid¶

With the posterior evaluated and normalized to sum to 1 on an equally spaced grid, a $95\%$ credible interval is built by growing a symmetric window around the posterior mode until it captures $95\%$ of the mass:

def PostProbAroundMax(m):
    est_ind = np.argmax(postvals)
    return np.sum(postvals[(est_ind-m):(est_ind+m)])

m = 0
while PostProbAroundMax(m) <= 0.95:
    m += 1

Since the posterior is unimodal away from the boundary singularities, a symmetric window around the mode stands in for a highest-posterior-density interval. On the simulated data this gives $m=10$ grid steps either side of the mode: $\hat f = 0.20356$ with $95\%$ credible interval $(0.20326,\ 0.20386)$ — matching the posterior mean ($0.20355$), the posterior mode ($0.20356$), and the least-squares estimate ($0.20355$) to several digits.

Posterior sampling for the whole parameter vector¶

To visualize uncertainty in the fitted curve, not just in $f$, the lab draws joint posterior samples of $(f,\sigma,\beta)$ by composition — draw $f$, then $\sigma$ given $f$, then $\beta$ given $(f,\sigma)$:

Sample $f$ from the discretized posterior: rng.choice(allfvals, N, p=postvals), treating the normalized grid values as a discrete approximation to the continuous posterior.
Sample $\sigma$ given $f$, using the pivotal fact (proved in Problem 5 of Homework 1, not re-derived here) that $$\frac{RSS(f)}{\sigma^2}\ \Big|\ \text{data}, f \ \sim\ \chi^2_{n-p},$$ so a draw $\chi\sim\chi^2_{n-p}$ gives the posterior draw $\sigma = \sqrt{RSS(f)/\chi}$.
Sample $\beta$ given $(f,\sigma)$ from $$\beta \mid \text{data}, \sigma, f \ \sim\ N_3\!\big(\hat\beta_f,\ \sigma^2 (X_f^\top X_f)^{-1}\big),$$ the same covariance shape as the sampling distribution of the OLS estimator, used here as a conditional posterior.

Repeating this $N=2000$ times and overlaying the resulting fitted curves $\hat\beta_0+\hat\beta_1\cos(2\pi f t)+\hat\beta_2\sin(2\pi f t)$ on the raw data gives a band of plausible curves — a direct visualization of parameter uncertainty rather than a single best fit. Summarizing the 2000 draws recovers the generating parameters well: posterior means $(\hat f,\hat\beta_0,\hat\beta_1,\hat\beta_2,\hat\sigma) \approx (0.2036,\ -0.048,\ 2.93,\ 5.49,\ 9.57)$ against the true $(0.2035,\ 0,\ 3,\ 5,\ 10)$.

Frequency aliasing¶

Now suppose the data are generated instead by $f = 2-0.2035$, still sampled at the equally spaced integer times $t=1,\dots,n$. The frequency estimate produced by the same procedure is still $\approx 0.2035$ — the two frequencies are aliases of one another. For integer $t$, $$\cos\big(2\pi(2-f)t\big) = \cos(4\pi t - 2\pi f t) = \cos(-2\pi f t) = \cos(2\pi f t),$$ $$\sin\big(2\pi(2-f)t\big) = \sin(4\pi t - 2\pi f t) = \sin(-2\pi f t) = -\sin(2\pi f t),$$ using that $4\pi t$ is an integer multiple of $2\pi$, and that $\cos$ and $\sin$ are $2\pi$-periodic. So a sinusoid at frequency $2-f$ agrees exactly with one at frequency $f$ on every integer sample, with the sine term flipping sign. No amount of data collected at equally spaced integer times can distinguish $f$ from $2-f$; correspondingly, fitting the model at the "wrong" alias recovers the true $\beta_0,\beta_1$ but a sign-flipped estimate of $\beta_2$, exactly as the identity predicts.

cos(2π·0.2·t) (solid) and cos(2π·1.8·t) (dashed) are different continuous curves but coincide at every integer t — the dots mark where equally spaced sampling cannot tell f = 0.2 apart from its alias 2 − 0.2 = 1.8.

Breaking aliasing with unequally spaced samples¶

The aliasing identity uses $t$ integer in an essential way — $4\pi t$ is a multiple of $2\pi$ only then. If the sampling times are not equally spaced, that cancellation no longer happens, and frequencies outside $[0,\tfrac12]$ become identifiable again. The lab demonstrates this by regenerating the data with $f = 2-0.2035\approx 1.80$, but with $t$ drawn as $n$ sorted uniform points on $[1,n]$ rather than the integers $1,\dots,n$. Repeating the same profile-RSS search, now over the wider grid $f\in[0,10]$ — there is no longer a reason to restrict to $[0,\tfrac12]$ — gives $\hat f = 1.7965$, matching the true value directly with no aliasing. The Bayesian posterior (computed the same way, again on the log scale and again excluding a neighborhood of $f=0$ to avoid the same near-singularity as before) concentrates tightly around $1.7965$ as well, with posterior mean $1.79654$ and mode $1.79652$.

Sources¶

Berkeley STAT 153, Fall 2025, Code Lab 3 — Inference in Sinusoid Models (Jupyter notebook), converted at docs/statistics/berkeley/stat153/fall-2025/CodeLabThree153248Fall2025.md. All code, numerical output, model statements and the aliasing discussion in this chapter are drawn from this single source; there is no separate slide deck or lecture transcript supplied for this session, and no exercises were supplied.
The lab refers to, but does not itself contain: (i) the derivation of the posterior density $p(f\mid\text{data})$ from a flat/reference prior by integrating $\beta$ and $\sigma$ out of the likelihood, attributed here to "the main lecture on Bayesian linear models"; (ii) a periodogram plotted earlier in the course, treated in the lab as equivalent in shape to the RSS-versus-frequency curve; (iii) "Problem 5 of Homework One," which proves the pivotal result $RSS(f)/\sigma^2 \mid \text{data}, f \sim \chi^2_{n-p}$ used when sampling $\sigma$ from its posterior. None of these three is reproduced here, since none is part of the supplied material.

← 33. Scaling and Initialization in PyTorch · Contents · 35. Estimating MA(1) Parameters →