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.

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

Compute the mean in-place in y.

mean(
    M::Manifold,
    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::Manifold,
    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::Manifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::GeodesicInterpolation;
    shuffle_rng=nothing,
    retraction::AbstractRetractionMethod = ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
    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::Manifold, 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 ManifoldM, 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::Manifold, x::AbstractVector, [w::AbstractWeights,] method::AbstractEstimationMethod; kwargs..., )

Compute the mean using the specified method.

mean(
    M::Manifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::GradientDescentEstimation;
    p0=x[1],
    stop_iter=100,
    retraction::AbstractRetractionMethod = ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
    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::Manifold, y, x::AbstractVector[, w::AbstractWeights]; kwargs...)
median!(
    M::Manifold,
    y,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::AbstractEstimationMethod;
    kwargs...,
)

computes the median in-place in y.

median(
    M::Manifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::CyclicProximalPointEstimation;
    p0=x[1],
    stop_iter=1000000,
    retraction::AbstractRetractionMethod = ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
    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::Manifold,
    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::Manifold, x::AbstractVector[, w::AbstractWeights]; kwargs...)
median(
    M::Manifold,
    x::AbstractVector,
    [w::AbstractWeights,]
    method::AbstractEstimationMethod;
    kwargs...,
)

Compute the (optionally weighted) Riemannian median of the vector x of points on the ManifoldM, 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 Vectorx of n data points on the ManifoldM, 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 Vectorx of n data points on the ManifoldM, 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::Manifold, 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 meanm.

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

Compute the mean and the standard deviation std simultaneously.

mean_and_std(
    M::Manifold,
    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::Manifold,
    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::Manifold,
    x::AbstractVector
    [w::AbstractWeights,]
    method::GeodesicInterpolation;
    shuffle_rng::Union{AbstractRNG,Nothing} = nothing,
    retraction::AbstractRetractionMethod = ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
    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::Manifold, 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::Manifold,
    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::Manifold, 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::Manifold, 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 meanm.