# macg: Matrix Angular Central Gaussian Distribution In Riemann: Learning with Data on Riemannian Manifolds

## 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

 1 2 3 4 5 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 maxitermaximum number of iterations to be run (default:50). epstolerance 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.

acg
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 # ------------------------------------------------------------------- # 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)