# list of functions
#' The univariate von Mises distribution
#' @inheritParams rvmsin
#' @param x vector of angles (in radians) where the densities are to be evaluated.
#' @param mu vector of means.
#' @param kappa vector of concentration (inverse-variance) parameters; \code{kappa} > 0.
#' @details If \code{mu} and \code{kappa} are not specified they assume the default values of \code{0} and \code{1} respectively.
#' @details The univariate von Mises distribution has density
#' \deqn{f(x) = 1/(2\pi I_0 (\kappa)) \exp(\kappa \cos(x - mu))}
#' where \eqn{I_0 (\kappa)} denotes the modified Bessel function of the first kind with order 0 evaluated at the point \eqn{\kappa}.
#' @return \code{dvm} gives the density and \code{rvm} generates random deviates.
#'
#' @examples
#'
#' kappa <- 1:3
#' mu <- 0:2
#' x <- 1:10
#' n <- 10
#'
#'
#' # when x and both parameters are scalars, dvm returns a single density
#' dvm(x[1], kappa[1], mu[1])
#'
#' # when x is a vector but both the parameters are scalars, dmv returns a vector of
#' # densities calculated at each entry of x with the same parameters
#' dvm(x, kappa[1], mu[1])
#'
#' # if x is scalar and at least one of the two paraemters is a vector, both parameters are
#' # recycled to the same length, and dvm returns a vector of with ith element being the
#' # density evaluated at x with parameter values kappa[i] and mu[i]
#' dvm(x[1], kappa, mu)
#'
#' # if x and at least one of the two paraemters is a vector, x and the two parameters are
#' # recycled to the same length, and dvm returns a vector of with ith element being the
#' # density at ith element of the (recycled) x with parameter values kappa[i] and mu[i]
#' dvm(x, kappa, mu)
#'
#' # when parameters are all scalars, number of observations generated by rvm is n
#' rvm(n, kappa[1], mu[1])
#'
#' # when at least one of the two parameters is a vector, both are recycled to the same length,
#' # n is ignored, and the number of observations generated by rvm is the same as the length of
#' # the recycled vectors
#' rvm(n, kappa, mu)
#'
#' @export
rvm <- function(n, kappa = 1, mu = 0)
{
if(any(kappa < 0)) stop("kappa must be non_negative")
if(any(mu < 0 | mu >= 2*pi)) mu <- prncp_reg(mu)
if(max(length(kappa), length(mu)) > 1) {
expanded <- expand_args(kappa, mu)
kappa <- expanded[[1]]; mu <- expanded[[2]]
m <- length(kappa)
out <- rep(0, m)
for(j in 1:m) {
if (kappa[j] > 1e-7) {
out[j] <- c(runivm_onepar(1, kappa[j], mu[j]))
} else {
out[j] <- runif(1, 0, 2*pi)
}
}
} else {
if (kappa> 1e-10) {
out <- c(runivm_onepar(n, kappa, mu))
} else {
out <- runif(n, 0, 2*pi)
}
}
out
}
#' @rdname rvm
#' @export
dvm <- function(x, kappa = 1, mu = 0, log = FALSE)
{
if(any(kappa < 0)) stop("kappa must be non-negative")
if(any(mu < 0 | mu >= 2*pi)) mu <- prncp_reg(mu)
# if(max(length(kappa), length(mu)) > 1) {
# expanded <- expand_args(kappa, mu)
# kappa <- expanded[[1]]; mu <- expanded[[2]]
#
# if(length(x) > 1) {
# x_set <- 1:length(x)
# par_set <- 1:length(kappa)
# expndn_set <- expand_args(x_set, par_set)
# x_set <- expndn_set[[1]]
# par_set <- expndn_set[[2]]
# } else{
# den <- as.vector(dunivm_onex_manypar(x, kappa, mu))
# }
#
# } else {
# if(length(x) > 1){
# den <- as.vector(dunivm_manyx_onepar(as.vector(x), kappa, mu))
# } else{
# den <- exp(ldunivmnum(as.vector(x), c(kappa, mu))) / const_univm(kappa)
# }
# }
if(max(length(kappa), length(mu)) > 1) {
expanded <- expand_args(kappa, mu)
kappa <- expanded[[1]]; mu <- expanded[[2]]
}
par.mat <- rbind(kappa, mu)
n_par <- ncol(par.mat)
n_x <- length(x)
if (n_par == 1) {
den <- c(dunivm_manyx_onepar(as.vector(x), kappa, mu))
} else if (n_x == 1) {
den <- c(dunivm_onex_manypar(x, kappa, mu))
} else {
x_set <- 1:n_x
par_set <- 1:n_par
expndn_set <- expand_args(x_set, par_set)
x_set <- expndn_set[[1]]
par_set <- expndn_set[[2]]
den <- c(dunivm_manyx_manypar(x[x_set], kappa[par_set], mu[par_set]))
}
if (log) den <- log(den)
den
}
#' The univariate von Mises mixtures
#' @inheritParams rvm
#' @param mu vector of component means.
#' @param kappa vector of component concentration (inverse-variance) parameters, \code{kappa > 0}.
#' @param pmix vector of mixing proportions.
#' @details \code{pmix}, \code{mu} and \code{kappa} must be of the same length, with \eqn{j}-th element corresponding to the \eqn{j}-th component of the mixture distribution.
#' @details The univariate von Mises mixture distribution with component size \code{K = \link{length}(pmix)} has density
#' \deqn{g(x) = p[1] * f(x; \kappa[1], \mu[1]) + ... + p[K] * f(x; \kappa[K], \mu[K])}
#' where \eqn{p[j], \kappa[j], \mu[j]} respectively denote the mixing proportion, concentration parameter and the mean parameter for the \eqn{j}-th component
#' and \eqn{f(. ; \kappa, \mu)} denotes the density function of the (univariate) von Mises distribution with mean parameter \eqn{\mu} and concentration parameter \eqn{\kappa}.
#' @return \code{dvmmix} computes the density and \code{rvmmix} generates random deviates from the mixture density.
#'
#' @examples
#' kappa <- 1:3
#' mu <- 0:2
#' pmix <- c(0.3, 0.3, 0.4)
#' x <- 1:10
#' n <- 10
#'
#' # mixture densities calculated at each point in x
#' dvmmix(x, kappa, mu, pmix)
#'
#' # number of observations generated from the mixture distribution is n
#' rvmmix(n, kappa, mu, pmix)
#'
#' @export
rvmmix <- function(n, kappa, mu, pmix)
{
allpar <- list(kappa=kappa, mu=mu, 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, 4) != 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(allpar$kappa < 0)) stop("kappa must be non-negative")
if(any(allpar$mu < 0 | allpar$mu >= 2*pi)) allpar$mu <- prncp_reg(allpar$mu)
out <- rep(0, n)
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] <- rvm(n_j, kappa[j], mu[j])
}
}
out
}
#' @rdname rvmmix
#' @export
dvmmix <- function(x, kappa, mu, pmix, log=FALSE)
{
allpar <- list(kappa=kappa, mu=mu, 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, 4) != 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(allpar$pmix < 0)) stop("pmix must be non-negative")
if(sum(allpar$pmix) != 1) {
allpar$pmix <- allpar$pmix/sum(allpar$pmix)
warning("\'pmix\' is rescaled to add up to 1")
}
if(any(allpar$kappa < 0)) stop("kappa must be non-negative")
if(any(allpar$mu < 0 | allpar$mu >= 2*pi)) allpar$mu <- prncp_reg(allpar$mu)
ncomp <- length(kappa)
allcompden <- vapply(1:ncomp,
function(j) dvm(x, kappa[j], mu[j], FALSE),
rep(0, length(x)))
mixden <- c(allcompden %*% pmix)
if (log) {
log(mixden)
} else {
mixden
}
}
#' Fitting univariate von Mises mixtures using MCMC
#' @inheritParams fit_vmsinmix
#'
#' @details
#' Wrapper for \link{fit_angmix} with \code{model = "vm"}.
#'
#'
#'
#' @examples
#' # illustration only - more iterations needed for convergence
#' fit.vm.20 <- fit_vmmix(wind$angle, ncomp = 3, n.iter = 20,
#' n.chains = 1)
#' fit.vm.20
#' @export
fit_vmmix <- function(...)
{
fit_angmix(model="vm", ...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.