BAMBI (version 2.3.5)

rvmcos: The bivariate von Mises cosine model

Description

The bivariate von Mises cosine model

Usage

rvmcos(
  n,
  kappa1 = 1,
  kappa2 = 1,
  kappa3 = 0,
  mu1 = 0,
  mu2 = 0,
  method = "naive"
)
dvmcos(
  x,
  kappa1 = 1,
  kappa2 = 1,
  kappa3 = 0,
  mu1 = 0,
  mu2 = 0,
  log = FALSE,
  ...
)

Value

dvmcos gives the density and rvmcos generates random deviates.

Arguments

n

number of observations. Ignored if at least one of the other parameters have length k > 1, in which case, all the parameters are recycled to length k to produce k random variates.

kappa1, kappa2, kappa3

vectors of concentration parameters; kappa1, kappa2 > 0.

mu1, mu2

vectors of mean parameters.

method

Rejection sampling method to be used. Available choices are "naive" (default) or "vmprop". See details.

x

bivariate vector or a two-column matrix with each row being a bivariate vector of angles (in radians) where the densities are to be evaluated.

log

logical. Should the log density be returned instead?

...

additional arguments to be passed to dvmcos. See details.

Details

The bivariate von Mises cosine model density at the point \(x = (x_1, x_2)\) is given by $$f(x) = C_c (\kappa_1, \kappa_2, \kappa_3) \exp(\kappa_1 \cos(T_1) + \kappa_2 \cos(T_2) + \kappa_3 \cos(T_1 - T_2))$$ where $$T_1 = x_1 - \mu_1; T_2 = x_2 - \mu_2$$ and \(C_c (\kappa_1, \kappa_2, \kappa_3)\) denotes the normalizing constant for the cosine model.

Because \(C_c\) involves an infinite alternating series with product of Bessel functions, if kappa3 < -5 or max(kappa1, kappa2, abs(kappa3)) > 50, \(C_c\) is evaluated numerically via (quasi) Monte carlo method for numerical stability. These (quasi) random numbers can be provided through the argument qrnd, which must be a two column matrix, with each element being a (quasi) random number between 0 and 1. Alternatively, if n_qrnd is provided (and qrnd is missing), a two dimensional sobol sequence of size n_qrnd is generated via the function sobol from the R package qrng. If none of qrnd or n_qrnd is available, a two dimensional sobol sequence of size 1e4 is used. By default Monte Carlo approximation is used only if kappa3 < -5 or max(kappa1, kappa2, abs(kappa3)) > 50. However, a forced Monte Carlo approximation can be made (irrespective of the choice of kappa1, kappa2 and kappa3) by setting force_approx_const = TRUE. See examples.

Examples

kappa1 <- c(1, 2, 3)
kappa2 <- c(1, 6, 5)
kappa3 <- c(0, 1, 2)
mu1 <- c(1, 2, 5)
mu2 <- c(0, 1, 3)
x <- diag(2, 2)
n <- 10

# when x is a bivariate vector and parameters are all scalars,
# dvmcos returns single density
dvmcos(x[1, ], kappa1[1], kappa2[1], kappa3[1], mu1[1], mu2[1])

# when x is a two column matrix and parameters are all scalars,
# dmvsin returns a vector of densities calculated at the rows of
# x with the same parameters
dvmcos(x, kappa1[1], kappa2[1], kappa3[1], mu1[1], mu2[1])

# if x is a bivariate vector and at least one of the parameters is
# a vector, all parameters are recycled to the same length, and
# dvmcos returns a vector with ith element being the density
# evaluated at x with parameter values kappa1[i], kappa2[i],
# kappa3[i], mu1[i] and mu2[i]
dvmcos(x[1, ], kappa1, kappa2, kappa3, mu1, mu2)

# if x is a two column matrix and at least one of the parameters is
# a vector, rows of x and the parameters are recycled to the same
# length, and dvmcos returns a vector with ith element being the
# density evaluated at ith row of x with parameter values kappa1[i],
# kappa2[i], # kappa3[i], mu1[i] and mu2[i]
dvmcos(x, kappa1, kappa2, kappa3, mu1, mu2)

# when parameters are all scalars, number of observations generated
# by rvmcos is n
rvmcos(n, kappa1[1], kappa2[1], kappa3[1], mu1[1], mu2[1])

# when at least one of the parameters is a vector, all parameters are
# recycled to the same length, n is ignored, and the number of
# observations generated by rvmcos is the same as the length of the
# recycled vectors
rvmcos(n, kappa1, kappa2, kappa3, mu1, mu2)


# \donttest{
## Visualizing (quasi) Monte Carlo based approximations of
## the normalizing constant through density evaluations.

# "good" setup, where the analytic formula for C_c can be
# calculated without numerical issues
# kappa1 = 1, kappa2 = 1, kappa3 = -2, mu1 = pi, mu2 = pi

n_qrnd <-  (1:500)*20
# analytic
good.a <- dvmcos(c(3,3), 1, 1, -2, pi, pi, log=TRUE)
# using quasi Monte Carlo
good.q <- sapply(n_qrnd,
                 function(j)
                   dvmcos(c(3,3), 1, 1, -2, pi, pi,
                          log=TRUE, n_qrnd = j,
                          force_approx_const = TRUE))
# using ordinary Monte Carlo
set.seed(1)
good.r <- sapply(n_qrnd,
                 function(j)
                   dvmcos(c(3,3), 1, 1, -2, pi, pi,
                          log=TRUE,
                          qrnd = matrix(runif(2*j), ncol = 2),
                          force_approx_const = TRUE))


plot(n_qrnd, good.q, ylim = range(good.a, good.q, good.r),
     col = "orange", type = "l",
     ylab = "",
     main = "dvmcos(c(3,3), 1, 1, -2, pi, pi, log = TRUE)")
points(n_qrnd, good.r, col = "skyblue", type = "l")
abline(h = good.a, lty = 2, col = "grey")
legend("topright",
       legend = c("Sobol", "Random", "Analytic"),
       col = c("orange", "skyblue", "grey"),
       lty = c(1, 1, 2))


# "bad" setup, where the calculating C_c
# numerically using the analytic formula is problematic
# kappa1 = 100, kappa2 = 100, kappa3 = -200, mu1 = pi, mu2 = pi

n_qrnd <-  (1:500)*20

# using quasi Monte Carlo
bad.q <- sapply(n_qrnd,
                function(j)
                  dvmcos(c(3,3), 100, 100, -200, pi, pi,
                         log=TRUE, n_qrnd = j,
                         force_approx_const = TRUE))
# using ordinary Monte Carlo
set.seed(1)
bad.r <- sapply(n_qrnd,
                function(j)
                  dvmcos(c(3,3), 100, 100, -200, pi, pi,
                         log=TRUE,
                         qrnd = matrix(runif(2*j), ncol = 2),
                         force_approx_const = TRUE))


plot(n_qrnd, bad.q, ylim = range(bad.q, bad.r),
     col = "orange", type = "l",
     ylab = "",
     main = "dvmcos(c(3,3), 100, 100, -200, pi, pi, log = TRUE)")
points(n_qrnd, bad.r, col = "skyblue", type = "l")
legend("topright",
       legend = c("Sobol", "Random"),
       col = c("orange", "skyblue"), lty = 1)
# }