Regression

fdatools.regression: fregress() for scalar and functional responses with scalar and functional covariates (model type read from the arguments, or an R-style formula), with predict(), stderr() and cv() on the result, and the scikit-learn FRegress estimator.

Functional linear regression: :func:fregress and the :class:FRegress estimator.

One entry point replaces R fda's fRegress family (fRegress, fRegress.formula, fRegress.fd, fRegress.double, predict.fRegress, fRegress.stderr and fRegress.CV). The kind of model is read off the arguments:

• scalar response (y a vector of n numbers). Each covariate enters as a scalar z_ij β_j or, when it is a curve, as ∫ x_ij(t) β_j(t) dt:

.. math:: y_i = \sum_j \int x_{ij}(t)\,\beta_j(t)\,dt + \varepsilon_i .

• functional response (y an :class:~fdatools.core.FData of n curves). Scalar covariates give the function-on-scalar model and curves give the concurrent model; the two mix freely:

.. math:: y_i(t) = \sum_j x_{ij}(t)\,\beta_j(t) + \varepsilon_i(t) .

Writing every coefficient function in its own basis, β_j = θ_jᵀ b_j, both models are penalised least-squares problems (Ramsay, Hooker & Graves 2009, chapters 9-10) whose normal equations C b = D stack one block per term:

.. math::

C_{jk} = \sum_i w_i \int x_{ij} x_{ik}\,\theta_j \theta_k^{T}\,dt
         + \delta_{jk}\,\lambda_j R_j ,
\qquad
D_j = \sum_i w_i \int x_{ij}\,y_i\,\theta_j\,dt

for a functional response, and C = ZᵀWZ + diag(λ_j R_j), D = ZᵀWy with Z_i = (∫ x_ij θ_j)_j for a scalar one. R_j = ∫ (Lθ_j)(Lθ_j)ᵀ is the roughness penalty of the j-th coefficient. Every integral is computed exactly (Gauss-Legendre on the break points of the bases involved), whereas R approximates them -- see tests/parity/test_regression.py for the measured consequences.

:func:linmod fits the fully functional model with a bivariate coefficient (R's linmod), y_i(t) = alpha(t) + ∫ x_i(s) β(s, t) ds + ε_i(t), where a response value depends on the whole covariate curve.

Both entry points compute in the array namespace of their inputs: PyTorch coefficients, covariates or weights give PyTorch results, and gradients flow back to the inputs.

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> rng = np.random.default_rng(0)
>>> z = rng.standard_normal(40)
>>> y = 1.0 + 2.0 * z + 0.01 * rng.standard_normal(40)
>>> model = fregress(y, [1.0, z])
>>> np.round([float(b.coefs[0, 0]) for b in model.beta], 2).tolist()
[1.0, 2.0]

FRegressStderr dataclass

FRegressStderr(
    beta: tuple[FData, ...], fitted: Any, cov: Array
)

Sampling variability of a fitted functional regression.

Replaces the output of R's fRegress.stderr.

Attributes:

Name	Type	Description
beta	tuple of FData	Pointwise standard error sqrt(θ_j(t)ᵀ V_jj θ_j(t)) of each coefficient function, projected (in L2) onto the coefficient's basis.
fitted	array or FData	Standard error of the fitted values: a vector for a scalar response, curves in the response basis for a functional one.
cov	array	Covariance V of the stacked coefficient vector b (R's bvar).

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> z = np.linspace(-1.0, 1.0, 20)
>>> y = 3.0 * z + np.sin(7.0 * z)
>>> se = fregress(y, [1.0, z]).stderr()
>>> se.cov.shape
(2, 2)

FRegressCV dataclass

FRegressCV(sse: float, errors: Any)

Leave-one-out cross-validation of a functional regression.

Replaces the output of R's fRegress.CV.

Attributes:

Name	Type	Description
sse	float	Sum over observations of the squared leave-one-out error (integrated over the domain for a functional response).
errors	array or FData	y_i - ŷ_i^{(-i)}, the error of the model fitted without observation i: a vector for a scalar response, curves in the response basis for a functional one.

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> z = np.linspace(-1.0, 1.0, 20)
>>> cv = fregress(2.0 * z, [1.0, z]).cv()
>>> bool(cv.sse < 1e-20)
True

FRegressResult dataclass

FRegressResult(
    beta: tuple[FData, ...],
    fitted: Any,
    names: tuple[str, ...],
    y: Any,
    terms: tuple[_Term, ...],
    weights: Array,
    cmat: Array,
    dmat: Array,
    df: float | None,
    gcv: float | None,
    ocv: float | None,
    y2c_map: Array | None = None,
)

A fitted functional linear model, as returned by :func:fregress.

Attributes:

Name	Type	Description
beta	tuple of FData	Estimated coefficient function of every term, one curve each (R's betaestlist). The coefficient of a scalar term in a scalar-response model lives in a :class:~fdatools.basis.Constant basis.
fitted	array or FData	Fitted values ŷ: a vector for a scalar response, curves expressed in the response basis for a functional one (R's yhatfdobj).
names	tuple of str	Term names, in model order.
y	array or FData	The response.
terms	tuple	The covariates with their coefficient bases, smoothing parameters and penalty operators, as used in the fit.
weights	array	Observation weights.
cmat, dmat	array	The normal equations C b = D (R's Cmat and Dmat).
df	float or None	Degrees of freedom tr H of the hat matrix (scalar response only).
gcv	float or None	SSE / (n - df)² (scalar response only).
ocv	float or None	Ordinary cross-validation score, the sum of squared leave-one-out residuals (scalar response only).
y2c_map	array or None	Data-to-coefficient map of the response smooth, kept when y was a :class:~fdatools.smoothing.SmoothResult; used by :meth:stderr.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.regression import fregress
>>> basis = fdt.BSpline(domain=(0.0, 1.0), n_basis=6)
>>> rng = np.random.default_rng(1)
>>> y = fdt.FData(rng.standard_normal((6, 12)), basis)
>>> group = np.repeat([0.0, 1.0], 6)
>>> model = fregress(y, {"const": 1.0, "group": group})
>>> model.names
('const', 'group')
>>> model.fitted.n_curves
12

functional_response property

functional_response: bool

Whether the response is a set of curves.

Returns:

Type	Description
bool	True for a functional response, False for a scalar one.

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> fregress(np.arange(5.0), [1.0]).functional_response
False

coefficients property

coefficients: Array

The stacked coefficient vector b solving C b = D.

Returns:

Type	Description
array	Vector of length sum_j n_basis_j.

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> fregress(np.array([1.0, 2.0, 3.0]), [1.0]).coefficients.tolist()
[2.0]

predict

predict(x: Any = None) -> Any

Predict the response for new covariates (R's predict.fRegress).

Parameters:

Name	Type	Description	Default
x	(mapping, sequence or covariate)	New covariates, laid out like those passed to :func:fregress: a mapping keyed by term name, or a sequence in term order. Numbers are broadcast. None returns the fitted values.	None

Returns:

Type	Description
array or FData	Predictions: a vector for a scalar response, curves in the response basis for a functional one.

Raises:

Type	Description
ValueError	If x does not supply one covariate per term, or a covariate has the wrong kind or size.

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> z = np.linspace(0.0, 1.0, 10)
>>> model = fregress(1.0 + 2.0 * z, [1.0, z])
>>> np.round(model.predict([1.0, np.array([0.0, 3.0])]), 8).tolist()
[1.0, 7.0]

stderr

stderr(
    sigma_e: Any = None, y2c_map: Any = None
) -> FRegressStderr

Return the standard errors of the coefficients (R's fRegress.stderr).

The coefficients are linear in the data, b = M y, so their covariance is V = M Σ Mᵀ. For a functional response y enters through the coefficients of its smooth, c_i = S y_i (S the y2c_map of the smoothing step) and the curves are independent, so V = C⁻¹ (Σ_i M_i S Σ Sᵀ M_iᵀ) C⁻¹.

Parameters:

Name	Type	Description	Default
sigma_e	float or array	Residual covariance Σ. A number means σ² I. For a scalar response it defaults to SSE / (n - df) · I; for a functional response it is the (n_t, n_t) covariance of the residuals at the observation points and must be given.	None
y2c_map	array	Map from observations to response coefficients. The identity for a scalar response. For a functional response it defaults to the map stored when y was a :class:~fdatools.smoothing.SmoothResult.	None

Returns:

Type	Description
FRegressStderr	Pointwise standard-error curves, fitted-value standard errors and the coefficient covariance.

Raises:

Type	Description
ValueError	If a functional response lacks sigma_e or y2c_map, or the shapes do not match.

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> z = np.linspace(-1.0, 1.0, 30)
>>> model = fregress(z + np.cos(9.0 * z), [1.0, z])
>>> se = model.stderr(sigma_e=0.25)
>>> bool(np.allclose(np.diag(se.cov), 0.25 * np.diag(np.linalg.inv(model.cmat))))
True

cv

cv() -> FRegressCV

Leave-one-out cross-validation (R's fRegress.CV).

Each observation is removed in turn, the model is refitted on the rest with the same bases and smoothing parameters, and the held-out observation is predicted.

Returns:

Type	Description
FRegressCV	The summed squared errors and the individual errors.

Examples:

>>> import numpy as np
>>> from fdatools.regression import fregress
>>> z = np.linspace(0.0, 1.0, 12)
>>> y = z + 0.1 * np.sin(20.0 * z)
>>> model = fregress(y, [1.0, z])
>>> bool(np.isclose(model.cv().sse, model.ocv))
True

LinmodResult dataclass

LinmodResult(
    alpha: FData,
    beta: BiFData,
    fitted: FData,
    y: FData,
    x: FData,
    weights: Array,
    lam: tuple[float, float, float],
    penalty: tuple[LDO, LDO, LDO],
    cmat: Array,
    dmat: Array,
)

A fitted fully functional linear model, as returned by :func:linmod.

Replaces the list R's linmod returns.

Attributes:

Name	Type	Description
alpha	FData	The intercept function alpha(t), one curve (R's beta0estfd).
beta	BiFData	The regression surface β(s, t) on the tensor basis θ_s ⊗ θ_t (R's beta1estbifd).
fitted	FData	Fitted curves ŷ_i(t) = alpha(t) + ∫ x_i(s) β(s, t) ds, projected onto the response basis (R's yhatfdobj).
y	FData	The response curves.
x	FData	The covariate curves.
weights	array	Observation weights.
lam	tuple of float	Smoothing parameters (λ_alpha, λ_s, λ_t).
penalty	tuple of LDO	Roughness operators (L_alpha, L_s, L_t).
cmat, dmat	array	The normal equations C b = D; b stacks the intercept coefficients and the surface coefficients row by row.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.regression import linmod
>>> basis = fdt.BSpline(domain=(0.0, 1.0), n_basis=6)
>>> rng = np.random.default_rng(0)
>>> x = fdt.FData(rng.standard_normal((6, 20)), basis)
>>> y = fdt.FData(rng.standard_normal((6, 20)), basis)
>>> model = linmod(y, x, lam_s=1e-4, lam_t=1e-4)
>>> model.beta.coefs.shape
(6, 6)
>>> model.residuals.n_curves
20

residuals property

residuals: FData

The residual curves y_i - ŷ_i in the response basis.

Returns:

Type	Description
FData	One residual curve per observation.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.regression import linmod
>>> basis = fdt.BSpline(domain=(0.0, 1.0), n_basis=5)
>>> rng = np.random.default_rng(2)
>>> x = fdt.FData(rng.standard_normal((5, 15)), basis)
>>> y = fdt.FData(rng.standard_normal((5, 15)), basis)
>>> model = linmod(y, x, lam_s=1e-3, lam_t=1e-3)
>>> bool(np.allclose((model.fitted + model.residuals).coefs, y.coefs))
True

predict

predict(x: Any = None) -> FData

Predict response curves for new covariate curves.

R has no predict method for linmod; this evaluates the fitted model alpha(t) + ∫ x(s) β(s, t) ds for the new curves.

Parameters:

Name	Type	Description	Default
x	FData or SmoothResult	New covariate curves on the domain of s, in any basis. None returns the fitted values.	None

Returns:

Type	Description
FData	Predicted curves in the response basis, one per new curve.

Raises:

Type	Description
TypeError	If x is not a set of curves.
ValueError	If x lives on another domain or carries a variable axis.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.regression import linmod
>>> basis = fdt.BSpline(domain=(0.0, 1.0), n_basis=5)
>>> rng = np.random.default_rng(1)
>>> x = fdt.FData(rng.standard_normal((5, 25)), basis)
>>> y = fdt.FData(rng.standard_normal((5, 25)), basis)
>>> model = linmod(y, x, lam_s=1e-3, lam_t=1e-3)
>>> model.predict(x[:3]).n_curves
3

FRegress

FRegress(
    beta: Any = None,
    *,
    lam: float = 0.0,
    penalty: int | LDO = 2,
    fit_intercept: bool = True,
)

Bases: RegressorMixin, BaseEstimator

Functional linear regression as a scikit-learn regressor.

X is either a numeric matrix of scalar covariates (one column per covariate), a single :class:FData covariate, or a list of covariates as accepted by :func:fregress. y is a vector (scalar response) or an :class:FData (functional response).

Parameters:

Name	Type	Description	Default
beta	mapping or sequence	Coefficient specification per covariate, as in :func:fregress. With fit_intercept the intercept is the first term, named const.	None
lam	float	Default smoothing parameter.	0.0
penalty	int or LDO	Default roughness operator.	2
fit_intercept	bool	Prepend the constant covariate 1. Default True.	True

Attributes:

Name	Type	Description
result_	FRegressResult	The fitted model.
coef_	tuple of FData	The estimated coefficient functions (result_.beta).
n_features_in_	int	Number of columns seen during fit when X was a matrix.

Examples:

>>> import numpy as np
>>> from fdatools.regression import FRegress
>>> rng = np.random.default_rng(0)
>>> X = rng.standard_normal((50, 2))
>>> y = 0.5 + X @ np.array([1.0, -2.0])
>>> model = FRegress().fit(X, y)
>>> np.round(model.predict(X[:2]) - y[:2], 10).tolist()
[0.0, 0.0]

fit

fit(X: Any, y: Any) -> FRegress

Fit the model.

Parameters:

Name	Type	Description	Default
X	array of shape (n_samples, n_features), FData or list	Covariates.	required
y	array of shape (n_samples,), FData or SmoothResult	Response. A numeric y must be finite; a column vector is flattened with a DataConversionWarning.	required

Returns:

Type	Description
FRegress	The fitted estimator.

Raises:

Type	Description
ValueError	If y is missing or not finite, X and y hold different numbers of samples, or there is only one sample.

Examples:

>>> import numpy as np
>>> from fdatools.regression import FRegress
>>> X = np.array([[0.0], [1.0], [2.0], [3.0]])
>>> model = FRegress().fit(X, 1.0 + 2.0 * X[:, 0])
>>> [round(float(b.coefs[0, 0]), 10) for b in model.coef_]
[1.0, 2.0]

predict

predict(X: Any) -> Any

Predict the response for new covariates.

Parameters:

Name	Type	Description	Default
X	array of shape (n_samples, n_features), FData or list	Covariates laid out as in :meth:fit.	required

Returns:

Type	Description
array or FData	Predicted response.

fregress

fregress(
    y: Any,
    x: Any = None,
    beta: Any = None,
    *,
    lam: float = 0.0,
    penalty: int | LDO = 2,
    weights: Any = None,
) -> FRegressResult

Fit a functional linear model; the model type follows from the arguments.

Replaces R's fRegress (including fRegress.formula, fRegress.fd and fRegress.double).

Parameters:

Name	Type	Description	Default
y	(array, FData, SmoothResult or str)	The response: n numbers (scalar response), an :class:FData of n curves (functional response), or a :class:~fdatools.smoothing.SmoothResult whose y2c_map is then kept for :meth:FRegressResult.stderr. A string is a formula "response ~ a + b", evaluated against the mapping x: the intercept const is included unless the formula has - 1 or + 0, and a covariate given as labels (strings) is expanded into treatment-coded indicators name.level, first sorted level dropped, as in R.	required
x	(mapping, sequence or covariate)	The covariates, each either an :class:FData of n curves, n numbers, or a single number broadcast to all observations (1.0 is an intercept). A mapping names the terms; a sequence names them x0, x1, .... The data mapping of a formula.	None
beta	mapping or sequence	Coefficient specification per term (by name or position): None, a :class:~fdatools.basis.Basis, (basis, lam) or (basis, lam, penalty). A mapping may also be keyed by the variable a formula term came from. The default basis is the response basis for a functional response, and for a scalar response the covariate's own basis (curves) or a :class:~fdatools.basis.Constant (numbers).	None
lam	float	Smoothing parameter applied to every term that does not give its own. Default 0.	0.0
penalty	int or LDO	Roughness operator for every term that does not give its own; an integer means D^penalty. Default 2.	2
weights	array	Positive observation weights, one per observation. The fit is then weighted least squares; see Notes.	None

Returns:

Type	Description
FRegressResult	Coefficient functions, fitted values and fit diagnostics, with :meth:~FRegressResult.predict, :meth:~FRegressResult.stderr and :meth:~FRegressResult.cv.

Raises:

Type	Description
ValueError	If the covariates disagree with the response in size or domain, a scalar-on-scalar term is given a non-constant basis, or the formula is malformed.
KeyError	If a formula names a variable missing from the data.

Notes

With weights the coefficients minimise the penalised weighted criterion :math:\sum_i w_i \|y_i - \hat y_i\|^2 + \text{penalties}, for a scalar and a functional response alike. This is what R's fRegress(..., wt = w) computes, and the two agree to the accuracy of R's numerical integration. R has a trap here: its weight argument is called wt, and fRegress(..., wtvec = w) (the name R's smoothing functions use) is swallowed by ... without a warning, so R returns the unweighted fit. Scripts that pass wtvec to fRegress must drop weights when ported, or they will get a different answer in fdatools. Only the relative weights matter when there is no penalty; with a penalty, scaling every weight by c acts like dividing lam by c. R's Fperm.fd ignores weights, while :func:fdatools.stats.f_test on a weighted model does not -- see its Notes.

Examples:

Scalar response on a functional covariate:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.regression import fregress
>>> rng = np.random.default_rng(3)
>>> basis = fdt.BSpline(domain=(0.0, 1.0), n_basis=7)
>>> x = fdt.FData(rng.standard_normal((7, 40)), basis)
>>> y = np.asarray(fdt.inprod(x, fdt.FData(np.linspace(-1, 1, 7), basis)))[:, 0]
>>> model = fregress(y, [1.0, x])
>>> np.round(model.beta[1].coefs[:, 0], 6).tolist()
[-1.0, -0.666667, -0.333333, 0.0, 0.333333, 0.666667, 1.0]

The formula interface with a categorical covariate:

>>> region = ["north", "south", "north", "east", "south", "east"]
>>> yf = fdt.FData(rng.standard_normal((7, 6)), basis)
>>> fregress("temp ~ region", {"temp": yf, "region": region}).names
('const', 'region.north', 'region.south')

linmod

linmod(
    y: Any,
    x: Any,
    *,
    alpha_basis: Basis | None = None,
    s_basis: Basis | None = None,
    t_basis: Basis | None = None,
    lam_alpha: float = 0.0,
    lam_s: float = 0.0,
    lam_t: float = 0.0,
    penalty_alpha: int | LDO = 2,
    penalty_s: int | LDO = 2,
    penalty_t: int | LDO = 2,
    weights: Any = None,
) -> LinmodResult

Fit the fully functional linear model with a bivariate coefficient.

Replaces R's linmod. Each response curve is explained by a whole covariate curve:

.. math:: y_i(t) = \alpha(t) + \int x_i(s)\,\beta(s, t)\,ds + e_i(t) ,

where s and t may live on different intervals. The intercept is alpha = φ_alphaᵀ a and the surface is β(s, t) = θ_s(s)ᵀ B θ_t(t); a and B minimise

.. math::

\sum_i w_i \int (y_i - \hat y_i)^2\,dt
+ \lambda_\alpha \int (L_\alpha \alpha)^2\,dt
+ \lambda_s \iint (L_s \beta)^2\,ds\,dt
+ \lambda_t \iint (L_t \beta)^2\,ds\,dt ,

L_s acting on s and L_t on t. This is R's linmod with betaList = list(fdPar(alpha_basis, penalty_alpha, lam_alpha), bifdPar(bifd(0, s_basis, t_basis), penalty_s, penalty_t, lam_s, lam_t)).

Parameters:

Name	Type	Description	Default
y	FData or SmoothResult	The n response curves.	required
x	FData or SmoothResult	The n covariate curves.	required
alpha_basis	Basis	Basis of the intercept alpha(t). Defaults to the response basis.	None
s_basis	Basis	Basis of β in s (the covariate's argument). Defaults to the covariate basis.	None
t_basis	Basis	Basis of β in t (the response's argument). Defaults to the response basis.	None
lam_alpha	float	Smoothing parameters of the intercept and of the surface in s and in t (R's lambda, lambdas and lambdat). Default 0.	0.0
lam_s	float	Smoothing parameters of the intercept and of the surface in s and in t (R's lambda, lambdas and lambdat). Default 0.	0.0
lam_t	float	Smoothing parameters of the intercept and of the surface in s and in t (R's lambda, lambdas and lambdat). Default 0.	0.0
penalty_alpha	int or LDO	Roughness operators matching the three smoothing parameters; an integer means D^penalty. Default 2.	2
penalty_s	int or LDO	Roughness operators matching the three smoothing parameters; an integer means D^penalty. Default 2.	2
penalty_t	int or LDO	Roughness operators matching the three smoothing parameters; an integer means D^penalty. Default 2.	2
weights	array	Positive observation weights, one per curve. Default: all 1.	None

Returns:

Type	Description
LinmodResult	The intercept, the surface, the fitted curves and the normal equations, with :meth:~LinmodResult.predict.

Raises:

Type	Description
TypeError	If x or y is not a set of curves.
ValueError	If x and y hold different numbers of curves, a basis lives on the wrong domain, a smoothing parameter is negative or not finite, or a weight is not positive.

Notes

Every integral is exact (Gauss-Legendre on the break points of the bases involved) and the fitted curves are the exact L2 projection onto the response basis. R's linmod integrates numerically and fits ŷ by least squares on 201 points, so the two agree to the accuracy of R's quadrature: to rounding when R's integrals are exact (polynomial bases), about 1e-6 on the weather data, 1e-4 on cubic B-splines. R's linmod(..., wtvec = w) stops with an error in fda 6.3.0; fdatools' weights give the weighted least-squares fit.

PyTorch coefficients (in x, y or weights) make the whole fit run in PyTorch, so gradients flow back to the inputs.

Examples:

A surface that the data determine exactly is recovered:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.regression import linmod
>>> sbasis = fdt.BSpline(domain=(0.0, 1.0), n_basis=5)
>>> tbasis = fdt.BSpline(domain=(0.0, 2.0), n_basis=4)
>>> rng = np.random.default_rng(3)
>>> x = fdt.FData(rng.standard_normal((5, 30)), sbasis)
>>> surface = rng.standard_normal((5, 4))
>>> z = np.asarray(fdt.inprod(x, sbasis))
>>> y = fdt.FData(1.0 + surface.T @ z.T, tbasis)
>>> model = linmod(y, x)
>>> bool(np.allclose(model.beta.coefs, surface))
True
>>> np.round(model.alpha(np.array([0.5, 1.5]))[:, 0], 8).tolist()
[1.0, 1.0]