Shape-constrained P-spline: monotone increasing or decreasing.
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)
Attributes
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
null_space_dimension()
Return the dimension of the basis’s unpenalized null space.
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.