spnorm: Spherical Normal Distribution

Description Usage Arguments Value Examples

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

1
2
3
4
5
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.i

...

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.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
# -------------------------------------------------------------------
#          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 June 20, 2021, 5:07 p.m.