scoreMatchWnBvm: Score and moment matching of a univariate or bivariate...

View source: R/bessel.R

scoreMatchWnBvmR Documentation

Score and moment matching of a univariate or bivariate wrapped normal by a von Mises

Description

Given a wrapped normal density, find the parameters of a von Mises that matches it according to two characteristics: moments and scores. Score matching estimators are available for univariate and bivariate cases and moment matching only for the former.

Usage

scoreMatchWnBvm(Sigma = NULL, invSigma)

scoreMatchWnVm(sigma, sigma2 = NULL)

momentMatchWnVm(sigma, sigma2 = NULL)

Arguments

Sigma, invSigma

covariance or precision matrix of the bivariate wrapped normal.

sigma, sigma2

standard deviation or variance of the wrapped normal.

Details

If the precision matrix is singular or if there are no solutions for the score matching estimator, c(0, 0, 0) is returned.

Value

Vector of parameters (\kappa_1,\kappa_2,\lambda), where (\kappa_1,\kappa_2,2\lambda) is a suitable input for kappa in dBvm.

References

Mardia, K. V., Kent, J. T., and Laha, A. K. (2016). Score matching estimators for directional distributions. arXiv:1604.0847. https://arxiv.org/abs/1604.08470

Examples

# Univariate WN approximation
sigma <- 0.5
curve(dWn1D(x = x, mu = 0, sigma = sigma), from = -pi, to = pi,
      ylab = "Density", ylim = c(0, 1))
curve(dVm(x = x, mu = 0, kappa = momentMatchWnVm(sigma = sigma)), from = -pi,
      to = pi, col = "red", add = TRUE)
curve(dVm(x = x, mu = 0, kappa = scoreMatchWnVm(sigma = sigma)), from = -pi,
      to = pi, col = "green", add = TRUE)

# Bivariate WN approximation

# WN
alpha <- c(2, 1, 1)
sigma <- c(1, 1)
mu <- c(pi / 2, pi / 2)
x <- seq(-pi, pi, l = 101)[-101]
plotSurface2D(x, x, f = function(x) dStatWn2D(x = x, alpha = alpha, mu = mu,
                                              sigma = sigma), fVect = TRUE)
A <- alphaToA(alpha = alpha, sigma = sigma)
S <- 0.5 * solve(A) %*% diag(sigma)

# Score matching
kappa <- scoreMatchWnBvm(Sigma = S)

# dBvm uses lambda / 2 in the exponent
plotSurface2D(x, x, f = function(x) dBvm(x = x, mu = mu,
                                        kappa = c(kappa[1:2], 2 * kappa[3])),
             fVect = TRUE)

# With singular Sigma
invSigma <- matrix(c(1, sqrt(0.999), sqrt(0.999), 1), nrow = 2, ncol = 2)
scoreMatchWnBvm(invSigma = invSigma)
invSigma <- matrix(1, nrow = 2, ncol = 2)
scoreMatchWnBvm(invSigma = invSigma)

egarpor/sdetorus documentation built on March 4, 2024, 1:23 a.m.