Sparse data (PACE)¶

fdatools.sparse: functional PCA for sparse or irregular longitudinal data (PACE, Yao, Müller & Wang 2005). sparse_mean fits the mean from pooled points, sparse_cov fits the covariance surface from within-curve cross-products, and the PACE estimator gives harmonics, eigenvalues and conditional-expectation scores.

Functional PCA for sparse, irregular longitudinal data (PACE).

When every curve is seen at only a handful of irregular times, smoothing each curve on its own fails: a curve with four points cannot carry a ten-function basis. PACE -- principal components analysis through conditional expectation (Yao, Müller & Wang, 2005) -- pools the data of all curves instead. It replaces R fda's smooth.sparse.mean, covPACE, pcaPACE and scoresPACE.

The input is the irregular per-curve form :func:fdatools.smoothing.smooth accepts: a sequence t of observation times, one array per curve, and a matching sequence y of observed values.

Mean. One penalised least-squares fit to the pooled points (t_ij, y_ij):

.. math:: \min_c \sum_{ij} (y_{ij} - \varphi(t_{ij})^{T} c)^2 + \lambda_\mu\, c^{T} R\, c .

Covariance. Every pair of distinct observations of the same curve gives a raw covariance r_ij r_ik at (t_ij, t_ik), with r_ij = y_ij - μ(t_ij). A tensor-product surface G(s, t) = ψ(s)ᵀ C ψ(t) is fitted to them by penalised least squares, with the roughness penalty λ_G vec(C)ᵀ (W ⊗ P + P ⊗ W) vec(C) (W the Gram matrix, P the roughness matrix of ψ). The diagonal pairs j = k are left out: they carry the measurement-error variance σ² on top of G(t, t).
Error variance. The left-out diagonal r_ij² is smoothed on its own into V(t), and σ² is the mean gap between V and G over the middle half of the domain, where both are best determined:

.. math:: \sigma^2 = \frac{2}{|T|} \int_{T_1} {V(t) - G(t, t)}\,dt , \qquad T_1 = [a + |T|/4,\; b - |T|/4] .

Eigenfunctions. The harmonics ξ = θ(t)ᵀ b in a basis θ solve the penalised eigenproblem of :class:~fdatools.decomposition.FPCA with the covariance operator of G:

.. math:: J C J^{T} b = \mu\,(W_\theta + \lambda R_\theta)\,b , \qquad J = \int \theta\,\psi^{T} .

Scores. The best linear unbiased predictor of the scores of curve i given its own observations is the conditional expectation

.. math:: \hat\xi_i = \Lambda \Xi_i^{T} \Sigma_i^{-1} (y_i - \mu_i), \qquad \Sigma_i = \Xi_i \Lambda \Xi_i^{T} + \sigma^2 I ,

where Ξ_i holds the harmonics at the times of curve i and Λ the retained eigenvalues. It works with any number of points per curve.

Every integral is exact (Gauss-Legendre on the basis break points), where R uses the Romberg approximation of inprod(); see tests/parity/test_pace.py for the measured consequences. R's scoresPACE does not compute the conditional expectation above (it evaluates the harmonics at a single, wrongly indexed point per curve); fdatools implements the published estimator.

References

Yao, F., Müller, H.-G. and Wang, J.-L. (2005). Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association, 100, 577-590.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.sparse import PACE
>>> rng = np.random.default_rng(0)
>>> t = [np.sort(rng.uniform(0.0, 1.0, 5)) for _ in range(200)]
>>> a = rng.normal(0.0, 1.0, 200)
>>> y = [ai * np.sin(np.pi * ti) + 0.1 * rng.normal(size=5) for ai, ti in zip(a, t)]
>>> pace = PACE(n=1, basis=BSpline(domain=(0.0, 1.0), n_basis=4)).fit(y, t=t)
>>> pace.scores.shape
(200, 1)
>>> bool(abs(np.corrcoef(pace.scores[:, 0], a)[0, 1]) > 0.95)
True

SparseCov `dataclass` ¶

SparseCov(
    cov: BiFData,
    mean: FData,
    sigma2: float,
    variance: FData,
)

Covariance estimate of sparse longitudinal data.

Returned by :func:sparse_cov; replaces the list R's covPACE returns.

Attributes:

Name	Type	Description
`cov`	`BiFData`	The smoothed covariance surface `G(s, t)`, fitted to the off-diagonal raw covariances only.
`mean`	`FData`	The mean function the residuals were taken from.
`sigma2`	`float`	Estimated measurement-error variance `σ²`.
`variance`	`FData`	Smooth `V(t)` of the diagonal raw covariances `r_ij²`, i.e. the variance of an observation, `G(t, t) + σ²`.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.sparse import sparse_cov
>>> rng = np.random.default_rng(2)
>>> t = [np.sort(rng.uniform(0.0, 1.0, 6)) for _ in range(40)]
>>> y = [rng.normal() + 0.2 * rng.normal(size=6) for _ in t]
>>> est = sparse_cov(y, t, basis=BSpline(domain=(0.0, 1.0), n_basis=3, order=2))
>>> est.cov.coefs.shape
(3, 3)

PACE ¶

PACE(
    n: SupportsIndex = 2,
    *,
    basis: Basis | None = None,
    mean_basis: Basis | None = None,
    harmonic_basis: Basis | None = None,
    lam_mean: float = 0.0,
    lam_cov: float = 0.0,
    lam: float = 0.0,
    penalty: int | LDO = 2,
    sigma2: float | None = None,
)

Bases: TransformerMixin, BaseEstimator

Functional PCA for sparse longitudinal data by conditional expectation.

Replaces R's pcaPACE and scoresPACE (and runs smooth.sparse.mean and covPACE on the way); see :mod:fdatools.sparse for the five steps.

Parameters:

Name	Type	Description	Default
`n`	`int or SupportsIndex`	Number of harmonics to keep. Default `2`.	`2`
`basis`	`Basis`	Basis of the covariance surface. Default: a cubic B-spline over the observed range with `min(6, number of distinct times)` functions.	`None`
`mean_basis`	`Basis`	Basis of the mean. Default: as `basis` but with up to 10 functions.	`None`
`harmonic_basis`	`Basis`	Basis of the harmonics. Default: the covariance basis.	`None`
`lam_mean`	`float`	Roughness penalties on the mean, the covariance surface and the harmonics. All default to `0.0`.	`0.0`
`lam_cov`	`float`	Roughness penalties on the mean, the covariance surface and the harmonics. All default to `0.0`.	`0.0`
`lam`	`float`	Roughness penalties on the mean, the covariance surface and the harmonics. All default to `0.0`.	`0.0`
`penalty`	`int or LDO`	Roughness operator used by all three penalties. Default `2`.	`2`
`sigma2`	`float`	Measurement-error variance to use in the scores. `None` (default) estimates it from the data; a given value must be positive.	`None`

Attributes:

Name	Type	Description
`mean_fd`	`FData`	The mean function.
`cov`	`BiFData`	The covariance surface, fitted without the diagonal.
`sigma2_`	`float`	The measurement-error variance used for the scores.
`harmonics`	`FData`	The `n` eigenfunctions, `(W + λR)`-orthonormal. Each is signed so that its coefficients in `(W + λR)`-orthonormal coordinates (`U b` with `W + λR = UᵀU`) sum to a positive number, R's rule.
`values`	`array`	The `n` leading eigenvalues of the covariance operator.
`varprop`	`array`	`values / values.sum()`, as R's `pcaPACE` reports it.
`scores`	`array`	Conditional-expectation scores of the fitted curves, `(n_curves, n)`.

Notes

fit and transform take ragged input (a different number of points per curve), which scikit-learn's generic estimator checks cannot generate, so :func:sklearn.utils.estimator_checks.check_estimator does not apply. Parameters, clone and pickling follow the estimator API. Negative eigenvalues of an indefinite surface estimate are kept in values but count as zero in the scores.

σ² is the difference of two estimated surfaces on their diagonal, so with few curves it is poorly determined and can come out negative; it is then floored (with a :class:RuntimeWarning) and the scores approach the least-squares fit of the harmonics to each curve. Pass sigma2 when the measurement error is known, or use more curves or smoother surfaces. With sigma2 given, the estimate is still kept in cov_estimate_ but is not used, so no warning is raised.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.sparse import PACE
>>> rng = np.random.default_rng(1)
>>> t = [np.sort(rng.uniform(0.0, 1.0, 4)) for _ in range(150)]
>>> y = [1.0 + rng.normal() * np.cos(np.pi * ti) + 0.05 * rng.normal(size=4) for ti in t]
>>> pace = PACE(n=2, basis=BSpline(domain=(0.0, 1.0), n_basis=4)).fit(y, t=t)
>>> pace.harmonics.n_curves
2
>>> bool(pace.varprop[0] > 0.8)
True

mean_fd `property` ¶

mean_fd: FData

The mean function.

cov `property` ¶

cov: BiFData

The covariance surface.

harmonics `property` ¶

harmonics: FData

The retained eigenfunctions.

values `property` ¶

values: Array

The retained eigenvalues, descending.

varprop `property` ¶

varprop: Array

Share of each retained eigenvalue in their sum.

scores `property` ¶

scores: Array

Conditional-expectation scores of the fitted curves.

fit ¶

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

Estimate mean, covariance, harmonics and scores from sparse curves.

Parameters:

Name	Type	Description	Default
`X`	`sequence of array_like`	Observed values, one 1-D array per curve; or, with `t=None`, R's list form of `(n_i, 2)` `(time, value)` arrays.	required
`y`	`ignored`	Present for API compatibility.	`None`
`t`	`sequence of array_like`	Observation times, one 1-D array per curve.	`None`

Returns:

Type	Description
`PACE`	The fitted estimator.

Raises:

Type	Description
`ValueError`	If the input or a parameter is invalid, or a basis is not determined by the data.

Examples:

>>> import numpy as np
>>> from fdatools.sparse import PACE
>>> rng = np.random.default_rng(5)
>>> t = [np.sort(rng.uniform(0.0, 10.0, 5)) for _ in range(30)]
>>> y = [rng.normal() * ti / 10.0 + 0.1 * rng.normal(size=5) for ti in t]
>>> PACE(n=1).fit(y, t=t).values.shape
(1,)

transform ¶

transform(X: Any, t: Any = None) -> Array

Return the conditional-expectation scores of new sparse curves.

Parameters:

Name	Type	Description	Default
`X`	`sequence of array_like`	Observed values, one 1-D array per curve (or R's list form when `t` is `None`).	required
`t`	`sequence of array_like`	Observation times, one 1-D array per curve.	`None`

Returns:

Type	Description
`array`	Scores of shape `(n_curves, n)`.

Examples:

>>> import numpy as np
>>> from fdatools.sparse import PACE
>>> rng = np.random.default_rng(6)
>>> t = [np.sort(rng.uniform(0.0, 1.0, 5)) for _ in range(40)]
>>> y = [rng.normal() * np.sin(np.pi * ti) + 0.1 * rng.normal(size=5) for ti in t]
>>> pace = PACE(n=2).fit(y, t=t)
>>> bool(np.allclose(pace.transform(y, t), pace.scores))
True

fit_transform ¶

fit_transform(
    X: Any, y: Any = None, *, t: Any = None
) -> Array

Fit to X and return the scores of the same curves.

Parameters:

Name	Type	Description	Default
`X`	`sequence of array_like`	Observed values, one array per curve.	required
`y`	`ignored`	Present for API compatibility.	`None`
`t`	`sequence of array_like`	Observation times, one array per curve.	`None`

Returns:

Type	Description
`array`	Scores of shape `(n_curves, n)`.

Examples:

>>> import numpy as np
>>> from fdatools.sparse import PACE
>>> rng = np.random.default_rng(7)
>>> t = [np.sort(rng.uniform(0.0, 1.0, 5)) for _ in range(40)]
>>> y = [rng.normal() + 0.1 * rng.normal(size=5) for ti in t]
>>> PACE(n=1).fit_transform(y, t=t).shape
(40, 1)

inverse_transform ¶

inverse_transform(X: Any, t: Any) -> Array

Evaluate the curves rebuilt from scores, μ(t) + Σ_k ξ_k φ_k(t).

Parameters:

Name	Type	Description	Default
`X`	`array_like of shape (n_samples, n)`	Scores, e.g. from :meth:`transform`.	required
`t`	`array_like`	Evaluation points (one shared grid).	required

Returns:

Type	Description
`array`	Values of shape `(n_samples, len(t))`.

Raises:

Type	Description
`ValueError`	If `X` does not have one column per harmonic.

Examples:

>>> import numpy as np
>>> from fdatools.sparse import PACE
>>> rng = np.random.default_rng(8)
>>> t = [np.sort(rng.uniform(0.0, 1.0, 6)) for _ in range(50)]
>>> y = [rng.normal() * ti + 0.05 * rng.normal(size=6) for ti in t]
>>> pace = PACE(n=1).fit(y, t=t)
>>> pace.inverse_transform(pace.scores, np.linspace(0.1, 0.9, 11)).shape
(50, 11)

sparse_mean ¶

sparse_mean(
    y: Any,
    t: Any = None,
    basis: Basis | None = None,
    *,
    lam: float = 0.0,
    penalty: int | LDO = 2,
) -> FData

Estimate the mean function of sparse, irregularly observed curves.

Replaces R's smooth.sparse.mean. All observations are pooled into one penalised least-squares fit

.. math:: (\Phi^{T}\Phi + \lambda R)\,c = \Phi^{T} y ,

where Φ is the basis evaluated at every observation time of every curve.

Parameters:

Name	Type	Description	Default
`y`	`sequence of array_like`	Observed values, one 1-D array per curve. With `t=None`, one `(n_i, 2)` array per curve holding `(time, value)` rows instead (R's list form).	required
`t`	`sequence of array_like`	Observation times, one 1-D array per curve, matching `y`.	`None`
`basis`	`Basis`	Basis of the mean. Default: a cubic B-spline over the observed range with `min(10, number of distinct times)` functions.	`None`
`lam`	`float`	Roughness penalty. Default `0.0`.	`0.0`
`penalty`	`int or LDO`	Roughness operator. Default `2` (the second derivative), as in R.	`2`

Returns:

Type	Description
`FData`	The mean as a single curve (coefficients of shape `(n_basis, 1)`).

Raises:

Type	Description
`ValueError`	If the input is malformed, `lam` is negative, or the basis is not determined by the pooled data.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.sparse import sparse_mean
>>> t = [np.array([0.0, 0.4, 0.9]), np.array([0.1, 0.5, 1.0]), np.array([0.2, 0.7])]
>>> y = [2.0 * ti + 1.0 for ti in t]
>>> mean = sparse_mean(y, t, BSpline(domain=(0.0, 1.0), n_basis=4))
>>> np.round(mean(np.array([0.0, 0.5, 1.0]))[:, 0], 10).tolist()
[1.0, 2.0, 3.0]

sparse_cov ¶

sparse_cov(
    y: Any,
    t: Any = None,
    *,
    mean: FData | None = None,
    basis: Basis | None = None,
    lam: float = 0.0,
    penalty: int | LDO = 2,
) -> SparseCov

Estimate the covariance surface of sparse, irregularly observed curves.

Replaces R's covPACE. The raw covariances r_ij r_ik of every pair of distinct observations of the same curve are fitted by the tensor product surface G(s, t) = ψ(s)ᵀ C ψ(t), minimising

.. math:: \sum_i \sum_{j \ne k} {r_{ij} r_{ik} - \psi(t_{ij})^{T} C\, \psi(t_{ik})}^2 + \lambda\, \mathrm{vec}(C)^{T} (W \otimes P + P \otimes W)\, \mathrm{vec}(C) .

The diagonal j = k is excluded because E r_ij² = G(t, t) + σ²; it is smoothed separately and gives the measurement-error variance σ² (see :mod:fdatools.sparse).

Parameters:

Name	Type	Description	Default
`y`	`sequence of array_like`	Observed values, one 1-D array per curve, or R's list form of `(n_i, 2)` `(time, value)` arrays when `t` is `None`.	required
`t`	`sequence of array_like`	Observation times, one 1-D array per curve.	`None`
`mean`	`FData`	Mean function to subtract. Default: :func:`sparse_mean` of the data with its default basis.	`None`
`basis`	`Basis`	Basis `ψ` of both arguments of the surface. Default: a cubic B-spline over the observed range with `min(6, number of distinct times)` functions.	`None`
`lam`	`float`	Roughness penalty. Default `0.0`.	`0.0`
`penalty`	`int or LDO`	Roughness operator. Default `2`.	`2`

Returns:

Type	Description
`SparseCov`	The surface, the mean, `σ²` and the smoothed diagonal.

Raises:

Type	Description
`ValueError`	If the input is malformed, no curve has two observations, or the surface is not determined by the data (use a smaller basis or `lam > 0`).

Warns:

Type	Description
`RuntimeWarning`	If the `σ²` estimate is not positive and is replaced by a small floor.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline
>>> from fdatools.sparse import sparse_cov
>>> rng = np.random.default_rng(3)
>>> t = [np.sort(rng.uniform(0.0, 1.0, 5)) for _ in range(200)]
>>> y = [rng.normal(0.0, 2.0) + rng.normal(0.0, 0.5, 5) for _ in t]
>>> est = sparse_cov(y, t, basis=BSpline(domain=(0.0, 1.0), n_basis=1, order=1))
>>> bool(abs(float(est.cov.coefs[0, 0]) - 4.0) < 1.0)
True
>>> bool(abs(est.sigma2 - 0.25) < 0.1)
True

Sparse data (PACE)¶

SparseCov dataclass ¶

PACE ¶

mean_fd property ¶

cov property ¶

harmonics property ¶

values property ¶

varprop property ¶

scores property ¶

fit ¶

transform ¶

fit_transform ¶

inverse_transform ¶

sparse_mean ¶

sparse_cov ¶

SparseCov `dataclass` ¶

mean_fd `property` ¶

cov `property` ¶

harmonics `property` ¶

values `property` ¶

varprop `property` ¶

scores `property` ¶