spnorm: Spherical Normal Distribution

In Riemann: Learning with Data on Riemannian Manifolds

Spherical Normal Distribution

Description

We provide tools for an isotropic spherical normal (SN) distributions on a (p-1)-sphere in \mathbf{R}^p for sampling, density evaluation, and maximum likelihood estimation of the parameters where the density is defined as

f_{SN}(x; μ, λ) = \frac{1}{Z(λ)} \exp ≤ft( -\frac{λ}{2} d^2(x,μ) \right)

for location and concentration parameters μ and λ respectively and the normalizing constant Z(λ).

Usage

dspnorm(data, mu, lambda, log = FALSE)

rspnorm(n, mu, lambda)

mle.spnorm(data, method = c("Newton", "Halley", "Optimize", "DE"), ...)

Arguments

data	data vectors in form of either an (n\times p) matrix or a length-n list. See wrap.sphere for descriptions on supported input types.
mu	a length-p unit-norm vector of location.
lambda	a concentration parameter that is positive.
log	a logical; TRUE to return log-density, FALSE for densities without logarithm applied.
n	the number of samples to be generated.
method	an algorithm name for concentration parameter estimation. It should be one of "Newton","Halley","Optimize", and "DE" (case sensitive).
...	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

dspnorm gives a vector of evaluated densities given samples. rspnorm generates unit-norm vectors in \mathbf{R}^p wrapped in a list. mle.spnorm computes MLEs and returns a list containing estimates of location (mu) and concentration (lambda) parameters.

References

\insertRef

hauberg_2018_DirectionalStatisticsSphericalRiemann

\insertRef

you_2022_ParameterEstimationModelbasedRiemann

Examples

# -------------------------------------------------------------------
#          Example with Spherical Normal Distribution
#
# Given a fixed set of parameters, generate samples and acquire MLEs.
# Especially, we will see the evolution of estimation accuracy.
# -------------------------------------------------------------------
## DEFAULT PARAMETERS
true.mu  = c(1,0,0,0,0)
true.lbd = 5

## GENERATE DATA N=1000
big.data = rspnorm(1000, true.mu, true.lbd)

## ITERATE FROM 50 TO 1000 by 10
idseq = seq(from=50, to=1000, by=10)
nseq  = length(idseq)

hist.mu  = rep(0, nseq)
hist.lbd = rep(0, nseq)

for (i in 1:nseq){
  small.data = big.data[1:idseq[i]]          # data subsetting
  small.MLE  = mle.spnorm(small.data) # compute MLE
  
  hist.mu[i]  = acos(sum(small.MLE$mu*true.mu)) # difference in mu
  hist.lbd[i] = small.MLE$lambda
}

## VISUALIZE
opar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
plot(idseq, hist.mu,  "b", pch=19, cex=0.5, main="difference in location")
plot(idseq, hist.lbd, "b", pch=19, cex=0.5, main="concentration param")
abline(h=true.lbd, lwd=2, col="red")
par(opar)

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