Feature: Hilbert Space Gaussian Process (HSGP) term for smooth nonparametric effects

[daimon-pymclabs]

Summary

pathmc currently has no Gaussian Process support of any kind — a grep for hsgp, gaussian process, pm.gp, Hilbert returns nothing across code, docs, and tests. This issue scopes adding a Hilbert Space Gaussian Process (HSGP) term to the DSL, giving pathmc a nonparametric smooth analogous to Bambi's hsgp() formula term (bambi HSGP 1D notebook).

Prior art worth reviewing for design cues:

• Bambi's closed HSGP PRs: https://github.com/bambinos/bambi/pulls?q=is%3Apr+is%3Aclosed+hsgp
• Bambi issue framing the feature: Bambi + HSGP + nutpie crashes during arviz conversion. bambinos/bambi#982

Related pathmc issue: #98 (time-varying coefficients for panel models) — an HSGP over time is one way to express smoothly-varying effects, so these should stay design-aware of each other.

Motivation

pathmc models are strictly linear-in-coefficients along each path. Any nonlinear-but-smooth dependence (a smooth trend over time, a smooth dose–response that isn't captured by logistic_saturation, a smooth confounder adjustment) currently has no expression. HSGP is the standard scalable way to add such a smooth: it approximates a GP with a finite Laplacian eigenfunction basis, so it reduces to a linear model in a fixed basis matrix — a natural fit for a library already built around design matrices and pm.do().

Why it isn't a drop-in today

The closest existing machinery is the transform system (pathmc/transforms.py): a registry (adstock, logistic_saturation) of parametrized, estimable, composable functions with automatic priors, panel/scan-awareness, and pm.do() compatibility. That's the right conceptual neighbor, but not a drop-in:

• Transforms map one input column → one output tensor.
• HSGP must expand one input column → an M-column basis matrix that multiplies M free coefficients.

So the bottleneck is the design-matrix layer, not the kernel math. build_design_matrix in compile.py emits one column per term, and the two resolvers (_make_cross_sectional_resolver, _make_scan_resolver) dispatch on slot.kind. HSGP needs a new slot kind plus basis expansion at that layer.

Good news: PyMC already ships pm.gp.HSGP, so the eigenfunction basis and its scaling are free — the work is wiring, not deriving.

Scope

Phase 1 (this issue): 1-D HSGP term, cross-sectional path only. Parity with Bambi's simplest hsgp_1d notebook. Estimated 3–5 focused days.

Follow-ups (separate issues):

• Panel / scan resolver path — HSGP interacts awkwardly with pytensor.scan (stateful, lagged). Coordinate with Feature: Time-varying coefficients for panel models #98.
• Multi-dimensional inputs (hsgp(x1, x2)).
• Grouped GPs (by=, Bambi-style).

Estimate widens to 1.5–2 weeks if panel-mode, multi-D, and grouped GPs are pulled into the first pass.

File-by-file breakdown

Layer	File	Work
Parse	pathmc/parse.py (Term, _parse_term)	Recognize hsgp(var, m=, c=, ...) syntax; produce a dedicated HSGP node (do not overload TransformCall)
Compile	pathmc/compile.py (PredictorSlot, build_design_matrix, build_mu_specs, both resolvers)	New kind="hsgp"; materialize the Laplacian eigenfunction basis; route basis × coeff in resolvers
Priors	pathmc/priors.py (default_priors, _collect_transform_defaults)	Hyperpriors for amplitude (eta) + lengthscale (ell); coefficient priors for the M basis weights
Introspect	pathmc/introspect.py (build_equations, build_dag_viz)	Render HSGP terms in equations / DAG
Tests	new tests/test_hsgp.py + additions to test_compile.py, test_priors.py	Basis shape, prior wiring, pm.do() recomputation
Docs	new example under docs/examples/	Mirror Bambi's hsgp_1d

Suggested implementation approach

1. DSL syntax. Add an hsgp(...) term recognized in _parse_term. Suggested signature mirroring Bambi/PyMC: hsgp(x, m=<num basis>, c=<boundary factor>, cov="expquad", by=None). Emit a new HSGPTerm/HSGPCall AST node rather than reusing TransformCall, since the downstream contract (basis matrix, not scalar) differs.

2. Design-matrix expansion. In build_design_matrix, when a slot is kind="hsgp", do not emit a single patsy column. Instead record the raw input column(s) and the HSGP config, and defer basis construction to compile time. Concretely, build the basis with phi, sqrt_psd = pm.gp.HSGP(m=[m], c=c, cov_func=...).prior_linearized(X) (or the current PyMC API equivalent) so the term becomes phi @ (sqrt_psd * beta_hsgp). Keep X as pm.Data so interventions flow through.

3. New slot kind. Add kind="hsgp" to PredictorSlot, carrying (input_name, m, c, cov, by). Both resolvers dispatch on this: the cross-sectional resolver returns phi @ (sqrt_psd * beta); the scan resolver raises NotImplementedError for now with a clear message ("HSGP not yet supported in panel/scan models — see follow-up").

4. Priors. Register defaults in priors.py: beta_hsgp ~ Normal(0, 1) (the standardized basis weights), ell ~ InverseGamma/LogNormal (lengthscale), eta ~ HalfNormal (amplitude). Wire them through the same customization path transforms use so users can override via the priors API.

5. pm.do() correctness. Because the basis is a deterministic function of X (a pm.Data), an intervention do(x=...) must recompute phi from the new X. Verify the basis is built inside the model graph from the pm.Data, not precomputed as a numpy constant — add an explicit test that phi changes under intervention.

6. Guardrails. HSGP must be a terminal term — reject nesting like adstock(hsgp(...)) and hsgp(...) inside interactions in Phase 1 with a clear parse-time error. Also guard against panel/scan usage until the follow-up lands.

7. Introspection. Render as e.g. f_{\text{hsgp}}(x) in build_equations, and mark the node in the DAG viz so users see the smooth is nonparametric.

Acceptance criteria (Phase 1)

• y ~ hsgp(x, m=..., c=...) parses, compiles, and samples on a cross-sectional model.
• Recovers a known smooth function on synthetic data (test), comparable to Bambi's hsgp_1d.
• pm.do() intervention on x recomputes the basis (test).
• Clear errors for panel/scan use, nesting, and interaction use.
• Docs example mirroring Bambi's hsgp_1d.