View source: R/all_vmcos_fns.R
rvmcos | R Documentation |
The bivariate von Mises cosine model
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, ... )
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. |
The bivariate von Mises cosine model density at the point x = (x_1, x_2) is given by
f(x) = C_c (κ_1, κ_2, κ_3) \exp(κ_1 \cos(T_1) + κ_2 \cos(T_2) + κ_3 \cos(T_1 - T_2))
where
T_1 = x_1 - μ_1; T_2 = x_2 - μ_2
and C_c (κ_1, κ_2, κ_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.
dvmcos
gives the density and rvmcos
generates random deviates.
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.