StreamingGAM

Streaming / online GAM.

Usage

Source

StreamingGAM(
    formula,
    *,
    family=None,
    decay=1.0,
    smoothing_params=None,
)

Fits a GAM incrementally by accumulating weighted sufficient statistics (X'WX and X'Wz) across data batches, rather than storing and re-fitting on the full dataset each time. This makes it suitable for data arriving continuously or in a stream too large to hold in memory at once: the memory footprint depends only on the number of basis coefficients p (via a p x p matrix), not on the number of observations seen. The model structure (formula, basis dimensions, penalties) is fixed at initialisation from a small pilot batch, and subsequent partial_fit() calls add data without storing the raw observations.

Two update modes are supported: accumulate mode (decay=1.0, the default), where every batch contributes equally regardless of when it arrived, and sliding-window mode (decay < 1.0), where older batches’ contributions to the accumulated statistics are exponentially downweighted, allowing the model to adapt to a slowly drifting data-generating process.

Use StreamingGAM for large or continuously-arriving datasets where holding the full data in memory (as ordinary GAM does) is impractical, or where the underlying relationship may drift over time and old data should be forgotten.

Parameters

formula: str | Formula

Model formula (e.g. "y ~ s(x1) + s(x2)").

family: Family | None = None

Response distribution family. Defaults to Gaussian().

decay: float = 1.0

Exponential decay factor for sliding window. 1.0 (default) means no decay (accumulate all data equally). Values less than 1.0 downweight older batches’ sufficient statistics by a factor of decay every time a new batch arrives.

smoothing_params: list[float] | None = None
Fixed smoothing parameters. If None, estimated from the pilot batch (via a one-off ordinary GAM fit with REML) and optionally re-estimated later via solve(reestimate_smoothing=True).

Notes

Each partial_fit() call treats the incoming batch as one step of iteratively reweighted least squares: given the current coefficients, it computes working responses z and IRLS weights W for the batch, forms the batch’s contribution to the weighted normal equations,

X_{\text{batch}}' W X_{\text{batch}}, \qquad X_{\text{batch}}' W z_{\text{batch}},

and adds these to the running totals (after applying the decay factor, if any, to the existing totals). Calling solve() then solves the penalized normal equations

\left(\sum_{\text{batches}} X'WX + S_\lambda\right) \beta = \sum_{\text{batches}} X'Wz

via a Cholesky factorization, where S_\lambda = \sum_j \lambda_j S_j is the weighted sum of penalty matrices. Because only the accumulated p x p matrix X'WX and length-p vector X'Wz are retained, this scales to arbitrarily many observations at fixed memory cost in p.

Examples

import numpy as np
from whittaker.streaming import StreamingGAM

rng = np.random.default_rng(0)
model = StreamingGAM("y ~ s(x)")

for _ in range(5):
    x = rng.uniform(0, 1, 200)
    y = np.sin(2 * np.pi * x) + rng.normal(scale=0.2, size=200)
    model.partial_fit({"x": x, "y": y})

model.solve()
print(model.summary())
StreamingGAM summary
============================================================
Formula:     y ~ s(x)
Family:      Gaussian
Decay:       1.0
N obs:       1000
N batches:   5
EDF total:   9.6
Scale:       0.3044
Deviance:    301.5

Smooth terms:
  s(x): edf=8.6

Snapshots:   1

Attributes

Name Description
coefficients Current coefficient vector from the most recent solve() (or the pilot fit).
edf_total Total effective degrees of freedom from the most recent solve().
family Response distribution family used for the link function and variance model.
formula Parsed model formula this streaming GAM was constructed with.
is_initialised Whether the model structure has been built from a pilot batch.
is_solved Whether coefficients are available from an initial fit or a solve() call.
n_batches Total number of batches passed to partial_fit() so far.
n_obs Accumulated observation count across all ingested batches.
scale Estimated dispersion (scale) parameter from the most recent solve().
smoothing_params Smoothing parameters currently in use, one per penalized smooth term.

coefficients

Current coefficient vector from the most recent solve() (or the pilot fit).

coefficients: NDArray

Returns a copy, so mutating the returned array does not affect the model.


edf_total

Total effective degrees of freedom from the most recent solve().

edf_total: float

Computed as the trace of the hat matrix implied by the accumulated X'WX and the current smoothing parameters.


family

Response distribution family used for the link function and variance model.

family: Family


formula

Parsed model formula this streaming GAM was constructed with.

formula: Formula

The formula (response, smooth and parametric terms) is fixed for the lifetime of the model; it determines the basis structure built from the first partial_fit() batch.


is_initialised

Whether the model structure has been built from a pilot batch.

is_initialised: bool

Becomes True after the first partial_fit() call, once the model matrix, basis dimensions, and penalties have been established.


is_solved

Whether coefficients are available from an initial fit or a solve() call.

is_solved: bool

predict(), coefficients, and other solved-state accessors raise RuntimeError when this is False.


n_batches

Total number of batches passed to partial_fit() so far.

n_batches: int

Unlike n_obs, this count is never decayed.


n_obs

Accumulated observation count across all ingested batches.

n_obs: int

Under sliding-window mode (decay < 1.0), this count is itself decayed at each partial_fit() call, so it reflects an effective rather than a raw cumulative count.


scale

Estimated dispersion (scale) parameter from the most recent solve().

scale: float

For families with a known scale (e.g. Poisson, binomial), this is fixed at 1.0; otherwise it is estimated from the accumulated deviance and edf_total.


smoothing_params

Smoothing parameters currently in use, one per penalized smooth term.

smoothing_params: list[float]

These come from the constructor’s smoothing_params argument if fixed, or from the pilot fit / most recent solve(reestimate_smoothing=True) call otherwise. Returns a copy, so mutating the returned list does not affect the model.

Methods

Name Description
partial_fit() Ingest a batch of data.
predict() Predict on new data.
reset() Reset accumulated statistics (keep model structure).
should_refit() Check if a refit (solve) is recommended.
smoothing_history() Return history of snapshots taken at each solve() call.
solve() Solve for coefficients from accumulated sufficient statistics.
summary() Build a human-readable text summary of the streaming GAM’s current state.

partial_fit()

Ingest a batch of data.

Usage

Source

partial_fit(data)

On the first call, builds the model matrix structure (basis dimensions, penalties, knot locations) from this batch via a pilot GAM fit, and initializes the coefficients and smoothing parameters from that fit (unless fixed smoothing parameters were supplied). Subsequent calls reuse this fixed structure: the batch’s working response and IRLS weights are computed using the current coefficients (from the most recent solve(), or the pilot fit), and its contribution to the running X'WX / X'Wz sufficient statistics is added (after applying the decay factor to existing totals, if decay < 1.0).

Parameters

data: InputData
Column-oriented batch data containing the response and all covariates in formula.

Returns

StreamingGAM
Returns self for method chaining.

predict()

Predict on new data.

Usage

Source

predict(
    new_data,
    *,
    se=False,
)

Builds the prediction design matrix using the fixed basis structure established at initialisation, and forms the linear predictor and (via the family’s inverse link) fitted values from the current coefficients. Requires that solve() has been called at least once.

Parameters

new_data: InputData

Column-oriented covariate data.

se: bool = False
If True, compute standard errors on the linear predictor scale from the Bayesian posterior covariance implied by the current accumulated X'WX and smoothing parameters.

Returns

PredictionResult

reset()

Reset accumulated statistics (keep model structure).

Usage

Source

reset()

Returns

StreamingGAM

should_refit()

Check if a refit (solve) is recommended.

Usage

Source

should_refit(
    *,
    min_batches=10,
)

Returns True when enough new data has accumulated since the last solve to justify re-estimating coefficients.

Parameters

min_batches: int = 10
Minimum batches since last solve before recommending refit.

smoothing_history()

Return history of snapshots taken at each solve() call.

Usage

Source

smoothing_history()

Returns

list[StreamingSnapshot]

solve()

Solve for coefficients from accumulated sufficient statistics.

Usage

Source

solve(
    *,
    reestimate_smoothing=False,
)

Forms the penalized normal equations (X'WX + S_lambda) beta = X'Wz from the currently accumulated statistics and the current smoothing parameters, and solves them via a Cholesky factorization (with a small ridge added for numerical stability). Also updates effective degrees of freedom, the scale estimate, and appends a StreamingSnapshot to the fit history.

Parameters

reestimate_smoothing: bool = False
If True, re-estimate smoothing parameters from the current sufficient statistics using a GCV line search over each smoothing parameter in turn (coordinate-wise). Default False uses the current (pilot or fixed) smoothing parameters. Ignored (has no effect) if the model was constructed with fixed smoothing_params.

Returns

StreamingGAM
Returns self for method chaining.

summary()

Build a human-readable text summary of the streaming GAM’s current state.

Usage

Source

summary()

Reports the formula, family, decay factor, accumulated observation and batch counts, total EDF, scale, and accumulated deviance, followed by a per-smooth-term EDF breakdown and the number of solve() snapshots recorded so far.

Returns

str
Multi-line summary text.