macg: Matrix Angular Central Gaussian Distribution

In Riemann: Learning with Data on Riemannian Manifolds

Matrix Angular Central Gaussian Distribution

Description

For Stiefel and Grassmann manifolds St(r,p) and Gr(r,p), the matrix variant of ACG distribution is known as Matrix Angular Central Gaussian (MACG) distribution MACG_{p,r}(Σ) with density

f(X\vert Σ) = |Σ|^{-r/2} |X^\top Σ^{-1} X|^{-p/2}

where Σ is a (p\times p) symmetric positive-definite matrix. Similar to vector-variate ACG case, we follow a convention that tr(Σ)=p.

Usage

dmacg(datalist, Sigma)

rmacg(n, r, Sigma)

mle.macg(datalist, ...)

Arguments

datalist	a list of (p\times r) orthonormal matrices.
Sigma	a (p\times p) symmetric positive-definite matrix.
n	the number of samples to be generated.
r	the number of basis.
...	extra parameters for computations, including maxiter maximum number of iterations to be run (default:50). eps tolerance level for stopping criterion (default: 1e-5).

Value

dmacg gives a vector of evaluated densities given samples. rmacg generates (p\times r) orthonormal matrices wrapped in a list. mle.macg estimates the SPD matrix Σ.

References

\insertRef

chikuse_matrix_1990Riemann

\insertRef

mardia_directional_1999Riemann

Kent JT, Ganeiber AM, Mardia KV (2013). "A new method to simulate the Bingham and related distributions in directional data analysis with applications." arXiv:1310.8110.

See Also

acg

Examples

# -------------------------------------------------------------------
#          Example with Matrix Angular Central Gaussian Distribution
#
# Given a fixed Sigma, generate samples and estimate Sigma via ML.
# -------------------------------------------------------------------
## GENERATE AND MLE in St(2,5)/Gr(2,5)
#  Generate data
Strue = diag(5)                  # true SPD matrix
sam1  = rmacg(n=50,  r=2, Strue) # random samples
sam2  = rmacg(n=100, r=2, Strue) # random samples

#  MLE
Smle1 = mle.macg(sam1)
Smle2 = mle.macg(sam2)

#  Visualize
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,3), pty="s")
image(Strue[,5:1], axes=FALSE, main="true SPD")
image(Smle1[,5:1], axes=FALSE, main="MLE with n=50")
image(Smle2[,5:1], axes=FALSE, main="MLE with n=100")
par(opar)

Riemann documentation built on March 18, 2022, 7:55 p.m.