Smoothing

fdatools.smoothing: smooth() with GCV or degrees-of-freedom selection of λ, positive / monotone / morph constraints, irregular designs, and the scikit-learn Smoother estimator.

Penalised smoothing: turn discrete observations into functional data.

A single entry point, :func:smooth, replaces the twelve R fda smoothing functions (smooth.basis, smooth.basisPar, Data2fd, smooth.monotone, smooth.pos, smooth.morph, lambda2df, lambda2gcv, df2lambda and friends). The unconstrained problem is the penalised least-squares fit of Ramsay & Silverman chapter 5,

.. math::

(\Phi^{T} W \Phi + \lambda R)\,c = \Phi^{T} W y ,
\qquad R_{ij} = \int (L\varphi_i)(L\varphi_j)\,dt ,

whose hat matrix H = Φ (ΦᵀWΦ + λR)⁻¹ ΦᵀW gives the equivalent degrees of freedom df = tr H and the generalised cross-validation criterion gcv_i = (SSE_i / n) / (1 - df/n)² for curve i.

Constrained fits reuse the same penalty on a latent function W: constraint="positive" fits x = exp W, constraint="monotone" fits x(t) = β₀ + β₁ ∫ₐᵗ exp W(u) du, and constraint="morph" is the monotone fit with β pinned so that the fitted curve maps the domain onto itself. Their derivatives of any order are exact: with w_j = D^j W, Faà di Bruno's formula gives Dⁿ exp W = exp(W) Bₙ(w_1, …, w_n) for the complete Bell polynomial Bₙ (B₀ = 1, B_{n+1} = Σ_{i=0}^{n} C(n, i) B_{n-i} w_{i+1}), and Dⁿ x = β₁ D^{n-1} exp W for a monotone fit.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.smoothing import smooth
>>> t = np.linspace(0.0, 1.0, 21)
>>> y = np.sin(2 * np.pi * t)
>>> result = smooth(y, t, basis=BSpline(domain=(0.0, 1.0), n_basis=12), lam=1e-8)
>>> bool(np.max(np.abs(result.fd(t)[:, 0] - y)) < 1e-3)
True

SmoothResult dataclass

SmoothResult(
    fd: FData,
    df: float,
    gcv: Array,
    sse: float,
    penalty_matrix: Array,
    lam: float,
    y2c_map: Array,
    beta: Array | None = None,
    constraint: str | None = None,
)

Outcome of a :func:smooth call.

Attributes:

Name	Type	Description
fd	FData	The fitted curves. For a constrained fit this holds the latent function W; call the result itself (result(t, deriv)) to evaluate the constrained curve exp W or β₀ + β₁ ∫ exp W and its derivatives of any order.
df	float	Equivalent degrees of freedom, tr H. For a constrained fit this is the trace of the hat matrix of the final Gauss-Newton linearisation, and for per-curve irregular designs it is the sum over the independent fits.
gcv	array	Generalised cross-validation score per curve, shaped like the curve axes of y (a scalar for a single unnamed curve, (n_curves,) for a matrix, (n_curves, n_vars) for a multivariate one).
sse	float	Residual sum of squares of the fit, summed over every point and curve. Unweighted even when weights is given, matching R's smooth.basis.
penalty_matrix	array	The roughness penalty R.
lam	float	The smoothing parameter actually used -- the value selected by GCV or by the degrees-of-freedom target when one of those was requested.
y2c_map	array or tuple of arrays	The linear map (ΦᵀWΦ + λR)⁻¹ ΦᵀW from data to coefficients, of shape (n_basis, n_t). A tuple with one entry per curve for an irregular design, and None for a constrained fit, which is not linear in y.
beta	array or None	(2, n_curves) location/scale pair of a monotone or morph fit; None otherwise.
constraint	str or None	The constraint that was applied, if any.

Examples:

>>> import numpy as np
>>> from fdatools.smoothing import smooth
>>> t = np.linspace(0.0, 1.0, 15)
>>> result = smooth(t**2, t, lam=1e-8)
>>> result.constraint is None and result.beta is None
True

__call__

__call__(t: Any, deriv: int = 0) -> Array

Evaluate the fitted curves or one of their derivatives.

For a constrained fit the derivatives are exact to rounding error, for any order: Dⁿ exp W = exp(W) Bₙ(DW, …, DⁿW) with the complete Bell polynomial Bₙ (Faà di Bruno's formula), and a monotone curve has Dⁿ x = β₁ D^{n-1} exp W for n >= 1. The value itself (deriv=0) of a monotone curve is β₀ + β₁ ∫ₐᵗ exp W, integrated by Gauss-Legendre quadrature between the basis break points.

Parameters:

Name	Type	Description	Default
t	array	Points at which to evaluate.	required
deriv	int	Derivative order, 0 (the default) or more.	0

Returns:

Type	Description
array	Fitted values with the shape of fd(t).

Raises:

Type	Description
ValueError	If deriv is negative, or a monotone fit carries no beta.

Examples:

>>> import numpy as np
>>> from fdatools.smoothing import smooth
>>> t = np.linspace(0.0, 1.0, 40)
>>> fit = smooth(np.exp(2 * t), t, constraint="positive", lam=1e-8)
>>> np.round(fit(np.array([0.5]), 3)[:, 0] / np.exp(1.0), 2)  # D³ e^{2t} = 8 e^{2t}
array([8.])

Smoother

Smoother(
    basis: Basis | None = None,
    *,
    t: Any = None,
    lam: float | str = "gcv",
    penalty: int | LDO = 2,
    weights: Any = None,
)

Bases: TransformerMixin, BaseEstimator

Penalised smoothing as a scikit-learn transformer.

Each row of X is one curve sampled on the shared grid t; transform returns the basis coefficients, so a smoother can head a pipeline whose later steps work on the coefficient representation.

Parameters:

Name	Type	Description	Default
basis	Basis	Expansion basis. None builds the default cubic B-spline of :func:smooth from the fitting grid.	None
t	array	Sampling grid of length n_features. When omitted, fit uses 0, 1, ..., n_features - 1; it may also be supplied per call as the t fit parameter.	None
lam	float or str	Smoothing parameter, as in :func:smooth.	'gcv'
penalty	int or LDO	Roughness operator, as in :func:smooth.	2
weights	array	Observation weights of length n_features.	None

Attributes:

Name	Type	Description
fd_	FData	The fitted curves, one per row of the training X.
result_	SmoothResult	The full fit, including df, gcv and the data-to-coefficient map.
t_	ndarray	The grid actually used.
n_features_in_	int	Number of sampling points seen during fit.

Examples:

>>> import numpy as np
>>> from fdatools.smoothing import Smoother
>>> rng = np.random.default_rng(0)
>>> t = np.linspace(0.0, 1.0, 25)
>>> X = np.sin(2 * np.pi * t) + 0.05 * rng.standard_normal((6, 25))
>>> coefs = Smoother(t=t, lam=1e-4).fit_transform(X)
>>> coefs.shape[0]
6

fit

fit(X: Any, y: Any = None, *, t: Any = None) -> Smoother

Smooth every row of X and store the resulting curves.

Parameters:

Name	Type	Description	Default
X	array of shape (n_samples, n_features)	One sampled curve per row.	required
y	ignored	Present for API compatibility.	None
t	array	Sampling grid for this call, overriding the one given at construction.	None

Returns:

Type	Description
Smoother	The fitted estimator.

transform

transform(X: Any) -> Array

Return the basis coefficients of every row of X.

Parameters:

Name	Type	Description	Default
X	array of shape (n_samples, n_features)	Curves sampled on the grid seen during fit.	required

Returns:

Type	Description
ndarray	Coefficients of shape (n_samples, n_basis).

get_feature_names_out

get_feature_names_out(input_features: Any = None) -> Array

Return the basis-function names produced by :meth:transform.

Parameters:

Name	Type	Description	Default
input_features	ignored	Present for API compatibility; the output names come from the basis.	None

Returns:

Type	Description
ndarray	One name per basis function.

lambda_to_df

lambda_to_df(
    t: Any,
    basis: Basis,
    lam: float,
    *,
    penalty: int | LDO = 2,
    weights: Any = None,
) -> float

Equivalent degrees of freedom of a smooth at a given lam.

Parameters:

Name	Type	Description	Default
t	array	Observation points.	required
basis	Basis	Expansion basis.	required
lam	float	Smoothing parameter, on the linear scale.	required
penalty	int or LDO	Roughness operator; default the curvature penalty D².	2
weights	array	Observation weights.	None

Returns:

Type	Description
float	tr H, the trace of the hat matrix.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.smoothing import lambda_to_df
>>> t = np.linspace(0.0, 1.0, 30)
>>> round(lambda_to_df(t, BSpline(domain=(0.0, 1.0), n_basis=10), 0.0), 6)
10.0

df_to_lambda

df_to_lambda(
    t: Any,
    basis: Basis,
    df: float,
    *,
    penalty: int | LDO = 2,
    weights: Any = None,
) -> float

Smoothing parameter whose smooth has df degrees of freedom.

Parameters:

Name	Type	Description	Default
t	array	Observation points.	required
basis	Basis	Expansion basis.	required
df	float	Target equivalent degrees of freedom.	required
penalty	int or LDO	Roughness operator; default the curvature penalty D².	2
weights	array	Observation weights.	None

Returns:

Type	Description
float	The lam solving lambda_to_df(t, basis, lam) == df.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.smoothing import df_to_lambda, lambda_to_df
>>> t = np.linspace(0.0, 1.0, 40)
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=12)
>>> lam = df_to_lambda(t, basis, 6.0)
>>> round(lambda_to_df(t, basis, lam), 8)
6.0

gcv_curve

gcv_curve(
    y: Any,
    t: Any,
    basis: Basis,
    lambdas: Any,
    *,
    penalty: int | LDO = 2,
    weights: Any = None,
) -> Array

GCV score of every curve over a grid of smoothing parameters.

One eigendecomposition serves the whole grid, so the cost is O(n_lambda · n_t · n_basis) rather than a factorisation per lambda.

Parameters:

Name	Type	Description	Default
y	array	Observations of shape (n_t,), (n_t, n_curves) or (n_t, n_curves, n_vars).	required
t	array	Observation points, shape (n_t,).	required
basis	Basis	Expansion basis.	required
lambdas	array	Smoothing parameters on the linear scale.	required
penalty	int or LDO	Roughness operator; default the curvature penalty D².	2
weights	array	Observation weights.	None

Returns:

Type	Description
array	Shape (n_lambda, *curve_shape): one GCV score per lambda per curve.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.smoothing import gcv_curve
>>> t = np.linspace(0.0, 1.0, 40)
>>> y = np.sin(2 * np.pi * t)
>>> scores = gcv_curve(y, t, BSpline(domain=(0.0, 1.0), n_basis=10), [1e-6, 1e-2])
>>> scores.shape
(2,)

smooth

smooth(
    y: Any,
    t: Any,
    *,
    basis: Basis | None = None,
    lam: float | str = "gcv",
    penalty: int | LDO = 2,
    constraint: str | None = None,
    weights: Any = None,
    df: float | None = None,
) -> SmoothResult

Fit smooth curves to discrete observations by penalised least squares.

Parameters:

Name	Type	Description	Default
y	array or sequence of arrays	Observations of shape (n_t,), (n_t, n_curves) or (n_t, n_curves, n_vars). For a per-curve irregular design, a sequence of 1-D arrays, one per curve.	required
t	array or sequence of arrays	Observation points shared by every curve, or one vector per curve for an irregular design (in which case y must be a matching sequence).	required
basis	Basis	Expansion basis. Defaults to a cubic B-spline on [min t, max t] with min(n_t, 40) + 2 functions -- deterministic, and independent of the observed values.	None
lam	float or str	Smoothing parameter. A number is used directly; "gcv" (the default) minimises the mean GCV score over log10 λ ∈ [-8, 8]; "df=12" solves for the λ with 12 degrees of freedom.	'gcv'
penalty	int or LDO	Roughness operator L. An integer means D^penalty; the default 2 is the curvature penalty.	2
constraint	(None, 'positive', 'monotone', 'morph')	Shape constraint. "positive" fits x = exp W, "monotone" fits x = β₀ + β₁ ∫ₐᵗ exp W, and "morph" is the monotone fit with β pinned so that x maps the domain onto itself.	None
weights	array	Observation weights of shape (n_t,).	None
df	float	Target degrees of freedom; equivalent to lam="df=<value>" and overriding lam.	None

Returns:

Type	Description
SmoothResult	The fit, its degrees of freedom, GCV scores and data-to-coefficient map.

Raises:

Type	Description
ValueError	If the shapes of y and t disagree, if constraint is unknown, or if a constrained fit is asked for on an irregular design.

Examples:

>>> import numpy as np
>>> from fdatools.smoothing import smooth
>>> t = np.linspace(0.0, 1.0, 60)
>>> y = np.exp(np.sin(4 * t))
>>> fit = smooth(y, t, constraint="positive", lam=1e-6)
>>> bool(np.min(fit(t)) > 0.0)
True