MonotonePSpline

Shape-constrained P-spline: monotone increasing or decreasing.

Usage

Source

MonotonePSpline(
    k=20,
    degree=3,
    m=2,
    decreasing=False,
)

Uses the same B-spline basis and difference penalty as ~whittaker.smooths.pspline.PSpline, but additionally requires the fitted curve f(x) = \sum_j \beta_j B_j(x) to be non-decreasing (or, with decreasing=True, non-increasing) over the whole domain. This is enforced as a linear inequality constraint on the coefficients rather than on f directly: because each B-spline basis function B_j is non-negative and has local support (a “bump” that overlaps only its neighbours), a non-decreasing sequence of coefficients \beta_1 \le \beta_2 \le \dots \le \beta_k guarantees a non-decreasing curve. Intuitively, moving x to the right shifts the basis functions’ weight away from earlier, smaller-or-equal coefficients and onto later, larger-or-equal ones, so the weighted sum can only stay flat or increase. Use MonotonePSpline (via s(x, bs="mpi") for increasing or s(x, bs="mpd") for decreasing in a formula) whenever domain knowledge says the relationship must be monotone — e.g. a dose-response curve, a cumulative distribution, or a growth curve — and an unconstrained smooth would otherwise wiggle non-monotonically due to noise.

Parameters

k: int = 20

Number of B-spline basis functions. The default is 20.

degree: int = 3

B-spline polynomial degree. The default is 3 (cubic).

m: int = 2

Difference penalty order. The default is 2.

decreasing: bool = False
If True, enforce monotone decreasing. The default is False (monotone increasing).

Notes

The monotonicity constraint is enforced during fitting, not by direct constrained optimization. Instead, at each penalized iteratively reweighted least squares (P-IRLS) iteration the ordinary (unconstrained) coefficient update is projected onto the monotone cone: whittaker.fitting.pirls detects any smooth term whose basis is a MonotonePSpline and passes its coefficient block through project_monotone before the next iteration’s linear predictor is formed. This projection uses the Pool Adjacent Violators Algorithm (PAVA), the standard algorithm for isotonic regression (Barlow, Bartholomew, Bremner & Brunk, 1972; Best & Chakravarti, 1990). Given an arbitrary vector, PAVA finds the closest non-decreasing vector to it in the least-squares sense by scanning for adjacent “violations” (a value followed by a smaller one) and replacing each violating block with its mean, merging blocks until no violations remain. Because this is an orthogonal projection onto the convex cone of non-decreasing sequences, iterating it alongside the P-IRLS coefficient update drives the fit toward the coefficient vector, within that cone, that best balances the penalized deviance and the constraint.

Examples

import numpy as np
import whittaker as wt

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 1, 200))
y = 3 * x + rng.normal(scale=0.15, size=200)

model = wt.GAM("y ~ s(x, bs='mpi')").fit({"x": x, "y": y})

new_x = np.linspace(0, 1, 50)
fitted = model.predict({"x": new_x}).values
np.all(np.diff(fitted) >= -1e-8)
np.True_

Attributes

Name Description
constraint_direction Sign convention used by the PAVA projection.

constraint_direction

Sign convention used by the PAVA projection.

constraint_direction: int

Returns -1 when decreasing=True (the coefficients are projected onto the non-increasing cone by negating, applying PAVA, and negating back) or 1 for the default non-decreasing case. Consulted by whittaker.fitting.pirls when applying project_monotone to this term’s coefficient block.

Methods

Name Description
null_space_dimension() Return the dimension of the basis’s unpenalized null space.

null_space_dimension()

Return the dimension of the basis’s unpenalized null space.

Usage

Source

null_space_dimension()

The monotonicity constraint is enforced by post-hoc projection rather than by removing degrees of freedom from the penalty, so the underlying P-spline penalty null space is treated as fully absorbed elsewhere in the model; this basis reports 0 so it does not additionally compete with a model intercept for identifiability.

Returns

int
Always 0.