Nothing
R/all_vmsin_fns.R
In BAMBI: Bivariate Angular Mixture Models
Defines functions
rvmsin_1par
vmsin_var_cor_singlepar
vmsin_var_cor_singlepar_numeric
fit_vmsinmix
dvmsinmix
rvmsinmix
dvmsin
rvmsin
Documented in
dvmsin
dvmsinmix
fit_vmsinmix
rvmsin
rvmsinmix
#' The bivariate von Mises sine model
#'
#' @param 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.
#' @param 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.
#' @param mu1,mu2 vectors of mean parameters.
#' @param kappa1,kappa2,kappa3 vectors of concentration parameters; \code{kappa1, kappa2 > 0}.
#' @param log logical. Should the log density be returned instead?
#' @param method Rejection sampling method to be used. Available choices are \code{"naive"} (default) or \code{"vmprop"}. See details.
#' @details
#' The bivariate von Mises sine model density at the point \eqn{x = (x_1, x_2)} is given by
#' \deqn{f(x) = C_s (\kappa_1, \kappa_2, \kappa_3) \exp(\kappa_1 \cos(T_1) + \kappa_2 \cos(T_2) + \kappa_3 \sin(T_1) \sin(T_2))}
#' where
#' \deqn{T_1 = x_1 - \mu_1; T_2 = x_2 - \mu_2}
#' and \eqn{C_s (\kappa_1, \kappa_2, \kappa_3)} denotes the normalizing constant for the sine model.
#'
#' Two different rejection sampling methods are implemented for random generation. If \code{method = "vmprop"}, then first the y-marginal
#' is drawn from the associated marginal density, and then x is generated from the conditional distributio of x given y. The marginal generation of
#' y is implemented in a rejection sampling scheme with proposal being either von Mises (if the target marginal density is unimodal), or a mixture of
#' von Mises (if bimodal), with optimally chosen concentration. This the method suggested in Mardia et al. (2007). On the other hand, when
#' \code{method = "naive"} (default) a (naive) bivariate rejection sampling scheme with (bivariate) uniform propsoal is used.
#'
#' Note that although method = \code{"vmprop"} may provide better efficiency when the density is highly concentrated, it does have
#' an (often substantial) overhead due to the optimziation step required to find a reasonable proposal concentration parameter.
#' This can compensate the efficiency gains of this method, especially when \code{n} is not large.
#'
#'
#'
#'
#' @return \code{dvmsin} gives the density and \code{rvmsin} 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,
#' # dvmsin returns single density
#' dvmsin(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
#' dvmsin(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
#' # dvmsin returns a vector of with ith element being the density
#' # evaluated at x with parameter values kappa1[i], kappa2[i],
#' # kappa3[i], mu1[i] and mu2[i]
#' dvmsin(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 dvmsin returns a vector of 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]
#' dvmsin(x[1, ], kappa1, kappa2, kappa3, mu1, mu2)
#'
#' # when parameters are all scalars, number of observations generated
#' # by rvmsin is n
#' rvmsin(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 rvmsin is the same as the length of the
#' # recycled vectors
#' rvmsin(n, kappa1, kappa2, kappa3, mu1, mu2)
#'
#' @export
rvmsin <- function(n, kappa1=1, kappa2=1,
kappa3=0, mu1=0, mu2=0, method="naive")
{
if(any(c(kappa1, kappa2) < 0)) stop("kappa1 and kappa2 must be nonnegative")
if(any(mu1 < 0 | mu1 >= 2*pi)) mu1 <- prncp_reg(mu1)
if(any(mu2 < 0 | mu2 >= 2*pi)) mu2 <- prncp_reg(mu2)
if(n < 0) stop("invalid n")
# if (is.null(method)) {
# if (n > 1e5) method <- "vmprop"
# else method <- "naive"
# }
if (!method %in% c("naive", "vmprop"))
stop("method must be either \'naive\' or \'vmprop\'")
if(max(length(kappa1), length(kappa2), length(kappa3),
length(mu1), length(mu2)) > 1) {
expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
k1 <- expanded[[1]]; k2 <- expanded[[2]]; k3 <- expanded[[3]]
mu1 <- expanded[[4]]; mu2 <- expanded[[5]]
t(vapply(1:length(k1),
function(j) rvmsin_1par(1, k1[j], k2[j], k3[j],
mu1[j], mu2[j], "naive"), c(0, 0)))
} else {
# if (is.null(method) & max(kappa1, kappa2, abs(kappa3)) < 0.1)
# method <- "vmprop"
rvmsin_1par(n, kappa1, kappa2, kappa3, mu1, mu2, method)
}
}
# rvmsin_naive <- function(n, kappa1=1, kappa2=1,
# kappa3=0, mu1=0, mu2=0) {
#
# opt_obj <- function(k1=1, k2=1, k3=0, mu1=0, mu2=0) {
# # for numerical stability, if k3 < 0, fabs(k1+k3) < 1e-5 or fabs(k2+k3) < 1e-5
# # make k3 = k3 * (1+1e-5)
# if (k3 < 0) {
# while (abs(k1 + k3) < 1e-5 || abs(k2 + k3) < 1e-5) {
# k3 = k3*(1+1e-5)
# }
# }
# obj <- optim(c(0,0), fn = function(x) -(k1*cos(x[1]-mu1)+k2*cos(x[2]-mu2)+k3*sin(x[1]-mu1)*sin(x[2]-mu2)),
# gr = function(x) -c(-k1*sin(x[1]-mu1)+k3*cos(x[1]-mu1)*sin(x[2]-mu2),
# -k2*sin(x[2]-mu2)+k3*sin(x[1]-mu1)*cos(x[2]-mu2)))
# -obj$value
#
# }
#
# if(max(length(kappa1), length(kappa2), length(kappa3), length(mu1), length(mu2)) > 1) {
# expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
# k1 <- expanded[[1]]; k2 <- expanded[[2]]; k3 <- expanded[[3]]
# mu1 <- expanded[[4]]; mu2 <- expanded[[5]]
# upper_bd_all <- vapply(1:length(k1),
# function(h) opt_obj(k1[h], k2[h], k3[h], mu1[h], mu2[h]),
# 0)
# rsin_manypar(k1, k2, k3, mu1, mu2, upper_bd_all)
# } else {
# upper_bd <- opt_obj(kappa1, kappa2, kappa3, mu1, mu2)
# rsin_onepar(n, kappa1, kappa2, kappa3, mu1, mu2, upper_bd)
# }
# }
#' @rdname rvmsin
#' @export
dvmsin <- function(x, kappa1=1, kappa2=1, kappa3=0, mu1=0, mu2=0, log = FALSE)
{
if(any(c(kappa1, kappa2) < 0)) stop("kappa1 and kappa2 must be non-negative")
if(any(mu1 < 0 | mu1 >= 2*pi)) mu1 <- prncp_reg(mu1)
if(any(mu2 < 0 | mu2 >= 2*pi)) mu2 <- prncp_reg(mu2)
if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
|| (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
if(max(length(kappa1), length(kappa2), length(kappa3), length(mu1), length(mu2)) > 1) {
expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
kappa1 <- expanded[[1]]
kappa2 <- expanded[[2]]
kappa3 <- expanded[[3]]
mu1 <- expanded[[4]]
mu2 <- expanded[[5]]
}
par.mat <- rbind(kappa1, kappa2, kappa3, mu1, mu2)
n_par <- ncol(par.mat)
if (length(x) == 2) x <- matrix(x, nrow = 1)
n_x <- nrow(x)
if (n_par == 1) {
log_den <- c(ldsin_manyx_onepar(x, kappa1, kappa2, kappa3,
mu1, mu2))
} else if (n_x == 1) {
log_den <- c(ldsin_onex_manypar(c(x), kappa1, kappa2, kappa3,
mu1, mu2))
} else {
x_set <- 1:nrow(x)
par_set <- 1:n_par
expndn_set <- expand_args(x_set, par_set)
x_set <- expndn_set[[1]]
par_set <- expndn_set[[2]]
log_den <- c(ldsin_manyx_manypar(x[x_set, ], kappa1[par_set],
kappa2[par_set], kappa3[par_set],
mu1[par_set], mu2[par_set]))
}
if (log) log_den
else exp(log_den)
}
#'
# dvmsin <- function(x, kappa1=1, kappa2=1, kappa3=0, mu1=0, mu2=0)
# {
# if(any(c(kappa1, kappa2) < 0)) stop("kappa1 and kappa2 must be non-negative")
# if(any(mu1 < 0 | mu1 >= 2*pi)) mu1 <- prncp_reg(mu1)
# if(any(mu2 < 0 | mu2 >= 2*pi)) mu2 <- prncp_reg(mu2)
# if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
# || (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
#
# if(max(length(kappa1), length(kappa2), length(kappa3), length(mu1), length(mu2)) > 1) {
# expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
# kappa1 <- expanded[[1]]
# kappa2 <- expanded[[2]]
# kappa3 <- expanded[[3]]
# mu1 <- expanded[[4]]
# mu2 <- expanded[[5]]
# }
#
# par.mat <- rbind(kappa1, kappa2, kappa3, mu1, mu2)
# n_par <- ncol(par.mat)
# if (length(x) == 2) x <- matrix(x, nrow = 1)
# n_x <- nrow(x)
#
#
# if (n_par == 1) {
# log_den <- c(ldsin_manyx_onepar(x, kappa1, kappa2, kappa3,
# mu1, mu2))
#
# } else if (n_x == 1) {
# log_den <- c(ldsin_onex_manypar(c(x), kappa1, kappa2, kappa3,
# mu1, mu2))
# } else {
# x_set <- 1:nrow(x)
# par_set <- 1:n_par
# expndn_set <- expand_args(x_set, par_set)
# x_set <- expndn_set[[1]]
# par_set <- expndn_set[[2]]
# log_den <- c(ldsin_manyx_manypar(x[x_set, ], kappa1[par_set],
# kappa2[par_set], kappa3[par_set],
# mu1[par_set], mu2[par_set]))
# }
#
# if (log) log_den
# else exp(log_den)
# }
#' The bivariate von Mises sine model mixtures
#' @inheritParams rvmsin
#' @param n number of observations.
#' @param x matrix of angles (in radians) where the density is to be evaluated, with each row being a
#' single bivariate vector of angles.
#'
#' @param pmix vector of mixture proportions.
#' @param mu1,mu2 vectors of mean parameters.
#' @param kappa1,kappa2,kappa3 vectors of concentration parameters; \code{kappa1, kappa2 > 0} for each component.
#'
#' @details All the argument vectors \code{pmix, kappa1, kappa2, kappa3, mu1} and \code{mu2} must be of
#' the same length ( = component size of the mixture model), with \eqn{j}-th element corresponding to the
#' \eqn{j}-th component of the mixture distribution.
#' @details The bivariate von Mises sine model mixture distribution with component size \code{K = \link{length}(p.mix)} has density
#' \deqn{g(x) = \sum p[j] * f(x; \kappa_1[j], \kappa_2[j], \kappa_3[j], \mu_1[j], \mu_2[j])}
#' where the sum extends over \eqn{j}; \eqn{p[j]; \kappa_1[j], \kappa_2[j], \kappa_3[j]}; and \eqn{\mu_1[j], \mu_2[j]} respectively denote the mixing proportion,
#' the three concentration parameters and the two mean parameter for the \eqn{j}-th component, \eqn{j = 1, ..., K},
#' and \eqn{f(. ; \kappa_1, \kappa_2, \kappa_3, \mu_1, \mu_2)} denotes the density function of the von Mises sine model
#' with concentration parameters \eqn{\kappa_1, \kappa_2, \kappa_3} and mean parameters \eqn{\mu_1, \mu_2}.
#'
#' @return \code{dvmsinmix} computes the density (vector if x is a two column matrix with more than one row)
#' and \code{rvmsinmix} generates random deviates from the mixture density.
#'
#' @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)
#' pmix <- c(0.3, 0.4, 0.3)
#' x <- diag(2, 2)
#' n <- 10
#'
#' # mixture densities calculated at the rows of x
#' dvmsinmix(x, kappa1, kappa2, kappa3, mu1, mu2, pmix)
#'
#' # number of observations generated from the mixture distribution is n
#' rvmsinmix(n, kappa1, kappa2, kappa3, mu1, mu2, pmix)
#'
#' @export
rvmsinmix <- function(n, kappa1, kappa2, kappa3, mu1, mu2, pmix, method="naive")
{
allpar <- list(kappa1=kappa1, kappa2=kappa2, kappa3=kappa3,
mu1=mu1, mu2=mu2, pmix=pmix)
allpar_len <- listLen(allpar)
if(min(allpar_len) != max(allpar_len))
stop("component size mismatch: number of components of the input parameter vectors differ")
if(any(allpar$pmix < 0)) stop("\'pmix\' must be non-negative")
sum_pmix <- sum(allpar$pmix)
if(signif(sum_pmix, 5) != 1) {
if(sum_pmix <= 0) stop("\'pmix\' must have at least one positive element")
allpar$pmix <- allpar$pmix/sum_pmix
warning("\'pmix\' is rescaled to add up to 1")
}
if(any(c(allpar$kappa1, allpar$kappa2) < 0)) stop("kappa1 and kappa2 must be non-negative")
if(any(allpar$mu1 < 0 | allpar$mu1 >= 2*pi)) allpar$mu1 <- prncp_reg(allpar$mu1)
if(any(allpar$mu2 < 0 | allpar$mu2 >= 2*pi)) allpar$mu2 <- prncp_reg(allpar$mu2)
out <- matrix(0, n, 2)
ncomp <- allpar_len[1] # number of components
comp_ind <- cID(tcrossprod(rep(1, n), allpar$pmix), ncomp, runif(n))
# n samples from multinom(ncomp, pmix)
for(j in seq_len(ncomp)) {
obs_ind_j <- which(comp_ind == j)
n_j <- length(obs_ind_j)
if(n_j > 0) {
out[obs_ind_j, ] <- rvmsin(n_j, kappa1[j], kappa2[j],
kappa3[j], mu1[j], mu2[j], method)
}
}
out
}
#' @rdname rvmsinmix
#' @export
dvmsinmix <- function(x, kappa1, kappa2, kappa3,
mu1, mu2, pmix, log=FALSE)
{
allpar <- list("kappa1"=kappa1, "kappa2"=kappa2, "kappa3"=kappa3,
"mu1"=mu1, "mu2"=mu2, "pmix"=pmix)
allpar_len <- listLen(allpar)
if(min(allpar_len) != max(allpar_len)) stop("component size mismatch: number of components of the input parameter vectors differ")
if(any(allpar$pmix < 0)) stop("\'pmix\' must be non-negative")
sum_pmix <- sum(allpar$pmix)
if(signif(sum_pmix, 5) != 1) {
if(sum_pmix <= 0) stop("\'pmix\' must have at least one positive element")
allpar$pmix <- allpar$pmix/sum_pmix
warning("\'pmix\' is rescaled to add up to 1")
}
if(any(c(allpar$kappa1, allpar$kappa2) < 0)) stop("kappa1 and kappa2 must be nonnegative")
if(any(allpar$mu1 < 0 | allpar$mu1 >= 2*pi)) allpar$mu1 <- prncp_reg(allpar$mu1)
if(any(allpar$mu2 < 0 | allpar$mu2 >= 2*pi)) allpar$mu2 <- prncp_reg(allpar$mu2)
if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
|| (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
ncomp <- length(kappa1)
if (length(x) == 2) x <- matrix(x, nrow=1)
allcompden <- vapply(1:ncomp,
function(j) dvmsin(x, kappa1[j], kappa2[j],
kappa3[j], mu1[j], mu2[j], FALSE),
rep(0, nrow(x)))
mixden <- c(allcompden %*% pmix)
if (log) {
log(mixden)
} else {
mixden
}
}
#' Fitting bivariate von Mises sine model mixtures using MCMC
#' @param ... arguments (other than \code{model}) passed to \link{fit_angmix}
#'
#' @details Wrapper for \link{fit_angmix} with \code{model = "vmsin"}
#'
#' @examples
#' # illustration only - more iterations needed for convergence
#' fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter = 20,
#' n.chains = 1)
#' fit.vmsin.20
#'
#' @export
fit_vmsinmix <- function(...)
{
fit_angmix(model = "vmsin", ...)
}
vmsin_var_cor_singlepar_numeric <- function(kappa1, kappa2, kappa3) {
fn_vmsin_const <- function(pars) {
const_vmsin(pars[1], pars[2], pars[3])
}
const <- fn_vmsin_const(c(kappa1, kappa2, kappa3))
grad <- numDeriv::grad(fn_vmsin_const, c(kappa1, kappa2, kappa3))
names(grad) <- c("k1", "k2", "k3")
hess <- numDeriv::hessian(fn_vmsin_const, c(kappa1, kappa2, kappa3))
dimnames(hess) <- list(c("k1", "k2", "k3"), c("k1", "k2", "k3"))
rho_fl <- unname(
(grad["k3"] * hess["k1", "k2"]) /
sqrt(
hess["k1", "k1"] * (const - hess["k1", "k1"])
* hess["k2", "k2"] * (const - hess["k2", "k2"])
)
)
rho_js <- unname(
grad["k3"] /
sqrt((const - hess["k1", "k1"]) * (const - hess["k2", "k2"]))
)
var1 <- unname(1 - grad["k1"]/const)
var2 <- unname(1 - grad["k2"]/const)
# if (rho_js <= -1) {
# rho_js <- -1
# } else if (rho_js >= 1) {
# rho_js <- 1
# }
#
# dat <- rvmsin(N, kappa1, kappa2, kappa3, 0, 0)
#
# ave_sin1sin2 <- sum(sin(dat[, 1]) * sin(dat[, 2]))/N
# ave_cos1cos2 <- sum(cos(dat[, 1]) * cos(dat[, 2]))/N
#
# ave_sin1sq <- sum(sin(dat[, 1])^2)/N
# ave_cos1sq <- 1-ave_sin1sq
# ave_cos1 <- sum(cos(dat[, 1]))/N
#
# ave_sin2sq <- sum(sin(dat[, 2])^2)/N
# ave_cos2sq <- 1-ave_sin2sq
# ave_cos2 <- sum(cos(dat[, 2]))/N
#
# rho_js <- ave_sin1sin2/sqrt(ave_sin1sq * ave_sin2sq)
# # ifelse(ave_sin1sin2 >= 0, 1, -1) *
# # min(abs(ave_sin1sin2)/sqrt(ave_sin1sq * ave_sin2sq), 1)
#
# rho_fl <- rho_js *
# ave_cos1cos2/sqrt(ave_cos1sq * ave_cos2sq)
# # ifelse(ave_cos1cos2 >= 0, 1, -1) *
# # min(abs(ave_cos1cos2)/sqrt(ave_cos1sq * ave_cos2sq), 1)
#
# var1 <- min(1 - ave_cos1, 1)
# var2 <- min(1 - ave_cos2, 1)
list(var1 = var1, var2 = var2, rho_fl = rho_fl, rho_js = rho_js)
}
vmsin_var_cor_singlepar <- function(kappa1, kappa2, kappa3, ...) {
if (max(kappa1, kappa2, abs(kappa3) > 150)) {
out <- vmsin_var_cor_singlepar_numeric(kappa1, kappa2, kappa3)
} else {
out <- vmsin_var_corr_anltc(kappa1, kappa2, kappa3)
}
for (rho in c("rho_js", "rho_fl")) {
if (out[[rho]] <= -1) {
out[[rho]] <- -1
} else if (out[[rho]] >= 1) {
out[[rho]] <- 1
}
}
for (var in c("var1", "var2")) {
if (out[[var]] <= 0) {
out[[var]] <- 0
} else if (out[[var]] >= 1) {
out[[var]] <- 1
}
}
out
}
rvmsin_1par <- function(n=1, kappa1=1, kappa2=1, kappa3=1,
mu1=0, mu2=0, method = "vmprop")
{
if (abs(abs(kappa3) - kappa1) < 1e-8 |
abs(abs(kappa3) - kappa2) < 1e-8)
kappa3 <- kappa3*(1+1e-6)
if (method == "vmprop") {
log_const_vmsin <- log(const_vmsin(kappa1, kappa2, kappa3))
log_2pi <- log(2*pi)
unimodal_x <- TRUE
phistar <- NULL
# first check if the joint density is bimodal
if (kappa3^2 >= kappa1*kappa2) {
# now check if the x marginal is bimodal
if (A_bessel(kappa2) > kappa1*kappa2/kappa3^2) {
unimodal_x <- FALSE
phistar_eqn <- function(phi) {
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi)^2)
cos(phi)*A_bessel(a_phi)/a_phi - kappa1/kappa3^2
}
# find the modes in (-pi, pi)
find_root <- uniroot.all(phistar_eqn, c(-pi, pi))
# if no root found, use the naive two dimensional rejection sampler
if (is.null(find_root)) {
method <- "naive"
} else {
phistar <- find_root[1]
}
}
}
if (method=="vmprop" & unimodal_x) {
grid_0_2pi <- seq(0, 2*pi, length.out = 100)
# do minimax to find optimum kappa and K
# abs difference between the target log marginal den
# and log proposal den
obj_kappa_abs <- function(kappa) {
obj_kappa_phi <- function(phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <- -log_2pi - log(besselI(kappa, 0)) + kappa*cos_phi_mu1
abs(log_target_marginal - log_proposal)
}
optimize(obj_kappa_phi, c(0, 2*pi), maximum=TRUE)$objective
}
# minimize obj_kappa_abs wrt kappa
kappa_opt <- optimize(obj_kappa_abs, c(0, max(50, kappa2*2)))$minimum
log_prop_den <- function(kappa, phi) {
cos_phi <- cos(phi)
-log_2pi - log(besselI(kappa, 0)) +
kappa*cos_phi
}
# find max of log_prop_den
max_log_prop_den <- max(log_prop_den(kappa_opt, grid_0_2pi))
# any point with density < exp(max_log_prop_den)*exp(-10)
# are less likely to appear, so make exp(max_log_prop_den)*exp(-10)
# as the lower bound for propsal density to avoid instability
# browser()
# now maximize the log ratio of the two densities
# w.r.t phi, given kappa = kappa_opt
obj_kappa_phi <- function(kappa, phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <- pmax(-log_2pi - log(besselI(kappa, 0)) +
kappa*cos_phi_mu1, max_log_prop_den-10)
# log_proposals are >= max_log_prop_den -10 to avoid instability
# also, proposal deviates with log density less than this bound are unlikely
log_target_marginal - log_proposal
}
# maximize the log-ratio over a grid
(logK <- max(obj_kappa_phi(kappa_opt, seq(0, 2*pi, length.out = 200))))
# cat(exp(-logK))
# rcos_unimodal(1, kappa1, kappa2, kappa3, mu1, mu2,
# kappa_opt, log(besselI(kappa_opt, 0, TRUE)) + kappa_opt,
# logK, log_const_vmsin)
rsin_unimodal(n, kappa1, kappa2, kappa3, mu1, mu2,
kappa_opt, log(BESSI0_C(kappa_opt)), logK,
log_const_vmsin)
}
else if (method=="vmprop" & !unimodal_x) {
# browser()
# do all otimizations is (-pi, pi) then revert back to [0, 2*pi]
# change the modes into [0, 2*pi]
mode_1 <- prncp_reg(mu1 + phistar)
mode_2 <- prncp_reg(mu1 - phistar)
# sin_phistar <- sin(phistar)
unifpropn <- 1e-10
vmpropn <- (1-1e-10)/2
grid_0_2pi <- seq(0, 2*pi, length.out = 100)
# proposal for y marginal = vmpropn-vmpropn-unifpropn mix
# of vm(kappa, mu2 +- phistar) and unif(0, 2*pi)
# and unif(0, 2*pi) (to avoid overflow of the ratio where
# the densities are flat)
# do minimax to find optimum kappa and K
# don't worry about mu1 mu2 while optimizing!
# abs difference between the target log marginal den
# and log proposal den
obj_kappa_abs <- function(kappa) {
obj_kappa_phi <- function(phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <-
-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
kappa*cos(phi-mode_1) +
log(1+exp(kappa*(cos(phi-mode_2)
- cos(phi-mode_1))))
) + unifpropn)
# simplified mixture density
abs(log_target_marginal - log_proposal)
}
# optimize(obj_kappa_phi, c(0, 2*pi), maximum=TRUE)$objective
max(obj_kappa_phi(grid_0_2pi))
}
# minimize obj_kappa_abs wrt kappa
(kappa_opt <- optimize(obj_kappa_abs,
c(0, max(kappa1, kappa2,
abs(kappa3))))$minimum)
# browser()
#
# target_den <- function(phi) {
# cos_phi_mu1 <- cos(phi-mu1)
# a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
# log_I_a_phi <- log(besselI(a_phi, 0))
# log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
# exp(log_target_marginal)
# }
log_prop_den <- function(kappa, phi) {
cos_phi <- cos(phi)
-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
kappa*cos(phi-mode_1) +
log(1+exp(kappa*(cos(phi-mode_2)
- cos(phi-mode_1))))
) + unifpropn)
}
# find max of log_prop_den
max_log_prop_den <- max(log_prop_den(kappa_opt, grid_0_2pi))
# any point with density < exp(max_log_prop_den)*exp(-10)
# are less likely to appear, so make exp(max_log_prop_den)*exp(-10)
# as the lower bound for propsal density to avoid instability
# now maximize the log ratio of the two densities
# w.r.t phi, given kappa = kappa_opt
obj_kappa_phi <- function(kappa, phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <-
pmax(
-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
kappa*cos(phi-mode_1) +
log(1+exp(kappa*(cos(phi-mode_2)
- cos(phi-mode_1))))
) + unifpropn),
max_log_prop_den-10)
# log_proposals are >= max_log_prop_den -10 to avoid instability
# also, proposal deviates with log density less than this bound are unlikely
log_target_marginal - log_proposal
}
#
# prop_den <- function(kappa, phi) {
# exp(log_prop_den(kappa, phi))
# }
#
# pp <- function(kappa, phi) {
# cos_phi <- cos(phi)
# log_proposal <-
# pmax(-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
# kappa*cos(phi-mode_1) +
# log(1+exp(kappa*(cos(phi-mode_2) - cos(phi-mode_1)))))
# + unifpropn), max_log_prop_den-10)
# exp(log_proposal)
# }
# #
# integrate(target_den, 0, 2*pi)
# integrate(function(x) prop_den(kappa_opt, x), 0, 2*pi)
# maximize the log-ratio over a grid
(logK <- max(obj_kappa_phi(kappa_opt, seq(0, 2*pi, length.out = 200))))
exp(logK)
# poin <- seq(0, 2*pi, length.out = 100)
# # ggp <- target_den(poin)
# # ppp <- prop_den(kappa_opt, poin)
# # #
# # plot(poin, ggp, type="l", ylim=range(ggp, ppp, ggp/ppp))
# # points(poin, ppp, type="l", col = "blue")
# # points(poin, ggp/ppp, type="l")
#
#
# lggp <- log(target_den(poin))
# lppp <- log(prop_den(kappa_opt, poin))
# difp <- obj_kappa_phi(kappa_opt, poin)
# #
# plot(poin, lggp, type="l", ylim=range(lggp, lppp, lggp-lppp))
# points(poin, lppp, type="l", col = "blue")
# points(poin, lggp-lppp, type="l")
# points(poin, difp, type="l", col = "red")
# abline(v=c(mode_1, mode_2))
# #
# #
# #
# # browser()
# cat(exp(-logK))
rsin_bimodal(n, kappa1, kappa2, kappa3, mu1, mu2,
kappa_opt, log(BESSI0_C(kappa_opt)), logK,
log_const_vmsin, mode_1, mode_2,
vmpropn, unifpropn)
}
}
else if (method == "naive") {
opt_obj <- function(k1=1, k2=1, k3=0, mu1=0, mu2=0) {
# for numerical stability, if k3 < 0, fabs(k1+k3) < 1e-5 or fabs(k2+k3) < 1e-5
# make k3 = k3 * (1+1e-5)
if (k3 < 0) {
while (abs(k1 + k3) < 1e-5 || abs(k2 + k3) < 1e-5) {
k3 = k3*(1+1e-5)
}
}
obj <- optim(c(0,0), fn = function(x) -(k1*cos(x[1]-mu1)+k2*cos(x[2]-mu2)+k3*sin(x[1]-mu1)*sin(x[2]-mu2)),
gr = function(x) -c(-k1*sin(x[1]-mu1)+k3*cos(x[1]-mu1)*sin(x[2]-mu2),
-k2*sin(x[2]-mu2)+k3*sin(x[1]-mu1)*cos(x[2]-mu2)))
-obj$value
}
upper_bd <- opt_obj(kappa1, kappa2, kappa3, mu1, mu2)
rsin_onepar(n, kappa1, kappa2, kappa3, mu1, mu2, upper_bd)
}
}
Try the BAMBI package in your browser
Any scripts or data that you put into this service are public.
BAMBI documentation built on March 31, 2023, 11:24 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions rvmsin_1par vmsin_var_cor_singlepar vmsin_var_cor_singlepar_numeric fit_vmsinmix dvmsinmix rvmsinmix dvmsin rvmsin
Documented in dvmsin dvmsinmix fit_vmsinmix rvmsin rvmsinmix
#' The bivariate von Mises sine model
#'
#' @param 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.
#' @param 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.
#' @param mu1,mu2 vectors of mean parameters.
#' @param kappa1,kappa2,kappa3 vectors of concentration parameters; \code{kappa1, kappa2 > 0}.
#' @param log logical. Should the log density be returned instead?
#' @param method Rejection sampling method to be used. Available choices are \code{"naive"} (default) or \code{"vmprop"}. See details.
#' @details
#' The bivariate von Mises sine model density at the point \eqn{x = (x_1, x_2)} is given by
#' \deqn{f(x) = C_s (\kappa_1, \kappa_2, \kappa_3) \exp(\kappa_1 \cos(T_1) + \kappa_2 \cos(T_2) + \kappa_3 \sin(T_1) \sin(T_2))}
#' where
#' \deqn{T_1 = x_1 - \mu_1; T_2 = x_2 - \mu_2}
#' and \eqn{C_s (\kappa_1, \kappa_2, \kappa_3)} denotes the normalizing constant for the sine model.
#'
#' Two different rejection sampling methods are implemented for random generation. If \code{method = "vmprop"}, then first the y-marginal
#' is drawn from the associated marginal density, and then x is generated from the conditional distributio of x given y. The marginal generation of
#' y is implemented in a rejection sampling scheme with proposal being either von Mises (if the target marginal density is unimodal), or a mixture of
#' von Mises (if bimodal), with optimally chosen concentration. This the method suggested in Mardia et al. (2007). On the other hand, when
#' \code{method = "naive"} (default) a (naive) bivariate rejection sampling scheme with (bivariate) uniform propsoal is used.
#'
#' Note that although method = \code{"vmprop"} may provide better efficiency when the density is highly concentrated, it does have
#' an (often substantial) overhead due to the optimziation step required to find a reasonable proposal concentration parameter.
#' This can compensate the efficiency gains of this method, especially when \code{n} is not large.
#'
#'
#'
#'
#' @return \code{dvmsin} gives the density and \code{rvmsin} 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,
#' # dvmsin returns single density
#' dvmsin(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
#' dvmsin(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
#' # dvmsin returns a vector of with ith element being the density
#' # evaluated at x with parameter values kappa1[i], kappa2[i],
#' # kappa3[i], mu1[i] and mu2[i]
#' dvmsin(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 dvmsin returns a vector of 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]
#' dvmsin(x[1, ], kappa1, kappa2, kappa3, mu1, mu2)
#'
#' # when parameters are all scalars, number of observations generated
#' # by rvmsin is n
#' rvmsin(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 rvmsin is the same as the length of the
#' # recycled vectors
#' rvmsin(n, kappa1, kappa2, kappa3, mu1, mu2)
#'
#' @export
rvmsin <- function(n, kappa1=1, kappa2=1,
kappa3=0, mu1=0, mu2=0, method="naive")
{
if(any(c(kappa1, kappa2) < 0)) stop("kappa1 and kappa2 must be nonnegative")
if(any(mu1 < 0 | mu1 >= 2*pi)) mu1 <- prncp_reg(mu1)
if(any(mu2 < 0 | mu2 >= 2*pi)) mu2 <- prncp_reg(mu2)
if(n < 0) stop("invalid n")
# if (is.null(method)) {
# if (n > 1e5) method <- "vmprop"
# else method <- "naive"
# }
if (!method %in% c("naive", "vmprop"))
stop("method must be either \'naive\' or \'vmprop\'")
if(max(length(kappa1), length(kappa2), length(kappa3),
length(mu1), length(mu2)) > 1) {
expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
k1 <- expanded[[1]]; k2 <- expanded[[2]]; k3 <- expanded[[3]]
mu1 <- expanded[[4]]; mu2 <- expanded[[5]]
t(vapply(1:length(k1),
function(j) rvmsin_1par(1, k1[j], k2[j], k3[j],
mu1[j], mu2[j], "naive"), c(0, 0)))
} else {
# if (is.null(method) & max(kappa1, kappa2, abs(kappa3)) < 0.1)
# method <- "vmprop"
rvmsin_1par(n, kappa1, kappa2, kappa3, mu1, mu2, method)
}
}
# rvmsin_naive <- function(n, kappa1=1, kappa2=1,
# kappa3=0, mu1=0, mu2=0) {
#
# opt_obj <- function(k1=1, k2=1, k3=0, mu1=0, mu2=0) {
# # for numerical stability, if k3 < 0, fabs(k1+k3) < 1e-5 or fabs(k2+k3) < 1e-5
# # make k3 = k3 * (1+1e-5)
# if (k3 < 0) {
# while (abs(k1 + k3) < 1e-5 || abs(k2 + k3) < 1e-5) {
# k3 = k3*(1+1e-5)
# }
# }
# obj <- optim(c(0,0), fn = function(x) -(k1*cos(x[1]-mu1)+k2*cos(x[2]-mu2)+k3*sin(x[1]-mu1)*sin(x[2]-mu2)),
# gr = function(x) -c(-k1*sin(x[1]-mu1)+k3*cos(x[1]-mu1)*sin(x[2]-mu2),
# -k2*sin(x[2]-mu2)+k3*sin(x[1]-mu1)*cos(x[2]-mu2)))
# -obj$value
#
# }
#
# if(max(length(kappa1), length(kappa2), length(kappa3), length(mu1), length(mu2)) > 1) {
# expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
# k1 <- expanded[[1]]; k2 <- expanded[[2]]; k3 <- expanded[[3]]
# mu1 <- expanded[[4]]; mu2 <- expanded[[5]]
# upper_bd_all <- vapply(1:length(k1),
# function(h) opt_obj(k1[h], k2[h], k3[h], mu1[h], mu2[h]),
# 0)
# rsin_manypar(k1, k2, k3, mu1, mu2, upper_bd_all)
# } else {
# upper_bd <- opt_obj(kappa1, kappa2, kappa3, mu1, mu2)
# rsin_onepar(n, kappa1, kappa2, kappa3, mu1, mu2, upper_bd)
# }
# }
#' @rdname rvmsin
#' @export
dvmsin <- function(x, kappa1=1, kappa2=1, kappa3=0, mu1=0, mu2=0, log = FALSE)
{
if(any(c(kappa1, kappa2) < 0)) stop("kappa1 and kappa2 must be non-negative")
if(any(mu1 < 0 | mu1 >= 2*pi)) mu1 <- prncp_reg(mu1)
if(any(mu2 < 0 | mu2 >= 2*pi)) mu2 <- prncp_reg(mu2)
if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
|| (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
if(max(length(kappa1), length(kappa2), length(kappa3), length(mu1), length(mu2)) > 1) {
expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
kappa1 <- expanded[[1]]
kappa2 <- expanded[[2]]
kappa3 <- expanded[[3]]
mu1 <- expanded[[4]]
mu2 <- expanded[[5]]
}
par.mat <- rbind(kappa1, kappa2, kappa3, mu1, mu2)
n_par <- ncol(par.mat)
if (length(x) == 2) x <- matrix(x, nrow = 1)
n_x <- nrow(x)
if (n_par == 1) {
log_den <- c(ldsin_manyx_onepar(x, kappa1, kappa2, kappa3,
mu1, mu2))
} else if (n_x == 1) {
log_den <- c(ldsin_onex_manypar(c(x), kappa1, kappa2, kappa3,
mu1, mu2))
} else {
x_set <- 1:nrow(x)
par_set <- 1:n_par
expndn_set <- expand_args(x_set, par_set)
x_set <- expndn_set[[1]]
par_set <- expndn_set[[2]]
log_den <- c(ldsin_manyx_manypar(x[x_set, ], kappa1[par_set],
kappa2[par_set], kappa3[par_set],
mu1[par_set], mu2[par_set]))
}
if (log) log_den
else exp(log_den)
}
#'
# dvmsin <- function(x, kappa1=1, kappa2=1, kappa3=0, mu1=0, mu2=0)
# {
# if(any(c(kappa1, kappa2) < 0)) stop("kappa1 and kappa2 must be non-negative")
# if(any(mu1 < 0 | mu1 >= 2*pi)) mu1 <- prncp_reg(mu1)
# if(any(mu2 < 0 | mu2 >= 2*pi)) mu2 <- prncp_reg(mu2)
# if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
# || (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
#
# if(max(length(kappa1), length(kappa2), length(kappa3), length(mu1), length(mu2)) > 1) {
# expanded <- expand_args(kappa1, kappa2, kappa3, mu1, mu2)
# kappa1 <- expanded[[1]]
# kappa2 <- expanded[[2]]
# kappa3 <- expanded[[3]]
# mu1 <- expanded[[4]]
# mu2 <- expanded[[5]]
# }
#
# par.mat <- rbind(kappa1, kappa2, kappa3, mu1, mu2)
# n_par <- ncol(par.mat)
# if (length(x) == 2) x <- matrix(x, nrow = 1)
# n_x <- nrow(x)
#
#
# if (n_par == 1) {
# log_den <- c(ldsin_manyx_onepar(x, kappa1, kappa2, kappa3,
# mu1, mu2))
#
# } else if (n_x == 1) {
# log_den <- c(ldsin_onex_manypar(c(x), kappa1, kappa2, kappa3,
# mu1, mu2))
# } else {
# x_set <- 1:nrow(x)
# par_set <- 1:n_par
# expndn_set <- expand_args(x_set, par_set)
# x_set <- expndn_set[[1]]
# par_set <- expndn_set[[2]]
# log_den <- c(ldsin_manyx_manypar(x[x_set, ], kappa1[par_set],
# kappa2[par_set], kappa3[par_set],
# mu1[par_set], mu2[par_set]))
# }
#
# if (log) log_den
# else exp(log_den)
# }
#' The bivariate von Mises sine model mixtures
#' @inheritParams rvmsin
#' @param n number of observations.
#' @param x matrix of angles (in radians) where the density is to be evaluated, with each row being a
#' single bivariate vector of angles.
#'
#' @param pmix vector of mixture proportions.
#' @param mu1,mu2 vectors of mean parameters.
#' @param kappa1,kappa2,kappa3 vectors of concentration parameters; \code{kappa1, kappa2 > 0} for each component.
#'
#' @details All the argument vectors \code{pmix, kappa1, kappa2, kappa3, mu1} and \code{mu2} must be of
#' the same length ( = component size of the mixture model), with \eqn{j}-th element corresponding to the
#' \eqn{j}-th component of the mixture distribution.
#' @details The bivariate von Mises sine model mixture distribution with component size \code{K = \link{length}(p.mix)} has density
#' \deqn{g(x) = \sum p[j] * f(x; \kappa_1[j], \kappa_2[j], \kappa_3[j], \mu_1[j], \mu_2[j])}
#' where the sum extends over \eqn{j}; \eqn{p[j]; \kappa_1[j], \kappa_2[j], \kappa_3[j]}; and \eqn{\mu_1[j], \mu_2[j]} respectively denote the mixing proportion,
#' the three concentration parameters and the two mean parameter for the \eqn{j}-th component, \eqn{j = 1, ..., K},
#' and \eqn{f(. ; \kappa_1, \kappa_2, \kappa_3, \mu_1, \mu_2)} denotes the density function of the von Mises sine model
#' with concentration parameters \eqn{\kappa_1, \kappa_2, \kappa_3} and mean parameters \eqn{\mu_1, \mu_2}.
#'
#' @return \code{dvmsinmix} computes the density (vector if x is a two column matrix with more than one row)
#' and \code{rvmsinmix} generates random deviates from the mixture density.
#'
#' @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)
#' pmix <- c(0.3, 0.4, 0.3)
#' x <- diag(2, 2)
#' n <- 10
#'
#' # mixture densities calculated at the rows of x
#' dvmsinmix(x, kappa1, kappa2, kappa3, mu1, mu2, pmix)
#'
#' # number of observations generated from the mixture distribution is n
#' rvmsinmix(n, kappa1, kappa2, kappa3, mu1, mu2, pmix)
#'
#' @export
rvmsinmix <- function(n, kappa1, kappa2, kappa3, mu1, mu2, pmix, method="naive")
{
allpar <- list(kappa1=kappa1, kappa2=kappa2, kappa3=kappa3,
mu1=mu1, mu2=mu2, pmix=pmix)
allpar_len <- listLen(allpar)
if(min(allpar_len) != max(allpar_len))
stop("component size mismatch: number of components of the input parameter vectors differ")
if(any(allpar$pmix < 0)) stop("\'pmix\' must be non-negative")
sum_pmix <- sum(allpar$pmix)
if(signif(sum_pmix, 5) != 1) {
if(sum_pmix <= 0) stop("\'pmix\' must have at least one positive element")
allpar$pmix <- allpar$pmix/sum_pmix
warning("\'pmix\' is rescaled to add up to 1")
}
if(any(c(allpar$kappa1, allpar$kappa2) < 0)) stop("kappa1 and kappa2 must be non-negative")
if(any(allpar$mu1 < 0 | allpar$mu1 >= 2*pi)) allpar$mu1 <- prncp_reg(allpar$mu1)
if(any(allpar$mu2 < 0 | allpar$mu2 >= 2*pi)) allpar$mu2 <- prncp_reg(allpar$mu2)
out <- matrix(0, n, 2)
ncomp <- allpar_len[1] # number of components
comp_ind <- cID(tcrossprod(rep(1, n), allpar$pmix), ncomp, runif(n))
# n samples from multinom(ncomp, pmix)
for(j in seq_len(ncomp)) {
obs_ind_j <- which(comp_ind == j)
n_j <- length(obs_ind_j)
if(n_j > 0) {
out[obs_ind_j, ] <- rvmsin(n_j, kappa1[j], kappa2[j],
kappa3[j], mu1[j], mu2[j], method)
}
}
out
}
#' @rdname rvmsinmix
#' @export
dvmsinmix <- function(x, kappa1, kappa2, kappa3,
mu1, mu2, pmix, log=FALSE)
{
allpar <- list("kappa1"=kappa1, "kappa2"=kappa2, "kappa3"=kappa3,
"mu1"=mu1, "mu2"=mu2, "pmix"=pmix)
allpar_len <- listLen(allpar)
if(min(allpar_len) != max(allpar_len)) stop("component size mismatch: number of components of the input parameter vectors differ")
if(any(allpar$pmix < 0)) stop("\'pmix\' must be non-negative")
sum_pmix <- sum(allpar$pmix)
if(signif(sum_pmix, 5) != 1) {
if(sum_pmix <= 0) stop("\'pmix\' must have at least one positive element")
allpar$pmix <- allpar$pmix/sum_pmix
warning("\'pmix\' is rescaled to add up to 1")
}
if(any(c(allpar$kappa1, allpar$kappa2) < 0)) stop("kappa1 and kappa2 must be nonnegative")
if(any(allpar$mu1 < 0 | allpar$mu1 >= 2*pi)) allpar$mu1 <- prncp_reg(allpar$mu1)
if(any(allpar$mu2 < 0 | allpar$mu2 >= 2*pi)) allpar$mu2 <- prncp_reg(allpar$mu2)
if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
|| (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
ncomp <- length(kappa1)
if (length(x) == 2) x <- matrix(x, nrow=1)
allcompden <- vapply(1:ncomp,
function(j) dvmsin(x, kappa1[j], kappa2[j],
kappa3[j], mu1[j], mu2[j], FALSE),
rep(0, nrow(x)))
mixden <- c(allcompden %*% pmix)
if (log) {
log(mixden)
} else {
mixden
}
}
#' Fitting bivariate von Mises sine model mixtures using MCMC
#' @param ... arguments (other than \code{model}) passed to \link{fit_angmix}
#'
#' @details Wrapper for \link{fit_angmix} with \code{model = "vmsin"}
#'
#' @examples
#' # illustration only - more iterations needed for convergence
#' fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter = 20,
#' n.chains = 1)
#' fit.vmsin.20
#'
#' @export
fit_vmsinmix <- function(...)
{
fit_angmix(model = "vmsin", ...)
}
vmsin_var_cor_singlepar_numeric <- function(kappa1, kappa2, kappa3) {
fn_vmsin_const <- function(pars) {
const_vmsin(pars[1], pars[2], pars[3])
}
const <- fn_vmsin_const(c(kappa1, kappa2, kappa3))
grad <- numDeriv::grad(fn_vmsin_const, c(kappa1, kappa2, kappa3))
names(grad) <- c("k1", "k2", "k3")
hess <- numDeriv::hessian(fn_vmsin_const, c(kappa1, kappa2, kappa3))
dimnames(hess) <- list(c("k1", "k2", "k3"), c("k1", "k2", "k3"))
rho_fl <- unname(
(grad["k3"] * hess["k1", "k2"]) /
sqrt(
hess["k1", "k1"] * (const - hess["k1", "k1"])
* hess["k2", "k2"] * (const - hess["k2", "k2"])
)
)
rho_js <- unname(
grad["k3"] /
sqrt((const - hess["k1", "k1"]) * (const - hess["k2", "k2"]))
)
var1 <- unname(1 - grad["k1"]/const)
var2 <- unname(1 - grad["k2"]/const)
# if (rho_js <= -1) {
# rho_js <- -1
# } else if (rho_js >= 1) {
# rho_js <- 1
# }
#
# dat <- rvmsin(N, kappa1, kappa2, kappa3, 0, 0)
#
# ave_sin1sin2 <- sum(sin(dat[, 1]) * sin(dat[, 2]))/N
# ave_cos1cos2 <- sum(cos(dat[, 1]) * cos(dat[, 2]))/N
#
# ave_sin1sq <- sum(sin(dat[, 1])^2)/N
# ave_cos1sq <- 1-ave_sin1sq
# ave_cos1 <- sum(cos(dat[, 1]))/N
#
# ave_sin2sq <- sum(sin(dat[, 2])^2)/N
# ave_cos2sq <- 1-ave_sin2sq
# ave_cos2 <- sum(cos(dat[, 2]))/N
#
# rho_js <- ave_sin1sin2/sqrt(ave_sin1sq * ave_sin2sq)
# # ifelse(ave_sin1sin2 >= 0, 1, -1) *
# # min(abs(ave_sin1sin2)/sqrt(ave_sin1sq * ave_sin2sq), 1)
#
# rho_fl <- rho_js *
# ave_cos1cos2/sqrt(ave_cos1sq * ave_cos2sq)
# # ifelse(ave_cos1cos2 >= 0, 1, -1) *
# # min(abs(ave_cos1cos2)/sqrt(ave_cos1sq * ave_cos2sq), 1)
#
# var1 <- min(1 - ave_cos1, 1)
# var2 <- min(1 - ave_cos2, 1)
list(var1 = var1, var2 = var2, rho_fl = rho_fl, rho_js = rho_js)
}
vmsin_var_cor_singlepar <- function(kappa1, kappa2, kappa3, ...) {
if (max(kappa1, kappa2, abs(kappa3) > 150)) {
out <- vmsin_var_cor_singlepar_numeric(kappa1, kappa2, kappa3)
} else {
out <- vmsin_var_corr_anltc(kappa1, kappa2, kappa3)
}
for (rho in c("rho_js", "rho_fl")) {
if (out[[rho]] <= -1) {
out[[rho]] <- -1
} else if (out[[rho]] >= 1) {
out[[rho]] <- 1
}
}
for (var in c("var1", "var2")) {
if (out[[var]] <= 0) {
out[[var]] <- 0
} else if (out[[var]] >= 1) {
out[[var]] <- 1
}
}
out
}
rvmsin_1par <- function(n=1, kappa1=1, kappa2=1, kappa3=1,
mu1=0, mu2=0, method = "vmprop")
{
if (abs(abs(kappa3) - kappa1) < 1e-8 |
abs(abs(kappa3) - kappa2) < 1e-8)
kappa3 <- kappa3*(1+1e-6)
if (method == "vmprop") {
log_const_vmsin <- log(const_vmsin(kappa1, kappa2, kappa3))
log_2pi <- log(2*pi)
unimodal_x <- TRUE
phistar <- NULL
# first check if the joint density is bimodal
if (kappa3^2 >= kappa1*kappa2) {
# now check if the x marginal is bimodal
if (A_bessel(kappa2) > kappa1*kappa2/kappa3^2) {
unimodal_x <- FALSE
phistar_eqn <- function(phi) {
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi)^2)
cos(phi)*A_bessel(a_phi)/a_phi - kappa1/kappa3^2
}
# find the modes in (-pi, pi)
find_root <- uniroot.all(phistar_eqn, c(-pi, pi))
# if no root found, use the naive two dimensional rejection sampler
if (is.null(find_root)) {
method <- "naive"
} else {
phistar <- find_root[1]
}
}
}
if (method=="vmprop" & unimodal_x) {
grid_0_2pi <- seq(0, 2*pi, length.out = 100)
# do minimax to find optimum kappa and K
# abs difference between the target log marginal den
# and log proposal den
obj_kappa_abs <- function(kappa) {
obj_kappa_phi <- function(phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <- -log_2pi - log(besselI(kappa, 0)) + kappa*cos_phi_mu1
abs(log_target_marginal - log_proposal)
}
optimize(obj_kappa_phi, c(0, 2*pi), maximum=TRUE)$objective
}
# minimize obj_kappa_abs wrt kappa
kappa_opt <- optimize(obj_kappa_abs, c(0, max(50, kappa2*2)))$minimum
log_prop_den <- function(kappa, phi) {
cos_phi <- cos(phi)
-log_2pi - log(besselI(kappa, 0)) +
kappa*cos_phi
}
# find max of log_prop_den
max_log_prop_den <- max(log_prop_den(kappa_opt, grid_0_2pi))
# any point with density < exp(max_log_prop_den)*exp(-10)
# are less likely to appear, so make exp(max_log_prop_den)*exp(-10)
# as the lower bound for propsal density to avoid instability
# browser()
# now maximize the log ratio of the two densities
# w.r.t phi, given kappa = kappa_opt
obj_kappa_phi <- function(kappa, phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <- pmax(-log_2pi - log(besselI(kappa, 0)) +
kappa*cos_phi_mu1, max_log_prop_den-10)
# log_proposals are >= max_log_prop_den -10 to avoid instability
# also, proposal deviates with log density less than this bound are unlikely
log_target_marginal - log_proposal
}
# maximize the log-ratio over a grid
(logK <- max(obj_kappa_phi(kappa_opt, seq(0, 2*pi, length.out = 200))))
# cat(exp(-logK))
# rcos_unimodal(1, kappa1, kappa2, kappa3, mu1, mu2,
# kappa_opt, log(besselI(kappa_opt, 0, TRUE)) + kappa_opt,
# logK, log_const_vmsin)
rsin_unimodal(n, kappa1, kappa2, kappa3, mu1, mu2,
kappa_opt, log(BESSI0_C(kappa_opt)), logK,
log_const_vmsin)
}
else if (method=="vmprop" & !unimodal_x) {
# browser()
# do all otimizations is (-pi, pi) then revert back to [0, 2*pi]
# change the modes into [0, 2*pi]
mode_1 <- prncp_reg(mu1 + phistar)
mode_2 <- prncp_reg(mu1 - phistar)
# sin_phistar <- sin(phistar)
unifpropn <- 1e-10
vmpropn <- (1-1e-10)/2
grid_0_2pi <- seq(0, 2*pi, length.out = 100)
# proposal for y marginal = vmpropn-vmpropn-unifpropn mix
# of vm(kappa, mu2 +- phistar) and unif(0, 2*pi)
# and unif(0, 2*pi) (to avoid overflow of the ratio where
# the densities are flat)
# do minimax to find optimum kappa and K
# don't worry about mu1 mu2 while optimizing!
# abs difference between the target log marginal den
# and log proposal den
obj_kappa_abs <- function(kappa) {
obj_kappa_phi <- function(phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <-
-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
kappa*cos(phi-mode_1) +
log(1+exp(kappa*(cos(phi-mode_2)
- cos(phi-mode_1))))
) + unifpropn)
# simplified mixture density
abs(log_target_marginal - log_proposal)
}
# optimize(obj_kappa_phi, c(0, 2*pi), maximum=TRUE)$objective
max(obj_kappa_phi(grid_0_2pi))
}
# minimize obj_kappa_abs wrt kappa
(kappa_opt <- optimize(obj_kappa_abs,
c(0, max(kappa1, kappa2,
abs(kappa3))))$minimum)
# browser()
#
# target_den <- function(phi) {
# cos_phi_mu1 <- cos(phi-mu1)
# a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
# log_I_a_phi <- log(besselI(a_phi, 0))
# log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
# exp(log_target_marginal)
# }
log_prop_den <- function(kappa, phi) {
cos_phi <- cos(phi)
-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
kappa*cos(phi-mode_1) +
log(1+exp(kappa*(cos(phi-mode_2)
- cos(phi-mode_1))))
) + unifpropn)
}
# find max of log_prop_den
max_log_prop_den <- max(log_prop_den(kappa_opt, grid_0_2pi))
# any point with density < exp(max_log_prop_den)*exp(-10)
# are less likely to appear, so make exp(max_log_prop_den)*exp(-10)
# as the lower bound for propsal density to avoid instability
# now maximize the log ratio of the two densities
# w.r.t phi, given kappa = kappa_opt
obj_kappa_phi <- function(kappa, phi) {
cos_phi_mu1 <- cos(phi-mu1)
a_phi <- sqrt(kappa2^2 + kappa3^2*sin(phi-mu1)^2)
log_I_a_phi <- log(besselI(a_phi, 0))
log_target_marginal <- -log_const_vmsin + log_2pi + log_I_a_phi + kappa2*cos_phi_mu1
log_proposal <-
pmax(
-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
kappa*cos(phi-mode_1) +
log(1+exp(kappa*(cos(phi-mode_2)
- cos(phi-mode_1))))
) + unifpropn),
max_log_prop_den-10)
# log_proposals are >= max_log_prop_den -10 to avoid instability
# also, proposal deviates with log density less than this bound are unlikely
log_target_marginal - log_proposal
}
#
# prop_den <- function(kappa, phi) {
# exp(log_prop_den(kappa, phi))
# }
#
# pp <- function(kappa, phi) {
# cos_phi <- cos(phi)
# log_proposal <-
# pmax(-log_2pi + log(exp(log(vmpropn) - log(besselI(kappa, 0)) +
# kappa*cos(phi-mode_1) +
# log(1+exp(kappa*(cos(phi-mode_2) - cos(phi-mode_1)))))
# + unifpropn), max_log_prop_den-10)
# exp(log_proposal)
# }
# #
# integrate(target_den, 0, 2*pi)
# integrate(function(x) prop_den(kappa_opt, x), 0, 2*pi)
# maximize the log-ratio over a grid
(logK <- max(obj_kappa_phi(kappa_opt, seq(0, 2*pi, length.out = 200))))
exp(logK)
# poin <- seq(0, 2*pi, length.out = 100)
# # ggp <- target_den(poin)
# # ppp <- prop_den(kappa_opt, poin)
# # #
# # plot(poin, ggp, type="l", ylim=range(ggp, ppp, ggp/ppp))
# # points(poin, ppp, type="l", col = "blue")
# # points(poin, ggp/ppp, type="l")
#
#
# lggp <- log(target_den(poin))
# lppp <- log(prop_den(kappa_opt, poin))
# difp <- obj_kappa_phi(kappa_opt, poin)
# #
# plot(poin, lggp, type="l", ylim=range(lggp, lppp, lggp-lppp))
# points(poin, lppp, type="l", col = "blue")
# points(poin, lggp-lppp, type="l")
# points(poin, difp, type="l", col = "red")
# abline(v=c(mode_1, mode_2))
# #
# #
# #
# # browser()
# cat(exp(-logK))
rsin_bimodal(n, kappa1, kappa2, kappa3, mu1, mu2,
kappa_opt, log(BESSI0_C(kappa_opt)), logK,
log_const_vmsin, mode_1, mode_2,
vmpropn, unifpropn)
}
}
else if (method == "naive") {
opt_obj <- function(k1=1, k2=1, k3=0, mu1=0, mu2=0) {
# for numerical stability, if k3 < 0, fabs(k1+k3) < 1e-5 or fabs(k2+k3) < 1e-5
# make k3 = k3 * (1+1e-5)
if (k3 < 0) {
while (abs(k1 + k3) < 1e-5 || abs(k2 + k3) < 1e-5) {
k3 = k3*(1+1e-5)
}
}
obj <- optim(c(0,0), fn = function(x) -(k1*cos(x[1]-mu1)+k2*cos(x[2]-mu2)+k3*sin(x[1]-mu1)*sin(x[2]-mu2)),
gr = function(x) -c(-k1*sin(x[1]-mu1)+k3*cos(x[1]-mu1)*sin(x[2]-mu2),
-k2*sin(x[2]-mu2)+k3*sin(x[1]-mu1)*cos(x[2]-mu2)))
-obj$value
}
upper_bd <- opt_obj(kappa1, kappa2, kappa3, mu1, mu2)
rsin_onepar(n, kappa1, kappa2, kappa3, mu1, mu2, upper_bd)
}
}
Try the BAMBI package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.