R/all_vm_fns.R

In BAMBI: Bivariate Angular Mixture Models

# 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", ...)
}