Statistics

AbstractEstimationMethod

Abstract type for defining statistical estimation methods.

CyclicProximalPointEstimation <: AbstractEstimationMethod

Method for estimation using the cyclic proximal point technique.

ExtrinsicEstimation <: AbstractEstimationMethod

Method for estimation in the ambient space and projecting to the manifold.

For mean estimation, GeodesicInterpolation is used for mean estimation in the ambient space.

GeodesicInterpolation <: AbstractEstimationMethod

Repeated weighted geodesic interpolation method for estimating the Riemannian center of mass.

The algorithm proceeds with the following simple online update:

\[\begin{aligned} μ_1 &= x_1\\ t_k &= \frac{w_k}{\sum_{i=1}^k w_i}\\ μ_{k} &= γ_{μ_{k-1}}(x_k; t_k), \end{aligned}\]

where $x_k$ are points, $w_k$ are weights, $μ_k$ is the $k$th estimate of the mean, and $γ_x(y; t)$ is the point at time $t$ along the shortest_geodesic between points $x,y ∈ \mathcal M$. The algorithm terminates when all $x_k$ have been considered. In the Euclidean case, this exactly computes the weighted mean.

The algorithm has been shown to converge asymptotically with the sample size for the following manifolds equipped with their default metrics when all sampled points are in an open geodesic ball about the mean with corresponding radius (see GeodesicInterpolationWithinRadius):

• All simply connected complete Riemannian manifolds with non-positive sectional curvature at radius $∞$ ^[Cheng2016], in particular:
  • Euclidean
  • SymmetricPositiveDefinite ^[Ho2013]
• Other manifolds:
  • Sphere: $\frac{π}{2}$ ^{[Salehian2015]}
  • Grassmann: $\frac{π}{4}$ ^{[Chakraborty2015]}
  • Stiefel/Rotations: $\frac{π}{2 \sqrt 2}$ ^{[Chakraborty2019]}

For online variance computation, the algorithm additionally uses an analogous recursion to the weighted Welford algorithm ^[West1979].

GeodesicInterpolationWithinRadius{T} <: AbstractEstimationMethod

Estimation of Riemannian center of mass using GeodesicInterpolation with fallback to GradientDescentEstimation if any points are outside of a geodesic ball of specified radius around the mean.

Constructor

GeodesicInterpolationWithinRadius(radius)

GradientDescentEstimation <: AbstractEstimationMethod

Method for estimation using gradient descent.

WeiszfeldEstimation <: AbstractEstimationMethod

Method for estimation using the Weiszfeld algorithm for the median

default_estimation_method(M::AbstractManifold, f)

Specify a default AbstractEstimationMethod for an AbstractManifold for a function f, e.g. the median or the mean.

Note that his function is decorated, so it can inherit from the embedding, for example for the IsEmbeddedSubmanifold trait.

Statistics.cov(
    M::AbstractManifold,
    x::AbstractVector;
    basis::AbstractBasis=DefaultOrthonormalBasis(),
    tangent_space_covariance_estimator::CovarianceEstimator=SimpleCovariance(;
        corrected=true,
    ),
    mean_estimation_method::AbstractEstimationMethod=GradientDescentEstimation(),
    inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method(
        M,
    ),
)

Estimate the covariance matrix of a set of points x on manifold M. Since the covariance matrix on a manifold is a rank 2 tensor, the function returns its coefficients in basis induced by the given tangent space basis. See Section 5 of ^[Pennec2006] for details.

The mean is calculated using the specified mean_estimation_method using [mean](@ref Statistics.mean(::AbstractManifold, ::AbstractVector, ::AbstractEstimationMethod), and tangent vectors at this mean are calculated using the provided inverse_retraction_method. Finally, the covariance matrix in the tangent plane is estimated using the Euclidean space estimator tangent_space_covariance_estimator. The type CovarianceEstimator is defined in StatsBase.jl and examples of covariance estimation methods can be found in CovarianceEstimation.jl.

mean!(M::AbstractManifold, y, x::AbstractVector[, w::AbstractWeights]; kwargs...)
mean!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::AbstractEstimationMethod;
    kwargs...,
)

Compute the mean in-place in y.

mean(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::ExtrinsicEstimation;
    kwargs...,
)

Estimate the Riemannian center of mass of x using ExtrinsicEstimation, i.e. by computing the mean in the embedding and projecting the result back. You can specify an extrinsic_method to specify which mean estimation method to use in the embedding, which defaults to GeodesicInterpolation.

See mean for a description of the remaining kwargs.

mean(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::GeodesicInterpolationWithinRadius;
    kwargs...,
)

Estimate the Riemannian center of mass of x using GeodesicInterpolationWithinRadius.

See mean for a description of kwargs.

mean(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::GeodesicInterpolation;
    shuffle_rng=nothing,
    retraction::AbstractRetractionMethod = default_retraction_method(M),
    inverse_retraction::AbstractInverseRetractionMethod = default_inverse_retraction_method(M),
    kwargs...,
)

Estimate the Riemannian center of mass of x in an online fashion using repeated weighted geodesic interpolation. See GeodesicInterpolation for details.

If shuffle_rng is provided, it is used to shuffle the order in which the points are considered for computing the mean.

Optionally, pass retraction and inverse_retraction method types to specify the (inverse) retraction.

mean(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...)

Compute the (optionally weighted) Riemannian center of mass also known as Karcher mean of the vector x of points on the AbstractManifold M, defined as the point that satisfies the minimizer

\[\argmin_{y ∈ \mathcal M} \frac{1}{2 \sum_{i=1}^n w_i} \sum_{i=1}^n w_i\mathrm{d}_{\mathcal M}^2(y,x_i),\]

where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian distance.

In the general case, the GradientDescentEstimation is used to compute the mean. mean( M::AbstractManifold, x::AbstractVector, [w::AbstractWeights,] method::AbstractEstimationMethod=defaultestimationmethod(M); kwargs..., )

Compute the mean using the specified method.

mean(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::GradientDescentEstimation;
    p0=x[1],
    stop_iter=100,
    retraction::AbstractRetractionMethod = default_retraction_method(M),
    inverse_retraction::AbstractInverseRetractionMethod = default_retraction_method(M),
    kwargs...,
)

Compute the mean using the gradient descent scheme GradientDescentEstimation.

Optionally, provide p0, the starting point (by default set to the first data point). stop_iter denotes the maximal number of iterations to perform and the kwargs... are passed to isapprox to stop, when the minimal change between two iterates is small. For more stopping criteria check the Manopt.jl package and use a solver therefrom.

Optionally, pass retraction and inverse_retraction method types to specify the (inverse) retraction.

The Theory stems from^{[Karcher1977]} and is also described in^{[PennecArsigny2013]} as the exponential barycenter. The algorithm is further described in^[Afsari2013].

median!(M::AbstractManifold, y, x::AbstractVector[, w::AbstractWeights]; kwargs...)
median!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::AbstractEstimationMethod;
    kwargs...,
)

computes the median in-place in y.

median(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::CyclicProximalPointEstimation;
    p0=x[1],
    stop_iter=1000000,
    retraction::AbstractRetractionMethod = default_retraction_method(M),
    inverse_retraction::AbstractInverseRetractionMethod = default_inverse_retraction_method(M),
    kwargs...,
)

Compute the median using CyclicProximalPointEstimation.

Optionally, provide p0, the starting point (by default set to the first data point). stop_iter denotes the maximal number of iterations to perform and the kwargs... are passed to isapprox to stop, when the minimal change between two iterates is small. For more stopping criteria check the Manopt.jl package and use a solver therefrom.

Optionally, pass retraction and inverse_retraction method types to specify the (inverse) retraction.

The algorithm is further described in ^{[Bačák2014]}.

median(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::ExtrinsicEstimation;
    extrinsic_method = CyclicProximalPointEstimation(),
    kwargs...,
)

Estimate the median of x using ExtrinsicEstimation, i.e. by computing the median in the embedding and projecting the result back. You can specify an extrinsic_method to specify which median estimation method to use in the embedding, which defaults to CyclicProximalPointEstimation.

See median for a description of kwargs.

median(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::WeiszfeldEstimation;
    α = 1.0,
    p0=x[1],
    stop_iter=2000,
    retraction::AbstractRetractionMethod = default_retraction_method(M),
    inverse_retraction::AbstractInverseRetractionMethod = default_inverse_retraction_method(M),
    kwargs...,
)

Compute the median using WeiszfeldEstimation.

Optionally, provide p0, the starting point (by default set to the first data point). stop_iter denotes the maximal number of iterations to perform and the kwargs... are passed to isapprox to stop, when the minimal change between two iterates is small. For more stopping criteria check the Manopt.jl package and use a solver therefrom.

The parameter $α\in (0,2]$ is a step size.

The algorithm is further described in ^{[FletcherVenkatasubramanianJoshi2008]}, especially the update rule in Eq. (6), i.e. Let $q_{k}$ denote the current iterate, $n$ the number of points $x_1,\ldots,x_n$, and

\[I_k = \bigl\{ i \in \{1,\ldots,n\} \big| x_i \neq q_k \bigr\}\]

all indices of points that are not equal to the current iterate. Then the update reads $q_{k+1} = \exp_{q_k}(αX)$, where

\[X = \frac{1}{s}\sum_{i\in I_k} \frac{w_i}{d_{\mathcal M}(q_k,x_i)}\log_{q_k}x_i \quad \text{ with } \quad s = \sum_{i\in I_k} \frac{w_i}{d_{\mathcal M}(q_k,x_i)},\]

and where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian distance.

Optionally, pass retraction and inverse_retraction method types to specify the (inverse) retraction, which by default use the exponential and logarithmic map, respectively.

median(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...)
median(
    M::AbstractManifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::AbstractEstimationMethod;
    kwargs...,
)

Compute the (optionally weighted) Riemannian median of the vector x of points on the AbstractManifold M, defined as the point that satisfies the minimizer

\[\argmin_{y ∈ \mathcal M} \frac{1}{\sum_{i=1}^n w_i} \sum_{i=1}^n w_i\mathrm{d}_{\mathcal M}(y,x_i),\]

where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian distance. This function is nonsmooth (i.e nondifferentiable).

In the general case, the CyclicProximalPointEstimation is used to compute the median. However, this default may be overloaded for specific manifolds.

Compute the median using the specified method.

std(M, x, m=mean(M, x); corrected=true, kwargs...)
std(M, x, w::AbstractWeights, m=mean(M, x, w); corrected=false, kwargs...)

compute the optionally weighted standard deviation of a Vector x of n data points on the AbstractManifold M, i.e.

\[\sqrt{\frac{1}{c} \sum_{i=1}^n w_i d_{\mathcal M}^2 (x_i,m)},\]

where c is a correction term, see Statistics.std. The mean of x can be specified as m, and the corrected variance can be activated by setting corrected=true.

var(M, x, m=mean(M, x); corrected=true)
var(M, x, w::AbstractWeights, m=mean(M, x, w); corrected=false)

compute the (optionally weighted) variance of a Vector x of n data points on the AbstractManifold M, i.e.

\[\frac{1}{c} \sum_{i=1}^n w_i d_{\mathcal M}^2 (x_i,m),\]

where c is a correction term, see Statistics.var. The mean of x can be specified as m, and the corrected variance can be activated by setting corrected=true. All further kwargs... are passed to the computation of the mean (if that is not provided).

kurtosis(M::AbstractManifold, x::AbstractVector, k::Int[, w::AbstractWeights], m=mean(M, x[, w]))

Compute the excess kurtosis of points in x on manifold M. Optionally provide weights w and/or a precomputed mean m.

mean_and_std(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...) -> (mean, std)

Compute the mean and the standard deviation std simultaneously.

mean_and_std(
    M::AbstractManifold,
    x::AbstractVector
    [w::AbstractWeights,]
    method::AbstractEstimationMethod;
    kwargs...,
) -> (mean, var)

Use the method for simultaneously computing the mean and standard deviation. To use a mean-specific method, call mean and then std.

mean_and_var(
    M::AbstractManifold,
    x::AbstractVector
    [w::AbstractWeights,]
    method::GeodesicInterpolationWithinRadius;
    kwargs...,
) -> (mean, var)

Use repeated weighted geodesic interpolation to estimate the mean. Simultaneously, use a Welford-like recursion to estimate the variance.

See GeodesicInterpolationWithinRadius and mean_and_var for more information.

mean_and_var(
    M::AbstractManifold,
    x::AbstractVector
    [w::AbstractWeights,]
    method::GeodesicInterpolation;
    shuffle_rng::Union{AbstractRNG,Nothing} = nothing,
    retraction::AbstractRetractionMethod = default_retraction_method(M),
    inverse_retraction::AbstractInverseRetractionMethod = default_inverse_retraction_method(M),
    kwargs...,
) -> (mean, var)

Use the repeated weighted geodesic interpolation to estimate the mean. Simultaneously, use a Welford-like recursion to estimate the variance.

If shuffle_rng is provided, it is used to shuffle the order in which the points are considered. Optionally, pass retraction and inverse_retraction method types to specify the (inverse) retraction.

See GeodesicInterpolation for details on the geodesic interpolation method.

mean_and_var(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...) -> (mean, var)

Compute the mean and the variance simultaneously. See those functions for a description of the arguments.

mean_and_var(
    M::AbstractManifold,
    x::AbstractVector
    [w::AbstractWeights,]
    method::AbstractEstimationMethod;
    kwargs...,
) -> (mean, var)

Use the method for simultaneously computing the mean and variance. To use a mean-specific method, call mean and then var.

moment(M::AbstractManifold, x::AbstractVector, k::Int[, w::AbstractWeights], m=mean(M, x[, w]))

Compute the kth central moment of points in x on manifold M. Optionally provide weights w and/or a precomputed mean.

skewness(M::AbstractManifold, x::AbstractVector, k::Int[, w::AbstractWeights], m=mean(M, x[, w]))

Compute the standardized skewness of points in x on manifold M. Optionally provide weights w and/or a precomputed mean m.