Add standard deviation as a function

skfda/_utils/__init__.py

@@ -20,6 +20,7 @@
             "_same_domain",
             "_to_grid",
             "_to_grid_points",
+            "function_to_fdatabasis",
             "nquad_vec",
         ],
         '_warping': [

skfda/_utils/_utils.py

@@ -4,7 +4,6 @@
 
 import functools
 import numbers
-from functools import singledispatch
 from typing import (
     TYPE_CHECKING,
     Any,
@@ -36,7 +35,7 @@
 ArrayDTypeT = TypeVar("ArrayDTypeT", bound="np.generic")
 
 if TYPE_CHECKING:
-    from ..representation import FData, FDataGrid
+    from ..representation import FData, FDataBasis, FDataGrid
     from ..representation.basis import Basis
     from ..representation.extrapolation import ExtrapolationLike
 
@@ -609,3 +608,36 @@ def _classifier_get_classes(
             f'one; got {classes.size} class',
         )
     return classes, y_ind
+
+
+def function_to_fdatabasis(
+    f: Callable[[NDArrayFloat], NDArrayFloat],
+    new_basis: Basis,
+) -> FDataBasis:
+    """Express a math function as a FDataBasis with a given basis.
+
+    Args:
+        f: math function.
+        new_basis: the basis of the output.
+
+    Returns:
+        FDataBasis: FDataBasis with calculated coefficients and the new
+        basis.
+    """
+    from .. import FDataBasis
+    from ..misc._math import inner_product_matrix
+
+    if isinstance(f, FDataBasis) and f.basis == new_basis:
+        return f
+
+    inner_prod = inner_product_matrix(
+        new_basis,
+        f,
+        _domain_range=new_basis.domain_range,
+    )
+
+    gram_matrix = new_basis.gram_matrix()
+
+    coefs = np.linalg.solve(gram_matrix, inner_prod)
+
+    return FDataBasis(new_basis, coefs.T)

skfda/exploratory/stats/__init__.py

@@ -19,6 +19,7 @@
             "gmean",
             "mean",
             "modified_epigraph_index",
+            "std",
             "trim_mean",
             "var",
         ],
@@ -37,6 +38,7 @@
         gmean as gmean,
         mean as mean,
         modified_epigraph_index as modified_epigraph_index,
+        std as std,
         trim_mean as trim_mean,
         var as var,
     )

skfda/exploratory/stats/_stats.py

@@ -4,12 +4,13 @@
 from builtins import isinstance
 from typing import Callable, TypeVar, Union
 
+import functools
 import numpy as np
 from scipy import integrate
 from scipy.stats import rankdata
 
 from ...misc.metrics._lp_distances import l2_distance
-from ...representation import FData, FDataGrid
+from ...representation import FData, FDataBasis, FDataGrid
 from ...typing._metric import Metric
 from ...typing._numpy import NDArrayFloat
 from ..depth import Depth, ModifiedBandDepth
@@ -76,7 +77,7 @@ def gmean(X: FDataGrid) -> FDataGrid:
 
 def cov(
     X: FData,
-    ddof: int = 1
+    ddof: int = 1,
 ) -> Callable[[NDArrayFloat, NDArrayFloat], NDArrayFloat]:
     """
     Compute the covariance.
@@ -99,6 +100,66 @@ def cov(
     return X.cov(ddof=ddof)
 
 
+@functools.singledispatch
+def std(X: F, ddof: int = 1) -> F:
+    r"""
+    Compute the standard deviation of all the samples in a FData object.
+
+    .. math::
+        \text{std}_X(t) = \sqrt{\frac{1}{N-\text{ddof}}
+        \sum_{n=1}^{N}{\left(X_n(t) - \overline{X}(t)\right)^2}}
+
+    Args:
+        X: Object containing all the samples whose standard deviation is
+            wanted.
+        ddof: Means "Delta Degrees of Freedom". The divisor used in
+            calculations is `N - ddof`, where `N` represents the number of
+            samples in `X`. By default ddof is 1.
+
+    Returns:
+        Standard deviation of all the samples in the original object, as a
+        :term:`functional data object` with just one sample.
+
+    """
+    raise NotImplementedError("Not implemented for this type")
+
+
+@std.register(FDataGrid)
+def std_fdatagrid(X: FDataGrid, ddof: int = 1) -> FDataGrid:
+    """Standard deviation of a FDataGrid."""
+    return X.copy(
+        data_matrix=np.std(X.data_matrix, axis=0, ddof=ddof)[np.newaxis, ...],
+        sample_names=("standard deviation",),
+    )
+
+
+@std.register(FDataBasis)
+def std_fdatabasis(X: FDataBasis, ddof: int = 1) -> FDataBasis:
+    """Standard deviation of a FDataBasis."""
+    from ..._utils import function_to_fdatabasis
+
+    if X.dim_domain != 1 or X.dim_codomain != 1:
+        raise NotImplementedError(
+            "Standard deviation only implemented "
+            "for univariate functions.",
+        )
+
+    basis = X.basis
+    coeff_matrix = np.cov(X.coefficients, rowvar=False, ddof=ddof)
+
+    def std_function(t_points: NDArrayFloat) -> NDArrayFloat:
+        assert len(t_points) == 1, (
+            "Standard deviation function only implemented for "
+            "one-point-at-a-time evaluations."
+        )
+        basis_evaluation = basis(t_points).reshape((-1, 1))
+        return np.sqrt(
+            basis_evaluation.T @ coeff_matrix @ basis_evaluation,
+        ).reshape((1, -1, 1))
+
+    return function_to_fdatabasis(f=std_function, new_basis=X.basis)
+
+
 def modified_epigraph_index(X: FDataGrid) -> NDArrayFloat:
     """
     Calculate the Modified Epigraph Index of a FDataGrid.

skfda/tests/test_stats.py

@@ -2,8 +2,15 @@
 
 import numpy as np
 
-from skfda.datasets import fetch_phoneme, fetch_tecator, fetch_weather
-from skfda.exploratory.stats import geometric_median, modified_epigraph_index
+from skfda.datasets import (
+    fetch_phoneme,
+    fetch_tecator,
+    fetch_weather,
+)
+from skfda.exploratory.stats import (
+    geometric_median,
+    modified_epigraph_index,
+)
 
 
 class TestGeometricMedian(unittest.TestCase):

skfda/tests/test_stats_std.py

@@ -0,0 +1,83 @@
+"""Test stats functions."""
+
+from __future__ import annotations
+
+from typing import Any
+
+import numpy as np
+import pytest
+
+from skfda import FDataGrid
+from skfda.datasets import make_gaussian_process
+from skfda.exploratory.stats import std
+from skfda.misc.covariances import Gaussian
+from skfda.representation.basis import FourierBasis
+
+
+@pytest.fixture(params=[61, 71])
+def n_basis(request: Any) -> int:
+    """Fixture for n_basis to test."""
+    return request.param
+
+
+@pytest.fixture
+def start() -> int:
+    """Fixture for the infimum of the domain."""
+    return 0
+
+
+@pytest.fixture
+def stop() -> int:
+    """Fixture for the supremum of the domain."""
+    return 1
+
+
+@pytest.fixture
+def n_features() -> int:
+    """Fixture for the number of features."""
+    return 1000
+
+
+@pytest.fixture
+def gaussian_process(start: int, stop: int, n_features: int) -> FDataGrid:
+    """Fixture for a Gaussian process."""
+    return make_gaussian_process(
+        start=start,
+        stop=stop,
+        n_samples=100,
+        n_features=n_features,
+        mean=0.0,
+        cov=Gaussian(variance=1, length_scale=0.1),
+        random_state=0,
+    )
+
+
+def test_std_gaussian_fourier(
+    start: int,
+    stop: int,
+    n_features: int,
+    n_basis: int,
+    gaussian_process: FDataGrid,
+) -> None:
+    """Test standard deviation: Gaussian processes and a Fourier basis."""
+    fourier_basis = FourierBasis(n_basis=n_basis, domain_range=(0, 1))
+    fd = gaussian_process.to_basis(fourier_basis)
+
+    std_fd = std(fd)
+    grid = np.linspace(start, stop, n_features)
+    almost_std_fd = std(fd.to_grid(grid)).to_basis(fourier_basis)
+
+    inner_grid_limit = n_features // 10
+    inner_grid = grid[inner_grid_limit:-inner_grid_limit]
+    np.testing.assert_allclose(
+        std_fd(inner_grid),
+        almost_std_fd(inner_grid),
+        rtol=1e-3,
+    )
+
+    outer_grid = grid[:inner_grid_limit] + grid[-inner_grid_limit:]
+    np.testing.assert_allclose(
+        std_fd(outer_grid),
+        almost_std_fd(outer_grid),
+        rtol=1e-2,
+    )