Series
Autocovariance and autocorrelation, rolling-window statistics, differencing, and the closed-form conjugate Bayesian updates.
This module covers two small toolkits. The first is for a single time series held as a vector: autocovariance and autocorrelation at one lag or at every lag, the partial autocorrelation at lags one and two, rolling sums, means and variances over a trailing window, and first differences. The second is the four standard conjugate prior updates, where observing data simply adds to the prior’s parameters.
A lag is expressed as a mask rather than a slice. The matrix \([\, j - i = k \,]\) picks out every pair of positions \(k\) apart, and contracting the centred series against it on both sides gives the autocovariance with no indexing. A rolling window is the same idea with a different mask. The masks are n×n (and n×n×n for the whole correlogram at once), so these functions suit series of hundreds to a few thousand points.
Conventions. A series x is a vector of length n, oldest value first. Lags k and window widths w are whole numbers passed as reals. The autocovariance divides by \(n\), not by \(n - k\): this biased estimator is the one correlograms and Box–Jenkins identification use, because it keeps the autocovariance sequence positive-definite. Partial autocorrelation beyond lag two is not provided.
open Prelude.Series
The shapes a lag needs
Three n×n matrices of positions underlie every lag and window: the row index \(i\), the column index \(j\), and their difference \(j - i\).
Building blocks
rowOf(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{rowOf}(x)_{ij} = i\) for positions \(i, j = 0, \dots, n-1\). An n×n matrix holding the row position at every entry.columnOf(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{columnOf}(x)_{ij} = j\) for positions \(i, j = 0, \dots, n-1\). An n×n matrix holding the column position at every entry.gapOf(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{gapOf}(x)_{ij} = j - i\) for positions \(i, j = 0, \dots, n-1\). The n×n matrix of gaps.
Autocovariance and autocorrelation
With \(c_i = x_i - \bar{x}\) the centred series, the autocovariance and autocorrelation at lag \(k\) are
\[\hat{\gamma}(k) = \frac{1}{n} \sum_{j - i = k} c_i\, c_j, \qquad \hat{\rho}(k) = \frac{\hat{\gamma}(k)}{\hat{\gamma}(0)}\]The partial autocorrelation at lag one equals \(\hat{\rho}(1)\); at lag two it is what remains of \(\hat{\rho}(2)\) once lag one has been accounted for:
\[\phi_{22} = \frac{\hat{\rho}(2) - \hat{\rho}(1)^2}{1 - \hat{\rho}(1)^2}\]autocovariance
autocovariance(x, k) · x series vector n · k lag scalar → autocovariance scalar
Let \(x = (x_0, \dots, x_{n-1}) \in \mathbb{R}^n\) be a series, indexed by positions \(0, 1, \dots, n-1\), let \(k \in \{0, 1, \dots, n-1\}\) be a lag, and let \(c_i = x_i - \bar{x}\) with \(\bar{x} = \frac{1}{n} \sum_i x_i\). Then
\[\operatorname{autocovariance}(x, k) = \hat{\gamma}(k) = \frac{1}{n} \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} c_i\, [\, j - i = k \,]\, c_j = \frac{1}{n} \sum_{i=0}^{n-1-k} c_i\, c_{i+k} .\]The sample autocovariance at lag k, divided by n. At lag 0 it is the population variance of the series.
Intuition How much the series moves with itself k steps later: every value is paired with the one k positions ahead, and the products of their deviations from the overall mean are averaged. Dividing by \(n\) rather than by the \(n-k\) pairs shrinks long lags toward zero, but it keeps the sequence over lags a valid (positive semi-definite) covariance function.
Use it when you need one lag’s covariance in the series’ squared units, for example an entry of the Yule–Walker equations for an AR fit, or the variance of a sum of correlated values. For a scale-free number use autocorrelation, and for every lag at once use acov. Each call builds an n×n mask.
time-seriesforecasting
def autocovariance(x, k) ≔ einsum("i,ij,j->", centre(x), lagged(x, k), centre(x)) / count(x)
autocorrelation
autocorrelation(x, k) · x series vector n · k lag scalar → autocorrelation scalar
Let \(x \in \mathbb{R}^n\) be a non-constant series, let \(k \in \{0, 1, \dots, n-1\}\), and let \(\hat{\gamma}\) be the autocovariance of autocovariance. Then
The sample autocorrelation at lag k: the autocovariance at k over the autocovariance at 0. It is 1 at lag 0.
Intuition The autocovariance divided by the variance: the correlation between the series and itself shifted k steps. A value near 1 means that what happened k steps ago strongly predicts the present. A spike at lag 7 in daily data points to a weekly cycle.
Use it when you want to check one specific lag: persistence at lag 1, a suspected seasonal period, or whether a model’s residuals are still correlated (for white noise they stay roughly within \(\pm 2/\sqrt{n}\)). For many lags use acf. The divisor is not guarded, so a constant series divides by zero.
time-seriesseasonalityresidual-diagnostics
def autocorrelation(x, k) ≔ autocovariance(x, k) / autocovariance(x, 0)
acov
acov(x) · x series vector n → autocovariances vector n
Let \(x \in \mathbb{R}^n\) be a series with centred values \(c_i = x_i - \bar{x}\), positions \(0, \dots, n-1\). Then
\[\operatorname{acov}(x)_k = \frac{1}{n} \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} c_i\, [\, j - i = k \,]\, c_j = \hat{\gamma}(k), \qquad k = 0, 1, \dots, n-1 .\]The autocovariance at every lag \(0, 1, \dots, n-1\) at once. It builds an n×n×n mask, so memory grows with \(n^3\).
Intuition autocovariance at every lag \(0, \dots, n-1\) in one contraction, against a single n×n×n mask rather than one n×n mask per call. Entry 0 is the population variance, and later entries fade toward zero because fewer pairs contribute while the divisor stays \(n\).
Use it when you need the whole autocovariance sequence, for example to set up Yule–Walker equations or a spectral estimate. The mask grows as \(n^3\): at \(n = 1000\) it already has a billion entries. For longer series, call autocovariance only for the lags you need.
time-seriesforecasting
def acov(x) ≔ einsum("i,kij,j->k", centre(x), lags(x), centre(x)) / count(x)
acf
acf(x) · x series vector n → autocorrelations vector n
Let \(x \in \mathbb{R}^n\) be a non-constant series. Then
\[\operatorname{acf}(x)_k = \frac{\hat{\gamma}(k)}{\hat{\gamma}(0)} = \hat{\rho}(k), \qquad k = 0, 1, \dots, n-1 .\]The whole correlogram: the autocorrelation at every lag \(0, 1, \dots, n-1\). The first entry is 1. Like acov, it needs \(n^3\) memory.
Intuition The correlogram: autocorrelation against lag. Its shape is a fingerprint of the process. Slow decay suggests a trend or non-stationarity, a sharp cutoff after lag \(q\) suggests a moving-average MA(\(q\)) model, and peaks repeating at a fixed spacing reveal seasonality.
Use it when you explore a new series, choose ARIMA orders (together with pacf1 and pacf2), or check a model’s residuals for leftover structure. High lags rest on few pairs and are shrunk toward 0, so read mainly the first quarter of the lags. It needs \(n^3\) memory like acov. For long series, call autocorrelation at chosen lags.
time-seriesseasonalitymodel-selectionresidual-diagnostics
def acf(x) ≔ acov(x) / autocovariance(x, 0)
pacf1
pacf1(x) · x series vector n → partial autocorrelation scalar
Let \(x \in \mathbb{R}^n\) be a non-constant series with \(n \ge 2\). Then
\[\operatorname{pacf1}(x) = \phi_{11} = \hat{\rho}(1) .\]The partial autocorrelation at lag one, which equals the autocorrelation at lag one.
Intuition Partial autocorrelation asks for the direct link at a lag once the shorter lags have been accounted for. At lag 1 there is nothing in between, so the direct link is just the lag-one autocorrelation.
Use it when you identify an autoregressive model together with pacf2: a clear \(\phi_{11}\) with \(\phi_{22}\) near zero suggests AR(1). It is also the Yule–Walker estimate of an AR(1) coefficient. Its value is identical to autocorrelation(x, 1); the name exists so that a PACF reads uniformly.
time-seriesmodel-selection
def pacf1(x) ≔ autocorrelation(x, 1)
pacf2
pacf2(x) · x series vector n → partial autocorrelation scalar
Let \(x \in \mathbb{R}^n\) be a non-constant series with \(n \ge 3\) and \(\lvert \hat{\rho}(1) \rvert < 1\). Then
\[\operatorname{pacf2}(x) = \phi_{22} = \frac{\hat{\rho}(2) - \hat{\rho}(1)^2}{1 - \hat{\rho}(1)^2} .\]The partial autocorrelation at lag two, in closed form. Deeper lags would need a Yule–Walker solve, which this module does not offer.
Intuition How much lag 2 still tells you once lag 1 is known. In an AR(1) process the lag-2 correlation is \(\rho^2\), carried entirely through lag 1, so \(\phi_{22}\) is 0. A clearly nonzero value means a direct two-step dependence.
Use it when you choose between AR(1) and AR(2), since a partial autocorrelation that cuts off after lag \(p\) points to AR(\(p\)). It equals the second coefficient of a Yule–Walker AR(2) fit. Judge it against roughly \(\pm 2/\sqrt{n}\). Deeper lags are not available in this module.
time-seriesmodel-selectionforecasting
def pacf2(x) ≔ (autocorrelation(x, 2) − autocorrelation(x, 1) · autocorrelation(x, 1)) / (1 − autocorrelation(x, 1) · autocorrelation(x, 1))
Building blocks
lagged(x, k)— Let \(x \in \mathbb{R}^n\) and \(k \in \mathbb{R}\). Then \(\operatorname{lagged}(x, k)_{ij} = [\, j - i = k \,]\). The n×n 0/1 mask of position pairs exactly k apart.lagIndex(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{lagIndex}(x)_{kij} = k\) for \(k, i, j = 0, \dots, n-1\). An n×n×n tensor holding the lag at every entry.gapIndex(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{gapIndex}(x)_{kij} = j - i\) for \(k, i, j = 0, \dots, n-1\). An n×n×n tensor holding the gap at every entry.lags(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{lags}(x)_{kij} = [\, j - i = k \,]\) for \(k, i, j = 0, \dots, n-1\). The n×n×n 0/1 mask of position pairs k apart, for every k at once.
Rolling windows
A trailing window of width \(w\) at position \(i\) covers the positions \(j\) with \(i - w < j \le i\). The first \(w - 1\) windows are shorter than \(w\), and each statistic divides by the number of values actually in the window, so an early value is the mean of what there was rather than a mean padded with zeros. The rolling variance is the population form (divided by the window count):
\[\bar{x}^{(w)}_i = \frac{1}{m_i}\sum_{j \in W_i} x_j, \qquad s^{2\,(w)}_i = \frac{1}{m_i}\sum_{j \in W_i} x_j^2 - \left(\bar{x}^{(w)}_i\right)^2\]rollingSum
rollingSum(x, w) · x series vector n · w window width scalar → sums vector n
Let \(x \in \mathbb{R}^n\) be a series with positions \(0, \dots, n-1\) and let \(w \in \{1, 2, \dots\}\) be the window width. Then
\[\operatorname{rollingSum}(x, w)_i = \sum_{j=0}^{n-1} [\, j \le i \,]\, [\, j > i - w \,]\, x_j = \sum_{j = \max(0,\, i - w + 1)}^{i} x_j .\]The sum of the last w values up to and including each position.
Intuition For each position, the total of the last w values ending there. An n×n mask marks which positions fall in each trailing window, and one contraction sums them all. Early positions sum over the shorter windows that exist, with no padding.
Use it when you aggregate recent activity: requests in the last 60 seconds, sales over the trailing 7 days, or events inside a sliding window for rate limits and alert rules. The first w − 1 entries are partial sums and so smaller by construction: ignore them, or use rollingMean, which divides by the actual window size. Memory grows as \(n^2\).
time-seriesaggregationmonitoring
def rollingSum(x, w) ≔ einsum("ij,j->i", window(x, w), x)
rollingMean
rollingMean(x, w) · x series vector n · w window width scalar → means vector n
Let \(x \in \mathbb{R}^n\) and \(w \in \{1, 2, \dots\}\), let \(W_i = \{ j : i - w < j \le i,\; 0 \le j \le n-1 \}\) be the trailing window at position \(i\), and let \(m_i = \min(w, i + 1)\) be its size. Then
\[\operatorname{rollingMean}(x, w)_i = \bar{x}^{(w)}_i = \frac{1}{m_i} \sum_{j \in W_i} x_j .\]The trailing moving average over w values. The first entries average over the shorter windows available.
Intuition A trailing moving average: each point is replaced by the mean of itself and the w − 1 values before it. Averaging cancels short-term noise. The wider the window, the smoother the line, and the further it lags behind turning points, by about \((w-1)/2\) steps.
Use it when you smooth a noisy series for display, estimate the local level, or build a baseline that new values are compared against (an anomaly is a value far from its rolling mean). It uses only current and past values, so it is safe as a forecasting feature without look-ahead leakage. The first entries average fewer values and are noisier.
smoothingtime-seriesanomaly-detection
def rollingMean(x, w) ≔ rollingSum(x, w) / windowCount(x, w)
rollingSpread
rollingSpread(x, w) · x series vector n · w window width scalar → variances vector n
Let \(x \in \mathbb{R}^n\) and \(w \in \{1, 2, \dots\}\), with trailing windows \(W_i\) of size \(m_i\) as in rollingMean. Then
The trailing moving variance over w values, in population form (divided by the number of values in the window, not by one fewer). The first entry is always 0.
Intuition The variance inside each trailing window, computed as the mean of squares minus the square of the mean. It tracks how volatile the series is now, rather than over its whole history.
Use it when you monitor changing volatility (asset returns, sensor noise), or pair it with rollingMean for a rolling z-score, (x − rollingMean(x, w)) / sqrt(rollingSpread(x, w)). It is the population form, smaller than a sample variance by \((m-1)/m\) for a window of \(m\) values. The first entry is 0, and the result is not clipped at 0, so rounding can leave tiny negatives: guard before a square root or a division.
time-seriesrisk-measurementanomaly-detection
def rollingSpread(x, w) ≔ einsum("ij,j,j->i", window(x, w), x, x) / windowCount(x, w) − rollingMean(x, w) · rollingMean(x, w)
previous
previous(x) · x series vector n → lagged series vector n
Let \(x \in \mathbb{R}^n\) be a series with positions \(0, \dots, n-1\). Then
\[\operatorname{previous}(x)_i = \sum_{j=0}^{n-1} [\, j - i = -1 \,]\, x_j = \begin{cases} x_{i-1}, & i \ge 1, \\ 0, & i = 0 . \end{cases}\]The series shifted by one step: entry i is \(x_{i-1}\), and the first entry is 0.
Intuition The lag operator: the series shifted one step, so entry \(i\) holds the value from step \(i-1\). It is built from the gap mask with no indexing, and position 0 has nothing before it, so it is filled with 0.
Use it when you build lag features for a forecasting or regression model (predict x from previous(x)), or compare each value with the one before it. The first entry is a fill value, not data: drop it before fitting, and note that a ratio such as x / previous(x) divides by zero there. For the change itself use differences.
time-seriesforecastingfeature-engineering
def previous(x) ≔ einsum("ij,j->i", gapOf(x) = 0 − 1, x)
differences
differences(x) · x series vector n → differences vector n
Let \(x \in \mathbb{R}^n\) be a series with positions \(0, \dots, n-1\). Then
\[\operatorname{differences}(x)_i = \bigl(x_i - \operatorname{previous}(x)_i\bigr)\, [\, i > 0 \,] = \begin{cases} x_i - x_{i-1}, & i \ge 1, \\ 0, & i = 0 . \end{cases}\]The first differences \(x_i - x_{i-1}\), the usual step towards stationarity. The result keeps length n; its first entry is 0, not \(x_0\).
Intuition Changes instead of levels: entry \(i\) records how far the series moved since step \(i-1\). Differencing turns a straight-line trend into a constant and a random walk into its independent steps, which is why it is the standard route to a stationary series.
Use it when the series trends or its acf decays slowly, before fitting an ARMA model (the “I” in ARIMA), or whenever the increments are the quantity of interest, such as daily changes. The first entry is 0 by convention, not a real difference: drop it before computing statistics, or it biases them toward zero.
time-seriestrend-removalforecasting
def differences(x) ≔ (x − previous(x)) · (positions(x) > 0)
Building blocks
window(x, w)— Let \(x \in \mathbb{R}^n\) and \(w \in \mathbb{R}\). Then \(\operatorname{window}(x, w)_{ij} = [\, j \le i \,]\, [\, j > i - w \,]\) for positions \(i, j = 0, \dots, n-1\). The n×n 0/1 mask whose row i marks the positions in the trailing window at i.windowCount(x, w)— Let \(x \in \mathbb{R}^n\) and \(w \in \{1, 2, \dots\}\). Then \(\operatorname{windowCount}(x, w)_i = \sum_{j} [\, j \le i \,]\, [\, j > i - w \,] = \min(w, i + 1)\). How many positions each window holds: w, or fewer at the start.
The conjugate updates
Under a conjugate prior the posterior belongs to the same family as the prior, and its parameters are the prior’s parameters plus what the data contributed. No sampling or optimisation is involved.
For \(k\) successes in \(n\) trials, \(n\) Poisson counts \(x\), or a vector of category counts \(c\):
\[\mathrm{Beta}(\alpha, \beta) \to \mathrm{Beta}(\alpha + k,\; \beta + n - k), \qquad \mathrm{Gamma}(\alpha, \beta) \to \mathrm{Gamma}\Bigl(\alpha + \sum_i x_i,\; \beta + n\Bigr), \qquad \mathrm{Dir}(\alpha) \to \mathrm{Dir}(\alpha + c)\]The Gamma prior is in shape–rate form. For a normal mean with known noise standard deviation \(\sigma\) and prior precision \(\tau_0\), the precisions add, \(\tau_n = \tau_0 + n/\sigma^2\), and the posterior mean is the precision-weighted average of the prior mean and the sample mean:
\[\mu_n = \frac{\tau_0\,\mu_0 + (n/\sigma^2)\,\bar{x}}{\tau_0 + n/\sigma^2}\]betaPosteriorA
betaPosteriorA(a, b, k, n) · a prior α scalar · b prior β scalar · k successes scalar · n trials scalar → posterior α scalar
Let \(\alpha, \beta > 0\) (the arguments a and b) be the parameters of a \(\mathrm{Beta}(\alpha, \beta)\) prior, and let \(k\) successes be observed in \(n\) trials, \(0 \le k \le n\). Then
The posterior α of a Beta prior after k successes in n trials. The arguments b and n are accepted for a uniform calling pattern and do not affect the result.
Intuition A Beta(\(\alpha, \beta\)) prior behaves like \(\alpha\) earlier successes and \(\beta\) earlier failures, as pseudo-counts. Observing k successes simply adds them to \(\alpha\). That addition is the whole Bayesian update for a probability under a Beta prior.
Use it when you track a success probability (conversion rate, click-through, defect rate) and need the full posterior, for example to draw samples for Thompson sampling or to compute a credible interval with a Beta quantile function. Pair it with betaPosteriorB; for a single estimate use betaPosteriorMean. Beta(1, 1) is the uniform prior, and updates can be chained batch after batch.
bayesian-updatingrate-estimationa-b-testing
def betaPosteriorA(a, b, k, n) ≔ a + k
betaPosteriorB
betaPosteriorB(a, b, k, n) · a prior α scalar · b prior β scalar · k successes scalar · n trials scalar → posterior β scalar
Let \(\alpha, \beta > 0\) (the arguments a and b) be the parameters of a Beta prior, and let \(k\) successes be observed in \(n\) trials, \(0 \le k \le n\). Then
The posterior β: the prior β plus the number of failures.
Intuition The failure side of the same bookkeeping: \(\beta\) grows by the \(n-k\) trials that failed. With betaPosteriorA it gives the posterior Beta(\(\alpha + k, \beta + n - k\)), and \(\alpha + \beta\) is how many real and pseudo observations the belief now rests on.
Use it when you need the second posterior parameter, for sampling, a credible interval, or the posterior variance \(\alpha\beta / ((\alpha+\beta)^2(\alpha+\beta+1))\) of the updated parameters. In an A/B/n test or a bandit, keep one pair per arm and update each with its own counts. n is the number of trials, not failures, and k must not exceed it.
bayesian-updatinga-b-testingrate-estimation
def betaPosteriorB(a, b, k, n) ≔ b + (n − k)
betaPosteriorMean
betaPosteriorMean(a, b, k, n) · a prior α scalar · b prior β scalar · k successes scalar · n trials scalar → posterior mean scalar
Let \(\alpha, \beta > 0\) (the arguments a and b) be the parameters of a Beta prior, and let \(k\) successes be observed in \(n\) trials, \(0 \le k \le n\). Then
The posterior mean of the success probability.
Intuition Successes plus prior pseudo-successes, over trials plus prior pseudo-trials. With little data it stays near the prior mean \(\alpha/(\alpha+\beta)\); with a lot it approaches the raw rate \(k/n\). This shrinkage is what stops one success in one trial from being read as a 100% rate.
Use it when you rank items by a rate despite very different sample sizes (products by approval, ads by click-through), smooth the rates of small groups, or estimate a conversion rate early in a test. Beta(1, 1) gives Laplace’s rule \((k+1)/(n+2)\). It is a point estimate only; use betaPosteriorA and betaPosteriorB when you need the uncertainty.
rate-estimationbayesian-updatingranking
def betaPosteriorMean(a, b, k, n) ≔ (a + k) / (a + b + n)
normalPosteriorPrecision
normalPosteriorPrecision(t0, s, n) · t0 prior precision of the mean scalar · s known noise standard deviation scalar · n number of observations scalar → posterior precision scalar
Let \(\tau_0 > 0\) (the argument t0) be the prior precision of a normal mean, let \(\sigma > 0\) (the argument s) be the known noise standard deviation, and let \(n \in \mathbb{N}\) be the number of observations. Then
The posterior precision of a normal mean with known noise.
Intuition Precision is one over variance: a measure of how much you know. Each observation with noise \(\sigma\) contributes \(1/\sigma^2\), and independent pieces of information simply add, so the posterior precision is the prior’s plus the data’s.
Use it when you need to know how certain an updated estimate of a mean is, for instance to decide whether more measurements are needed or to fuse readings from sensors of known noise. It depends only on n and s, never on the observed values. s must be the known noise level, not a standard deviation estimated from the same data. Its reciprocal is normalPosteriorSpread.
bayesian-updatingsensor-fusionuncertainty-quantification
def normalPosteriorPrecision(t0, s, n) ≔ t0 + n / (s · s)
normalPosteriorMean
normalPosteriorMean(m0, t0, s, x) · m0 prior mean scalar · t0 prior precision scalar · s known noise standard deviation scalar · x observations vector n → posterior mean scalar
Let \(\mu_0 \in \mathbb{R}\) (the argument m0) be the prior mean, \(\tau_0 > 0\) (the argument t0) the prior precision, \(\sigma > 0\) (the argument s) the known noise standard deviation, and let \(x \in \mathbb{R}^n\) be the observations with mean \(\bar{x} = \frac{1}{n} \sum_i x_i\). Then
The posterior mean: the average of the prior mean and the sample mean, weighted by their precisions.
Intuition A tug of war between the prior mean and the sample mean, each pulling with its precision. With a vague prior (small t0) the result is almost the sample mean; with few, noisy observations it stays near m0.
Use it when you estimate a quantity from noisy measurements plus prior knowledge: calibrating a sensor, shrinking a small group’s average toward a global one, or a one-step Kalman update of a static state. It assumes independent observations with known noise s. t0 is a precision (one over a variance), not a standard deviation, which is an easy slip. For its uncertainty use normalPosteriorSpread.
bayesian-updatingsensor-fusionshrinkage-estimation
def normalPosteriorMean(m0, t0, s, x) ≔ (t0 · m0 + count(x) · mean(x) / (s · s)) / (t0 + count(x) / (s · s))
normalPosteriorSpread
normalPosteriorSpread(t0, s, n) · t0 prior precision scalar · s known noise standard deviation scalar · n number of observations scalar → posterior variance scalar
Let \(\tau_0 > 0\) (the argument t0) be the prior precision, \(\sigma > 0\) (the argument s) the known noise standard deviation, and \(n \in \mathbb{N}\) the number of observations. Then
The posterior variance of the mean, the reciprocal of the posterior precision. Take its square root for a standard deviation.
Intuition The posterior variance of the mean, the reciprocal of the added precisions. Once the data dominate it behaves like \(\sigma^2/n\), so quadrupling the observations halves the posterior standard deviation.
Use it when you put a credible interval on normalPosteriorMean (about ±1.96 times its square root for 95%), or plan how many measurements reach a target uncertainty. Despite the name it returns a variance, so take sqrt for a standard deviation. Like normalPosteriorPrecision, it depends only on t0, s and n, not on the observed values.
uncertainty-quantificationbayesian-updatingsample-size-planning
def normalPosteriorSpread(t0, s, n) ≔ 1 / (t0 + n / (s · s))
gammaPosteriorShape
gammaPosteriorShape(a, x) · a prior shape scalar · x counts vector n → posterior shape scalar
Let \(\alpha > 0\) (the argument a) be the shape of a Gamma prior on a Poisson rate, and let \(x \in \mathbb{N}^n\) be observed counts. Then
The posterior shape of a Gamma prior on a Poisson rate: the prior shape plus the total count.
Intuition A Gamma(\(\alpha, \beta\)) prior on a Poisson rate reads as “\(\alpha\) events seen over \(\beta\) intervals”. New data add the events counted, so the shape grows by the total count.
Use it when you update a belief about an event rate (arrivals per hour, defects per batch, goals per match) and need the full posterior for sampling or intervals. Each entry of x is the count in one interval, and all intervals must have the same length. Pair it with gammaPosteriorRate; for a point estimate use gammaPosteriorMean.
bayesian-updatingrate-estimationcount-data
def gammaPosteriorShape(a, x) ≔ a + total(x)
gammaPosteriorRate
gammaPosteriorRate(b, x) · b prior rate scalar · x counts vector n → posterior rate scalar
Let \(\beta > 0\) (the argument b) be the rate of a Gamma prior on a Poisson rate, and let \(x \in \mathbb{N}^n\) be observed counts. Then
The posterior rate: the prior rate plus the number of observations.
Intuition The exposure side of the update: the rate parameter grows by the number of intervals observed, whatever counts they held. The posterior variance is shape over rate squared, so for a given mean every extra interval narrows the belief.
Use it when you need the second posterior parameter for sampling, an interval, or the posterior variance. It counts the entries of x, which is correct only when every entry covers the same exposure. With unequal exposures (different hours, areas or batch sizes) this form does not apply: the rate should grow by the total exposure instead.
bayesian-updatingrate-estimationcount-data
def gammaPosteriorRate(b, x) ≔ b + count(x)
gammaPosteriorMean
gammaPosteriorMean(a, b, x) · a prior shape scalar · b prior rate scalar · x counts vector n → posterior mean scalar
Let \(\alpha, \beta > 0\) (the arguments a and b) be the shape and rate of a Gamma prior on a Poisson rate, and let \(x \in \mathbb{N}^n\) be observed counts. Then
The posterior mean of the Poisson rate.
Intuition Prior events plus observed events, over prior intervals plus observed intervals: the observed rate shrunk toward the prior rate \(\alpha/\beta\), by an amount that fades as data accumulate.
Use it when you estimate rates for many small units that need stabilising, such as incidents per team, failures per machine or visits per store, so that a unit with one quiet day is not declared to have a rate of zero. It is a point estimate; for the uncertainty use gammaPosteriorShape and gammaPosteriorRate.
rate-estimationbayesian-updatingshrinkage-estimation
def gammaPosteriorMean(a, b, x) ≔ (a + total(x)) / (b + count(x))
dirichletPosterior
dirichletPosterior(a, c) · a prior concentrations vector K · c category counts vector K → posterior concentrations vector K
Let \(\alpha \in (0, \infty)^K\) (the argument a) be the concentrations of a Dirichlet prior and let \(c \in \mathbb{N}^K\) be category counts. Then
The posterior concentrations of a Dirichlet prior after observing counts c: the two vectors added.
Intuition The many-category version of the Beta update: each concentration \(\alpha_j\) is a pseudo-count for category \(j\), and the observed counts are added to it category by category. The total of the concentrations says how firmly the prior is held.
Use it when you track a probability vector over \(K\) categories (word frequencies, market shares, die faces, class priors) and need the full posterior for sampling, or want to chain updates batch after batch. a and c must list the categories in the same order. For the point estimate use dirichletMean.
bayesian-updatingcategorical-datacount-data
def dirichletPosterior(a, c) ≔ a + c
dirichletMean
dirichletMean(a, c) · a prior concentrations vector K · c category counts vector K → posterior mean probabilities vector K
Let \(\alpha \in (0, \infty)^K\) (the argument a) be the concentrations of a Dirichlet prior and let \(c \in \mathbb{N}^K\) be category counts. Then
a point of the simplex \(\Delta^{K-1}\).
The posterior mean of the category probabilities: the posterior concentrations divided by their sum. The result sums to 1.
Intuition Smoothed category shares: counts plus pseudo-counts, normalised to sum to 1. With a all ones this is Laplace (add-one) smoothing, so a category never observed still gets a small, nonzero probability.
Use it when you estimate category probabilities from sparse counts, as in naive Bayes word probabilities, n-gram models or rare-category frequencies, where a zero probability would break a product or a logarithm. Smaller concentrations (0.5, or 0.01) smooth less, larger ones pull the shares toward the prior’s proportions.
categorical-databayesian-updatingsparse-data
def dirichletMean(a, c) ≔ (a + c) / einsum("g->", a + c)
Examples
A correlogram
open Prelude.Series
fn correlogram(x) ↦ acf(x)
fn partialTwo(x) ↦ pacf2(x)
using Examples;
var x = new[] { 1d, 3, 2, 5, 4 };
var acf = TimeSeries.Correlogram(x);
var pacf2 = TimeSeries.PartialTwo(x);
Output (verified):
acf = [1, 0, 0.1, -0.4, -0.2]
pacf2 = 0.1
The centred series is −2, 0, −1, 2, 1, with \(\hat{\gamma}(0) = 10/5 = 2\). At lag 2 the products are \((-2)(-1) + 0 \cdot 2 + (-1)(1) = 1\), so \(\hat{\rho}(2) = (1/5)/2 = 0.1\). Since \(\hat{\rho}(1) = 0\), the lag-two partial autocorrelation equals \(\hat{\rho}(2)\).
A three-point moving average and variance
open Prelude.Series
fn movingAverage(x, w) ↦ rollingMean(x, w)
fn movingVariance(x, w) ↦ rollingSpread(x, w)
using Examples;
var x = new[] { 1d, 3, 2, 5, 4 };
var mean = TimeSeries.MovingAverage(x, 3.0);
var variance = TimeSeries.MovingVariance(x, 3.0);
Output (verified, rounded to 6 decimal places):
roll = [1, 2, 2, 3.333333, 3.666667]
rollvar = [0, 1, 0.666667, 1.555556, 1.555556]
The second entry averages only 1 and 3, because the window there holds two values.
Differencing a series
open Prelude.Series
fn firstDifferences(x) ↦ differences(x)
fn lagged1(x) ↦ previous(x)
using Examples;
var x = new[] { 1d, 3, 2, 5, 4 };
var diff = TimeSeries.FirstDifferences(x);
var prev = TimeSeries.Lagged1(x);
Output (verified):
diff = [0, 2, -1, 3, -1]
prev = [0, 1, 3, 2, 5]
Update a conversion rate with a Beta prior
A Beta(2, 2) prior, then 7 conversions in 10 visits.
open Prelude.Series
fn betaMean(a, b, k, n) ↦ betaPosteriorMean(a, b, k, n)
fn betaA(a, b, k, n) ↦ betaPosteriorA(a, b, k, n)
using Examples;
var mean = TimeSeries.BetaMean(2.0, 2.0, 7.0, 10.0);
var alpha = TimeSeries.BetaA(2.0, 2.0, 7.0, 10.0);
Output (verified, rounded to 6 decimal places):
betaMean = 0.642857
betaA = 9
The posterior is Beta(9, 5), whose mean is \(9/14\).
A normal mean with known noise
Prior mean 0 with precision 1, noise standard deviation 2, and four observations with mean 2.
open Prelude.Series
fn normalMean(m0, t0, s, x) ↦ normalPosteriorMean(m0, t0, s, x)
fn normalVar(t0, s, n) ↦ normalPosteriorSpread(t0, s, n)
using Examples;
var mean = TimeSeries.NormalMean(0.0, 1.0, 2.0, new[] { 1d, 2, 3, 2 });
var variance = TimeSeries.NormalVar(1.0, 2.0, 4.0);
Output (verified):
normalMean = 1
normalVar = 0.5
The data contribute precision \(4/2^2 = 1\), equal to the prior’s, so the posterior mean lies halfway between 0 and 2 and the posterior variance is \(1/(1 + 1)\).
Count data: Gamma–Poisson and Dirichlet–Multinomial
open Prelude.Series
fn gammaMean(a, b, x) ↦ gammaPosteriorMean(a, b, x)
fn dirichletShares(a, c) ↦ dirichletMean(a, c)
using Examples;
var rate = TimeSeries.GammaMean(2.0, 1.0, new[] { 3d, 1, 4, 2 });
var shares = TimeSeries.DirichletShares(new[] { 1d, 1, 1 }, new[] { 5d, 2, 1 });
Output (verified, rounded to 6 decimal places):
gammaMean = 2.4
dir = [0.545455, 0.272727, 0.181818]
The Gamma posterior mean is \((2 + 10)/(1 + 4) = 2.4\); the Dirichlet posterior concentrations are 6, 3, 2 out of 11.