Fitting
A two-component normal mixture by expectation–maximisation, a reproducible Metropolis chain, AR(1) and AR(2) fits with forecasts, and exponential smoothing.
This module fits models that need an iteration rather than a closed form: a two-component normal mixture by expectation–maximisation (EM), the mean and spread of a normal target estimated by a Metropolis chain, autoregressive models of order one and two by the Yule–Walker equations, and simple exponential smoothing. Each iterative method is a definition that names itself, which Hilbert compiles to a bounded loop.
Randomness is always an input. The Metropolis chain does not draw random numbers itself: you pass it vectors of uniform draws in [0, 1), and it consumes one of each per turn. The same draws always give the same chain, so a Monte Carlo result can be reproduced exactly. In C#, MLambda.Hilbert.Runtime.Random produces seeded draws (new MLambda.Hilbert.Runtime.Random(seed).Uniform(n)).
Conventions. Every input series or sample x is a vector n; positions are counted from 0. The mixture functions each run the full EM fit (80 iterations), so calling several of them repeats the work; the results agree because the fit is deterministic. EM finds a local maximum of the likelihood, and it starts from the lower and upper quartiles of the data, so “first component” normally means the one with the smaller mean.
Honest limitations. Autoregressive fits stop at order two; higher orders need a Yule–Walker solve on a leading block of the autocorrelation matrix, which the language cannot slice. The Metropolis chain targets a normal density only, keeps no history (it returns the mean and spread of its states, not the states) and has no burn-in, so start it near the bulk of the target. The chain is slow, and its cost grows faster than the number of draws: every turn reads its draws by a contraction over the whole draw vector. In our runs a walk of 500 draws took under five seconds and a walk of 2 000 draws about four minutes (chainSpread makes three walks).
open Prelude.Fitting
Reading one element of a carried vector
einsum has no notion of position, so “the k-th element” is a contraction against a one-hot mask built from the vector of positions \(0, 1, \dots, n-1\). Loops build that position vector once and carry it, because building it costs \(O(n^2)\).
elementAt
elementAt(v, k) · v vector n · k position, from 0 scalar → element scalar
Let \(v = (v_0, v_1, \dots, v_{n-1}) \in \mathbb{R}^n\) be a vector indexed by positions \(0, 1, \dots, n-1\), and let \(k \in \mathbb{R}\). Then
\[\operatorname{elementAt}(v, k) = \sum_{i=0}^{n-1} v_i\, [\,i = k\,],\]which is \(v_k\) when \(k \in \{0, 1, \dots, n-1\}\) and \(0\) otherwise.
The element of v at position k, counting from 0.
Intuition Hilbert has no indexing, so “the k-th element” is a dot product with a vector that is 1 at position k and 0 elsewhere. The mask selects exactly one entry; an out-of-range k selects nothing and gives 0 rather than an error.
Use it when a formula needs a specific observation, such as the second-to-last value in ar2Forecast. Each call builds the position vector, which costs \(O(n^2)\): fine once, ruinous inside a loop. In a loop, build positions(v) once, carry it, and call pickAt as the chain and smoothing loops do. For the final element, lastOf reads more clearly.
element-accessnumerical-performance
def elementAt(v, k) ≔ pickAt(v, positions(v), k)
lastOf
lastOf(v) · v vector n → element scalar
Let \(v = (v_0, \dots, v_{n-1}) \in \mathbb{R}^n\) with \(n \ge 1\). Then
\[\operatorname{lastOf}(v) = \sum_{i=0}^{n-1} v_i\, [\,i = n - 1\,] = v_{n-1} .\]The last element of v.
Intuition The most recent observation of a series, reached with the same one-hot mask as elementAt at position \(n-1\). For forecasting, it is the point the next value grows from.
Use it when a one-step forecast starts from the latest value, as ar1Forecast and ar2Forecast do, or when you report the current state of a series. It carries the same \(O(n^2)\) cost as elementAt, so do not call it on every turn of a loop. An empty vector has no last element.
element-accessforecasting
def lastOf(v) ≔ elementAt(v, count(v) − 1)
Building blocks
pickAt(v, ix, k)— Let \(v \in \mathbb{R}^n\), \(\mathit{ix} \in \mathbb{R}^n\) and \(k \in \mathbb{R}\). Then \(\operatorname{pickAt}(v, \mathit{ix}, k) = \sum_{i} v_i\,[\,\mathit{ix}_i = k\,]\); with \(\mathit{ix} = (0, 1, \dots, n-1)\) this is \(v_k\). The element ofvat positionk, given a precomputed position vectorix; what loops call.
Expectation–maximisation, over a two-component normal mixture
The model is \(x \sim (1 - w)\,\mathcal N(\mu_1, \sigma_1^2) + w\,\mathcal N(\mu_2, \sigma_2^2)\). The E-step computes the responsibility of the second component for each point,
\[r_i = \frac{w\,\varphi(x_i; \mu_2, \sigma_2)}{(1 - w)\,\varphi(x_i; \mu_1, \sigma_1) + w\,\varphi(x_i; \mu_2, \sigma_2)}\]and the M-step re-estimates the parameters as weighted averages: \(w = \bar r\), \(\mu_2 = \sum_i r_i x_i / \sum_i r_i\), \(\sigma_2^2 = \sum_i r_i (x_i - \mu_2)^2 / \sum_i r_i\), and the same with \(1 - r_i\) for the first component. The fit starts from \(\mu_1 = Q_1\), \(\mu_2 = Q_3\), both spreads equal to the sample standard deviation and \(w = 0.5\), and runs 80 iterations. Spreads are floored at \(10^{-6}\).
On this page \(\varphi(x; \mu, \sigma) = \exp\!\big(-(x - \mu)^2 / (2\sigma^2)\big) / (\sigma\sqrt{2\pi})\) is the normal density of Prelude.Distributions, applied entrywise, and \(\gamma(v) = v + 10^{-300}\,[\,\lvert v \rvert < 10^{-300}\,]\) is its divisor guard.
responsibilities
responsibilities(x) · x sample n → responsibilities n
Let \(x = (x_1, \dots, x_n) \in \mathbb{R}^n\) be a sample, \(n \ge 2\). For parameters \(\theta = (\mu_1, \mu_2, \sigma_1, \sigma_2, w)\) let the E-step be
\[E(\theta)_i = \frac{w\,\varphi(x_i; \mu_2, \sigma_2)}{\gamma\big((1 - w)\,\varphi(x_i; \mu_1, \sigma_1) + w\,\varphi(x_i; \mu_2, \sigma_2)\big)}, \qquad i = 1, \dots, n,\]and for weights \(\rho \in \mathbb{R}^n\) let \(c(\rho) = \sum_i \rho_i x_i / \gamma(\sum_i \rho_i)\) and \(s(\rho) = \sqrt{\max\big(\sum_i \rho_i (x_i - c(\rho))^2 / \gamma(\sum_i \rho_i),\; 10^{-12}\big)}\). Let the M-step be
\[M(r) = \Big(c(\mathbf{1}_n - r),\; c(r),\; s(\mathbf{1}_n - r),\; s(r),\; \frac{1}{n}\sum_{i=1}^{n} r_i\Big).\]Let \(\theta^{(0)} = (Q_{1/4}, Q_{3/4}, s_x, s_x, \tfrac12)\), where \(s_x = \sqrt{\frac{1}{n-1}\sum_i (x_i - \bar x)^2}\) and \(Q_q = x_{(\lfloor h \rfloor)}\,(1 - (h - \lfloor h \rfloor)) + x_{(\lceil h \rceil)}\,(h - \lfloor h \rfloor)\) with \(h = q(n - 1)\) and order statistics \(x_{(0)} \le \dots \le x_{(n-1)}\); and let \(\theta^{(j+1)} = M(E(\theta^{(j)}))\). Then
\[\operatorname{responsibilities}(x) = E(\theta^{(80)}) .\]The fitted probability that each point belongs to the second component, after 80 EM iterations. The first component’s responsibilities are one minus these.
Intuition Soft cluster labels. EM alternates between two easy questions: given the current bell curves, how likely is each point to belong to each one (the E-step), and given those soft memberships, where should each curve sit and how wide should it be (the M-step)? Each round never lowers the likelihood.
Use it when you suspect data come from two overlapping populations, such as two machine modes, healthy versus faulty readings, or a bimodal histogram, and want each point’s membership probability. Threshold at 0.5 for hard labels. It runs a fixed 80 iterations from the quartiles and finds a local optimum; heavily unbalanced mixtures may start badly. It fits exactly two components, and a single normal population still gets split in two, so check the result makes sense.
clusteringmixture-modellinganomaly-detection
def responsibilities(x) ≔ emLoop(80, x, q1(x), q3(x), std(x), std(x), 0.5)
mixtureWeight
mixtureWeight(x) · x sample n → weight scalar
Let \(x \in \mathbb{R}^n\) be a sample and \(r = \operatorname{responsibilities}(x)\). Then
\[\operatorname{mixtureWeight}(x) = \frac{1}{n}\sum_{i=1}^{n} r_i .\]The mixing weight \(w\) of the second component; the first component’s weight is \(1 - w\).
Intuition The share of the population that belongs to the second component: the soft memberships added up and divided by n. It counts each point fractionally, so a point halfway between the groups adds half to each.
Use it when you estimate prevalence: the fraction of faulty units, of a second species, of a high-activity regime. A weight near 0 or 1 suggests one component is barely used and a single normal may describe the data as well. Every mixture function reruns the full 80-iteration fit, so if you need several parameters, compute responsibilities once and derive them from it.
mixture-modellingrate-estimation
def mixtureWeight(x) ≔ mean(responsibilities(x))
mixtureFirstMean
mixtureFirstMean(x) · x sample n → mean scalar
Let \(x \in \mathbb{R}^n\) be a sample and \(r = \operatorname{responsibilities}(x)\). Then
\[\operatorname{mixtureFirstMean}(x) = \frac{\sum_i (1 - r_i)\, x_i}{\gamma\big(\sum_i (1 - r_i)\big)} .\]The mean \(\mu_1\) of the first component.
Intuition An average of all points in which each counts by how strongly it belongs to the first component. Points clearly in the other group contribute almost nothing, so the mean is not dragged toward them the way a plain mean would be.
Use it when you want the centre of the lower population, for example the normal operating level beneath occasional high readings. “First” is set by the starting point, the lower quartile, so it is normally the component with the smaller mean, but EM does not enforce that order: compare with mixtureSecondMean before interpreting. It reruns the full fit on every call.
mixture-modellingparameter-estimation
def mixtureFirstMean(x) ≔ centreOf(1 − responsibilities(x), x)
mixtureSecondMean
mixtureSecondMean(x) · x sample n → mean scalar
Let \(x \in \mathbb{R}^n\) be a sample and \(r = \operatorname{responsibilities}(x)\). Then
\[\operatorname{mixtureSecondMean}(x) = \frac{\sum_i r_i\, x_i}{\gamma\big(\sum_i r_i\big)} .\]The mean \(\mu_2\) of the second component.
Intuition The responsibility-weighted average for the second component, the one that started at the upper quartile. The gap between the two component means, relative to their spreads, says how clearly the data split into two groups.
Use it when you need the centre of the upper population, such as the elevated level of a faulty or high-activity regime. If the two means come out nearly equal, the data do not support two separate groups. Like its siblings it recomputes the whole 80-iteration fit.
mixture-modellingparameter-estimation
def mixtureSecondMean(x) ≔ centreOf(responsibilities(x), x)
mixtureFirstSpread
mixtureFirstSpread(x) · x sample n → standard deviation scalar
Let \(x \in \mathbb{R}^n\) be a sample, \(r = \operatorname{responsibilities}(x)\), \(\rho = \mathbf{1}_n - r\) and \(\mu_1 = \operatorname{mixtureFirstMean}(x)\). Then
\[\operatorname{mixtureFirstSpread}(x) = \sqrt{\max\!\left(\frac{\sum_i \rho_i (x_i - \mu_1)^2}{\gamma\big(\sum_i \rho_i\big)},\; 10^{-12}\right)} .\]The standard deviation \(\sigma_1\) of the first component.
Intuition How wide the first bell curve is, measured only over the points it claims, each weighted by its membership. It is floored at \(10^{-6}\) so a component that collapses onto a single point cannot drive the density to infinity.
Use it when you need the within-group variability of the first population, for example to set a control limit (“more than 3 spreads above the first mean is unusual”) or to judge overlap with the second component. A value near the floor means the component has shrunk onto one or a few identical points, a known EM failure; treat that fit with suspicion.
mixture-modellingparameter-estimationnumerical-stability
def mixtureFirstSpread(x) ≔ componentSpread(1 − responsibilities(x), x)
mixtureSecondSpread
mixtureSecondSpread(x) · x sample n → standard deviation scalar
Let \(x \in \mathbb{R}^n\) be a sample, \(r = \operatorname{responsibilities}(x)\) and \(\mu_2 = \operatorname{mixtureSecondMean}(x)\). Then
\[\operatorname{mixtureSecondSpread}(x) = \sqrt{\max\!\left(\frac{\sum_i r_i (x_i - \mu_2)^2}{\gamma\big(\sum_i r_i\big)},\; 10^{-12}\right)} .\]The standard deviation \(\sigma_2\) of the second component.
Intuition The width of the second bell curve, weighted by the second component’s memberships. The two spreads are fitted separately, so one population may be tight and the other diffuse.
Use it when you describe the second population’s variability or check separation: when the gap between the means is small compared with the spreads, the groups overlap heavily and individual memberships from responsibilities will be uncertain. The same \(10^{-6}\) floor and collapse warning apply as for mixtureFirstSpread.
mixture-modellingparameter-estimation
def mixtureSecondSpread(x) ≔ componentSpread(responsibilities(x), x)
Building blocks
weighed(w, m, s, x)— Let \(w, m \in \mathbb{R}\), \(s > 0\) and \(x \in \mathbb{R}^n\). Then \(\operatorname{weighed}(w, m, s, x)_i = w\,\varphi(x_i; m, s)\). A normal density atxscaled by the weightw.blend(m1, m2, s1, s2, w, x)— Let \(\mu_1, \mu_2 \in \mathbb{R}\), \(\sigma_1, \sigma_2 > 0\), \(0 \le w \le 1\) and \(x \in \mathbb{R}^n\). Then \(\operatorname{blend}(\mu_1, \mu_2, \sigma_1, \sigma_2, w, x)_i = (1 - w)\,\varphi(x_i; \mu_1, \sigma_1) + w\,\varphi(x_i; \mu_2, \sigma_2)\). The mixture density at every point.owed(m1, m2, s1, s2, w, x)— With the same symbols, \(\operatorname{owed}(\mu_1, \mu_2, \sigma_1, \sigma_2, w, x)_i = w\,\varphi(x_i; \mu_2, \sigma_2) / \gamma(\operatorname{blend}(\mu_1, \mu_2, \sigma_1, \sigma_2, w, x)_i)\). The E-step: the second component’s responsibility for every point.mass(r)— Let \(r \in \mathbb{R}^n\). Then \(\operatorname{mass}(r) = \sum_{i=1}^{n} r_i\). The sum of a responsibility vector.centreOf(r, x)— Let \(r, x \in \mathbb{R}^n\). Then \(\operatorname{centreOf}(r, x) = \sum_i r_i x_i / \gamma(\sum_i r_i)\). The responsibility-weighted mean ofx.componentSpread(r, x)— Let \(r, x \in \mathbb{R}^n\) and \(c = \operatorname{centreOf}(r, x)\). Then \(\operatorname{componentSpread}(r, x) = \sqrt{\max(\sum_i r_i (x_i - c)^2 / \gamma(\sum_i r_i),\, 10^{-12})}\). The responsibility-weighted standard deviation ofx, floored at 1e-6.emLoop(k, x, m1, m2, s1, s2, w)— Let \(k \in \mathbb{N}\), \(x \in \mathbb{R}^n\), \(\theta = (\mu_1, \mu_2, \sigma_1, \sigma_2, w)\), and \(E\), \(M\) as inresponsibilities. Then \(\operatorname{emLoop}(0, x, \theta) = E(\theta)\) and \(\operatorname{emLoop}(k, x, \theta) = \operatorname{emLoop}(k - 1, x, M(E(\theta)))\). RunskEM iterations and returns the responsibilities of the final parameters.
A Metropolis chain
The chain samples the normal density \(\mathcal N(m, s^2)\) with a uniform random-walk proposal. At each turn it reads one proposal draw \(p\) and one acceptance draw \(u\), both uniform on [0, 1):
\[y = x + \text{step}\cdot(2p - 1), \qquad x \leftarrow \begin{cases} y & u < \varphi(y; m, s)/\varphi(x; m, s) \\ x & \text{otherwise} \end{cases}\]The acceptance is computed as a 0/1 mask rather than a branch. The number of turns is the length of the draw vectors, which must be equal. The draws are consumed from the last position to the first: turn \(j\) reads the draws at position \(n - j\). The chain carries only the running sum (or sum of squares) of its states, not the states themselves.
chainMean
chainMean(m, s, u, p, step, x0) · m target mean scalar · s target standard deviation scalar · u acceptance draws n · p proposal draws n · step half-width of the proposal scalar · x0 starting state scalar → mean scalar
Let \(m \in \mathbb{R}\) and \(s > 0\) be the target mean and standard deviation, let \(u, p \in [0, 1)^n\) be independent uniform draws at positions \(0, 1, \dots, n-1\), let \(\delta > 0\) (the argument step) be the half-width of the proposal, and let \(x^{(0)} \in \mathbb{R}\) (the argument x0) be the starting state. For \(j = 1, \dots, n\) let \(q = n - j\),
Then
\[\operatorname{chainMean}(m, s, u, p, \delta, x^{(0)}) = \frac{1}{n}\sum_{j=1}^{n} x^{(j)} .\]The average of the chain’s states over its n turns (the starting state is not included): a Monte Carlo estimate of m.
Intuition A random walker that always accepts a step uphill in density and accepts a downhill step only sometimes, in proportion to how much lower it goes. Over time it spends time in each region in proportion to the target density, so averaging its positions estimates the target’s mean.
Use it when you learn, test or demonstrate MCMC: the target is a normal whose mean you already know, so the result checks sampler behaviour, step size and chain length against the true answer, reproducibly from a seed. It does not sample an arbitrary density. Choose step around the target’s standard deviation; too small and the walk crawls, too large and most moves are rejected. There is no burn-in, so start x0 near the bulk. Cost grows faster than linearly in the number of draws.
simulationsampler-validationreproducibility
def chainMean(m, s, u, p, step, x0) ≔ chainLoop(count(u), m, s, u, p, positions(u), step, x0, 0) / count(u)
chainSpread
chainSpread(m, s, u, p, step, x0) · m target mean scalar · s target standard deviation scalar · u acceptance draws n · p proposal draws n · step half-width of the proposal scalar · x0 starting state scalar → standard deviation scalar
Let \(m\), \(s\), \(u\), \(p\), \(\delta\) and the states \(x^{(0)}, x^{(1)}, \dots, x^{(n)}\) be as in chainMean. Then
The standard deviation of the chain’s states, \(\sqrt{\overline{x^2} - \bar x^2}\): a Monte Carlo estimate of s. It walks the same chain again for the squares and for the mean, so it costs three walks.
Intuition How widely the walker roams. If the chain samples the target correctly, the spread of its visited positions matches the target’s standard deviation; a chain with too small a step, or too short a run, explores a narrower region and underestimates it.
Use it when you check that a chain has explored the whole target, not just its mean: a result well below s signals poor mixing, a start far from the bulk, or too few draws. It divides by n and consecutive states are correlated, so expect Monte Carlo error larger than for independent draws. It walks the chain three times, and calling chainMean as well adds a fourth walk, which matters because long chains are slow.
simulationsampler-validationuncertainty-quantification
def chainSpread(m, s, u, p, step, x0) ≔ sqrt(max(chainSquares(count(u), m, s, u, p, positions(u), step, x0, 0) / count(u) − chainMean(m, s, u, p, step, x0) · chainMean(m, s, u, p, step, x0), 0))
Building blocks
proposed(x, p, ix, k, step)— Let \(x \in \mathbb{R}\), \(p \in [0, 1)^n\), \(\mathit{ix} = (0, \dots, n-1)\), \(k \in \{0, \dots, n-1\}\) and \(\delta > 0\) (the argumentstep). Then \(\operatorname{proposed}(x, p, \mathit{ix}, k, \delta) = x + \delta\,(2p_k - 1)\). The proposal from statexusing the proposal draw at positionk.accepts(m, s, x, y, u, ix, k)— Let \(m \in \mathbb{R}\), \(s > 0\), \(x, y \in \mathbb{R}\), \(u \in [0, 1)^n\), and \(\mathit{ix}\), \(k\) as above. Then \(\operatorname{accepts}(m, s, x, y, u, \mathit{ix}, k) = [\, u_k < \varphi(y; m, s) / \gamma(\varphi(x; m, s)) \,]\). 1 if the move fromxtoyis accepted by the draw at positionk, else 0.stepped(m, s, u, p, ix, k, step, x)— With the same symbols, let \(y = \operatorname{proposed}(x, p, \mathit{ix}, k, \delta)\) and \(a = \operatorname{accepts}(m, s, x, y, u, \mathit{ix}, k)\). Then \(\operatorname{stepped}(m, s, u, p, \mathit{ix}, k, \delta, x) = a\,y + (1 - a)\,x\). The next state: the proposal if accepted, otherwisex.chainLoop(k, m, s, u, p, ix, step, x, run)— Let \(k \in \mathbb{N}\), \(\mathit{run} \in \mathbb{R}\), and the rest as above; write \(x' = \operatorname{stepped}(m, s, u, p, \mathit{ix}, k - 1, \delta, x)\). Then \(\operatorname{chainLoop}(0, \dots, x, \mathit{run}) = \mathit{run}\) and \(\operatorname{chainLoop}(k, \dots, x, \mathit{run}) = \operatorname{chainLoop}(k - 1, \dots, x', \mathit{run} + x')\). Walks the chain forkturns, reading draws at positions \(k-1, k-2, \dots, 0\), accumulating the sum of states.chainSquares(k, m, s, u, p, ix, step, x, run)— With the same symbols and \(x'\), \(\operatorname{chainSquares}(0, \dots, x, \mathit{run}) = \mathit{run}\) and \(\operatorname{chainSquares}(k, \dots, x, \mathit{run}) = \operatorname{chainSquares}(k - 1, \dots, x', \mathit{run} + x'^2)\). Walks the same chain, accumulating the sum of squared states.
Autoregressive fitting
The Yule–Walker estimates use the sample autocorrelations \(r_1, r_2\) of Prelude.Series (autocovariances divided by \(n\)). For AR(1), \(x_t = c + \varphi x_{t-1} + \varepsilon_t\), the estimate is \(\varphi = r_1\). For AR(2), \(x_t = c + \varphi_1 x_{t-1} + \varphi_2 x_{t-2} + \varepsilon_t\):
\[\varphi_1 = \frac{r_1 (1 - r_2)}{1 - r_1^2}, \qquad \varphi_2 = \frac{r_2 - r_1^2}{1 - r_1^2}\]In both cases the constant is \(c = \bar x\,(1 - \sum \varphi)\), so the process mean equals the sample mean. Forecasts are one step ahead. The section ends with simple exponential smoothing, the forecast an ARIMA(0,1,1) model makes.
In the statements below, a series is \(x = (x_1, \dots, x_n) \in \mathbb{R}^n\) with mean \(\bar x = \frac{1}{n}\sum_i x_i\), autocovariances \(C_k = \frac{1}{n}\sum_{i=1}^{n-k} (x_i - \bar x)(x_{i+k} - \bar x)\) and autocorrelations \(r_k = C_k / C_0\).
ar1
ar1(x) · x series n → coefficient scalar
Let \(x \in \mathbb{R}^n\) be a series that is not constant. Then
\[\operatorname{ar1}(x) = r_1 = \frac{\sum_{i=1}^{n-1} (x_i - \bar x)(x_{i+1} - \bar x)}{\sum_{i=1}^{n} (x_i - \bar x)^2} .\]The AR(1) coefficient \(\varphi = r_1\), the lag-one autocorrelation.
Intuition Memory of one step: how much of today’s deviation from the mean survives into tomorrow. Near 1 the series drifts slowly and shocks linger; near 0 each value is almost fresh noise; negative values mean the series tends to flip above and below its mean.
Use it when you model a mean-reverting quantity (temperatures, inventory levels, spreads, sensor drift) or need a quick measure of persistence. A shock decays roughly like \(\varphi^h\) after h steps, so the half-life is about \(\ln 0.5 / \ln \varphi\) for positive \(\varphi\). It assumes a stationary series, so remove trends first. When the lag-two dependence matters, use ar2First and ar2Second.
time-seriespersistence-measurementparameter-estimation
def ar1(x) ≔ autocorrelation(x, 1)
ar1Constant
ar1Constant(x) · x series n → constant scalar
Let \(x \in \mathbb{R}^n\) be a series that is not constant. Then
\[\operatorname{ar1Constant}(x) = \bar x\,(1 - r_1) .\]The AR(1) intercept \(c = \bar x (1 - \varphi)\).
Intuition The constant is not the mean itself but the amount that, together with \(\varphi\) times the previous value, holds the series at its mean. At equilibrium \(x = c + \varphi x\), so \(c = \bar x(1-\varphi)\): the stronger the memory, the smaller the push needed.
Use it when you write the fitted model out as \(x_t = c + \varphi x_{t-1} + \varepsilon_t\), simulate from it, or forecast by hand. It is set so the model’s long-run mean equals the sample mean. For the forecast itself, ar1Forecast combines it for you.
time-seriesparameter-estimation
def ar1Constant(x) ≔ mean(x) · (1 − ar1(x))
ar1Forecast
ar1Forecast(x) · x series n → forecast scalar
Let \(x \in \mathbb{R}^n\) be a series that is not constant. Then
\[\operatorname{ar1Forecast}(x) = \bar x\,(1 - r_1) + r_1\, x_n .\]The one-step-ahead forecast \(c + \varphi\, x_n\) from the last observation.
Intuition Start from the latest value and pull it toward the mean by a fraction \(1 - \varphi\). With strong memory the forecast stays near the last value; with none it is simply the mean.
Use it when you need a quick, explainable next-step prediction for a stationary, mean-reverting series, or a baseline to beat with a richer model. Its typical error variance is ar1Residual. It forecasts one step only; further steps continue to decay toward the mean by \(\varphi\) each step. For a series whose level wanders instead of reverting, smoothForecast is usually better.
forecastingtime-series
def ar1Forecast(x) ≔ ar1Constant(x) + ar1(x) · lastOf(x)
ar1Residual
ar1Residual(x) · x series n → variance scalar
Let \(x \in \mathbb{R}^n\) be a series that is not constant. Then
\[\operatorname{ar1Residual}(x) = \hat\sigma^2\,(1 - r_1^2), \qquad \hat\sigma^2 = \frac{1}{n}\sum_{i=1}^{n} (x_i - \bar x)^2 = C_0 .\]The innovation variance implied by the fit, \(\hat\sigma^2 (1 - \varphi^2)\), where \(\hat\sigma^2\) is the population variance of the series.
Intuition The part of the series’ variance that the previous value cannot predict. A fraction \(\varphi^2\) of the variance is explained by yesterday, and the rest, \(1-\varphi^2\), is new noise arriving every step.
Use it when you need the one-step forecast error variance: its square root is the typical miss of ar1Forecast, and it sets prediction intervals (roughly forecast plus or minus 2 of those) and the noise level for simulating the fitted process. It divides by n, so it is slightly low for short series. A value close to the full variance means the AR(1) model adds little over the mean.
uncertainty-quantificationforecastingtime-series
def ar1Residual(x) ≔ populationVariance(x) · (1 − ar1(x) · ar1(x))
ar2First
ar2First(x) · x series n → coefficient scalar
Let \(x \in \mathbb{R}^n\) be a series that is not constant, with \(r_1^2 \ne 1\). Then
\[\operatorname{ar2First}(x) = \varphi_1 = \frac{r_1\,(1 - r_2)}{1 - r_1^2} .\]The AR(2) coefficient \(\varphi_1\) on lag one.
Intuition The weight on yesterday once the day before is also in the model. It is not the lag-one autocorrelation any more: part of yesterday’s apparent influence was really the day before acting through it, and the Yule–Walker formula separates the two.
Use it when a series shows dynamics one lag cannot capture, such as damped oscillations or cycles (which appear when \(\varphi_2 < 0\) and \(\varphi_1^2 + 4\varphi_2 < 0\)). Use it with ar2Second; the pair only makes sense together. The series should be stationary. If ar2Second is near zero, ar1 is the simpler adequate model.
time-seriesparameter-estimation
def ar2First(x) ≔ autocorrelation(x, 1) · (1 − autocorrelation(x, 2)) / (1 − autocorrelation(x, 1) · autocorrelation(x, 1))
ar2Second
ar2Second(x) · x series n → coefficient scalar
Let \(x \in \mathbb{R}^n\) be a series that is not constant, with \(r_1^2 \ne 1\). Then
\[\operatorname{ar2Second}(x) = \varphi_2 = \frac{r_2 - r_1^2}{1 - r_1^2} .\]The AR(2) coefficient \(\varphi_2\) on lag two, which equals the partial autocorrelation at lag two.
Intuition What the day before yesterday adds after yesterday is already known. That is exactly the lag-two partial autocorrelation: the correlation left at lag two once the chain effect through lag one is removed.
Use it when you decide between AR(1) and AR(2): for a true AR(1) process this value is near zero (roughly within \(\pm 2/\sqrt{n}\)), and a clearly non-zero value says the second lag earns its place. Use it together with ar2First in the model and forecast. The series should be stationary and not constant.
model-selectiontime-seriesparameter-estimation
def ar2Second(x) ≔ pacf2(x)
ar2Constant
ar2Constant(x) · x series n → constant scalar
Let \(x \in \mathbb{R}^n\) be a series, and \(\varphi_1\), \(\varphi_2\) as in ar2First and ar2Second. Then
The AR(2) intercept \(c = \bar x (1 - \varphi_1 - \varphi_2)\).
Intuition The same balancing constant as in AR(1), now with two lags pulling: at equilibrium \(x = c + \varphi_1 x + \varphi_2 x\), so the constant supplies whatever the two coefficients leave out, holding the long-run mean at \(\bar x\).
Use it when you write down, simulate or hand-forecast the fitted AR(2) model. It recomputes both coefficients each time it is called. For the one-step forecast use ar2Forecast.
time-seriesparameter-estimation
def ar2Constant(x) ≔ mean(x) · (1 − ar2First(x) − ar2Second(x))
ar2Forecast
ar2Forecast(x) · x series n, at least 2 → forecast scalar
Let \(x \in \mathbb{R}^n\), \(n \ge 2\), be a series, and \(\varphi_1\), \(\varphi_2\) as in ar2First and ar2Second. Then
The one-step-ahead forecast \(c + \varphi_1 x_n + \varphi_2 x_{n-1}\) from the last two observations.
Intuition Extrapolate from the last two values with the fitted weights, then add the constant that keeps the long-run mean in place. Using two points lets the forecast carry momentum or turn back, which one lag cannot express.
Use it when the series oscillates or has short cycles, or ar2Second is clearly non-zero; otherwise ar1Forecast is simpler and nearly as good. It needs at least two observations and forecasts one step only. Compare the two models on held-out steps rather than on in-sample fit.
forecastingtime-series
def ar2Forecast(x) ≔ ar2Constant(x) + ar2First(x) · lastOf(x) + ar2Second(x) · elementAt(x, count(x) − 2)
smoothed
smoothed(a, x) · a smoothing weight in (0, 1] scalar · x series n → level scalar
Let \(x = (x_1, \dots, x_n) \in \mathbb{R}^n\) be a series and \(0 < a \le 1\). Let \(\ell_0 = x_1\) and
\[\ell_t = a\, x_t + (1 - a)\, \ell_{t-1}, \qquad t = 1, \dots, n .\]Then \(\operatorname{smoothed}(a, x) = \ell_n\). The first update leaves the level at \(x_1\), so the level starts at the first observation and is updated through the series in order.
The final level of simple exponential smoothing. A larger a follows recent observations more closely.
Intuition A running average with fading memory: each new value moves the level a fraction a of the way toward itself. Past values count with weights \(a, a(1-a), a(1-a)^2, \dots\), so recent data matter most and old data fade out gradually.
Use it when you want the current underlying level of a noisy series with no trend or seasonality: demand, latency, sensor readings. Small a (such as 0.1–0.3) smooths heavily and reacts slowly; a = 1 just returns the last value. The level starts at the first observation, so with small a a short series still leans on it. Choose a by one-step forecast error on past data. It does not model trends.
denoisingtime-seriesforecasting
def smoothed(a, x) ≔ smoothLoop(count(x), a, x, positions(x), elementAt(x, 0))
smoothForecast
smoothForecast(a, x) · a smoothing weight in (0, 1] scalar · x series n → forecast scalar
Let \(x \in \mathbb{R}^n\) be a series, \(0 < a \le 1\), and \(\ell_n\) as in smoothed. Then
The one-step-ahead forecast of simple exponential smoothing, which is the final smoothed level.
Intuition If the series is a level that wanders randomly plus noise, the best guess for the next value is the current smoothed level. That is why the forecast equals smoothed, and why the forecast stays flat however far ahead you look.
Use it when you forecast a series whose level drifts without a clear trend or pull back to a fixed mean, the case where AR forecasts revert to a stale mean. It is the forecast an ARIMA(0,1,1) model makes. For mean-reverting data use ar1Forecast; for trending or seasonal data this method will lag behind, since it has no trend term.
forecastingtime-series
def smoothForecast(a, x) ≔ smoothed(a, x)
Building blocks
smoothLoop(k, a, x, ix, level)— Let \(x = (x_0, \dots, x_{n-1}) \in \mathbb{R}^n\) indexed from zero, \(\mathit{ix} = (0, \dots, n-1)\), \(k \in \{0, \dots, n\}\), \(0 < a \le 1\) and \(\ell \in \mathbb{R}\). Then \(\operatorname{smoothLoop}(0, a, x, \mathit{ix}, \ell) = \ell\) and \(\operatorname{smoothLoop}(k, a, x, \mathit{ix}, \ell) = \operatorname{smoothLoop}(k - 1, a, x, \mathit{ix}, a\,x_{n-k} + (1 - a)\,\ell)\). Applies the smoothing update to the lastkobservations, positions \(n-k, \dots, n-1\).
Examples
Inputs are passed as Tensor, because every function here reduces a vector to a number.
Separate a sample into two normal components
Ten values in two tight groups around 1 and 5. EM recovers the two means, an even weight, and each group’s standard deviation \(\sqrt{0.02} = 0.141421\).
open Prelude.Fitting
fn firstMean(x) ↦ mixtureFirstMean(x)
fn secondMean(x) ↦ mixtureSecondMean(x)
fn secondShare(x) ↦ mixtureWeight(x)
fn firstSpread(x) ↦ mixtureFirstSpread(x)
fn secondSpread(x) ↦ mixtureSecondSpread(x)
using MLambda.Hilbert.Runtime;
using Examples;
var data = (Tensor)new[] { 1.0, 1.2, 0.8, 1.1, 0.9, 5.0, 5.2, 4.8, 5.1, 4.9 };
var m1 = Mixture.FirstMean(data);
var m2 = Mixture.SecondMean(data);
var w = Mixture.SecondShare(data);
var s1 = Mixture.FirstSpread(data);
var s2 = Mixture.SecondSpread(data);
Output (verified, six decimals):
m1 = 1
m2 = 5
w = 0.5
s1 = 0.141421
s2 = 0.141421
Estimate a mean and spread with a seeded Metropolis chain
The target is \(\mathcal N(3, 2^2)\). The chain starts at 0, proposes moves of up to ±2, and runs for 2 000 turns on seeded draws. The estimates are close to the true mean 3 and standard deviation 2, within Monte Carlo error; the same seed always gives the same digits. Expect this example to take several minutes (see the limitations above).
open Prelude.Fitting
fn walkMean(m, s, u, p, step, x0) ↦ chainMean(m, s, u, p, step, x0)
fn walkSpread(m, s, u, p, step, x0) ↦ chainSpread(m, s, u, p, step, x0)
using MLambda.Hilbert.Runtime;
using Examples;
var rng = new MLambda.Hilbert.Runtime.Random(7);
var u = rng.Uniform(2000); // acceptance draws
var p = rng.Uniform(2000); // proposal draws
var mean = Chain.WalkMean(3.0, 2.0, u, p, 2.0, 0.0);
var spread = Chain.WalkSpread(3.0, 2.0, u, p, 2.0, 0.0);
Output (verified, six decimals):
mean = 3.04745
spread = 2.072533
Fit an AR(1) model and forecast the next value
The series 1, 3, 1, 3, 1, 3 has mean 2 and lag-one autocorrelation \(-5/6 = -0.833333\). The constant is \(2 \cdot (1 + 5/6) = 3.666667\), so after the last value 3 the forecast is \(3.666667 - 0.833333 \cdot 3 = 1.166667\). The innovation variance is \(1 - 25/36 = 0.305556\).
open Prelude.Fitting
fn lagOne(x) ↦ ar1(x)
fn nextByAr1(x) ↦ ar1Forecast(x)
fn innovation(x) ↦ ar1Residual(x)
fn nextByAr2(x) ↦ ar2Forecast(x)
fn level(a, x) ↦ smoothed(a, x)
fn secondLast(x) ↦ elementAt(x, 1)
using MLambda.Hilbert.Runtime;
using Examples;
var series = (Tensor)new[] { 1d, 3d, 1d, 3d, 1d, 3d };
var phi = Forecast.LagOne(series);
var next = Forecast.NextByAr1(series);
var variance = Forecast.Innovation(series);
Output (verified, six decimals):
ar1 = -0.833333
next = 1.166667
variance = 0.305556
Forecast the same series with AR(2)
For the same series \(r_2 = 2/3\), giving \(\varphi_1 = -10/11\), \(\varphi_2 = -1/11\) and \(c = 4\). The forecast is \(4 - \tfrac{10}{11}\cdot 3 - \tfrac{1}{11}\cdot 1 = 13/11 = 1.181818\). This uses nextByAr2 from the Forecast.hb file above.
var next2 = Forecast.NextByAr2(series);
Output (verified, six decimals):
1.181818
Smooth a short series, and read an element
With weight 0.5 the level of 1, 2, 3 goes 1, then \(0.5 \cdot 2 + 0.5 \cdot 1 = 1.5\), then \(0.5 \cdot 3 + 0.5 \cdot 1.5 = 2.25\). elementAt counts from 0, so position 1 of 10, 20, 30 is 20. Both use the Forecast.hb file above.
var level = Forecast.Level(0.5, (Tensor)new[] { 1d, 2d, 3d });
var second = Forecast.SecondLast((Tensor)new[] { 10d, 20d, 30d });
Output (verified):
smoothed = 2.25
elementAt = 20