| rvmcos {BAMBI} | R Documentation |
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,
...
)
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; |
mu1, mu2 |
vectors of mean parameters. |
method |
Rejection sampling method to be used. Available choices are |
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.
Value
dvmcos gives the density and rvmcos generates random deviates.
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)
## 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)