R/all_wnorm_fns.R

In c7rishi/BAMBI: Bivariate Angular Mixture Models

# list of functions


#' The univariate Wrapped Normal distribution
#' @inheritParams rvm
#' @param int.displ integer displacement. If \code{int.displ =} M, then the infinite sum in the
#' density is approximated by a sum over 2*M + 1 elements. (See Details.) The allowed values are 1, 2, 3, 4 and 5. Default is 3.
#' @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 wrapped normal distribution has density
#' \deqn{f(x) = \sqrt(\kappa/(2\pi)) \sum \exp(-\kappa/2 (x - \mu(2\pi\omega))^2)}
#' where the sum extends over all integers \eqn{\omega},
#' and is approximated by a sum over \eqn{\omega} in \eqn{\{-M, -M+1, ..., M-1, M \}} if \code{int.displ = } \eqn{M}.
#' @return \code{dwnorm} gives the density  and \code{rwnorm} generates random deviates.
#'
#' @examples
#'
#' kappa <- 1:3
#' mu <- 0:2
#' x <- 1:10
#' n <- 10
#'
#'
#' # when x and both parameters are scalars, dwnorm returns a single density
#' dwnorm(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
#' dwnorm(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 dwnorm returns a vector of with ith element being the
#' # density evaluated at x with parameter values kappa[i] and mu[i]
#' dwnorm(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 dwnorm 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]
#' dwnorm(x, kappa, mu)
#'
#' # when parameters are all scalars, number of observations generated by rwnorm is n
#' rwnorm(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 rwnorm is the same as the length
#' # of the recycled vectors
#' rwnorm(n, kappa, mu)
#'
#' @export

rwnorm <- function(n=1, 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 (all(kappa > 1e-10)) {
    samp <- rnorm(n, mean = mu, sd = sqrt(1/kappa))
    prncp_reg(samp)
  } else {
    expndn_set <- expand_args(kappa, mu, 1:n)
    kappa <- expndn_set[[2]]
    mu <- expndn_set[[3]]
    # den <- rep(0, n_x)
    n_final <- length(kappa)
    which_unif <- which(kappa < 1e-10)
    n_unif <- length(which_unif)
    samp <- rep(0, n_final)
    samp[which_unif] <- runif(n_unif, 0, 2*pi)
    if (n_unif < n_final)
      samp[-which_unif] <- prncp_reg(rnorm(n_final - n_unif, mu[-which_unif],
                                           sqrt(1/kappa[-which_unif])))


    # for (i in 1:n_x) {
    #   if (kappa1[par_set[i]] < 1e-10 | kappa2[par_set[i]] < 1e-10) {
    #     den[i] <- c(dwnorm2_onex_manypar(x[x_set[i], ], kappa[par_set[i]],
    #                                      kappa2[par_set[i]], kappa3[par_set[i]],
    #                                      mu1[par_set[i]], mu2[par_set[i]], omega.2pi))
    #   } else {
    #     den[i] <- 1/(4*pi^2)
    #   }
    # }
    samp
  }

}

#' @rdname rwnorm
#' @export
dwnorm <- function(x, kappa = 1, mu = 0, int.displ,
                   log=FALSE)
{
  if(any(kappa < 0)) stop("kappa must be non-negative")
  if(any(mu < 0 | mu >= 2*pi)) mu <- prncp_reg(mu)

  if(missing(int.displ)) int.displ <- 3
  else if(int.displ >= 5) int.displ <- 5
  else if(int.displ <= 1) int.displ <- 1
  displ <- floor(int.displ)
  omega.2pi.1d <- (-displ):displ * (2*pi) # 2pi * 1d integer displacements


  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)

  # browser()

  if (all(kappa > 1e-10)) {
    if (n_par == 1) {
      den <- c(duniwnorm_manyx_onepar(as.vector(x), kappa, mu, omega.2pi.1d))

    } else if (n_x == 1) {
      den <- c(duniwnorm_onex_manypar(x, kappa, mu, omega.2pi.1d))
    } 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(duniwnorm_manyx_manypar(x[x_set], kappa[par_set],
                                       mu[par_set], omega.2pi.1d))
    }
  } 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]]
    x_long <- x[x_set]
    kappa_long <- kappa[par_set]
    mu_long <- mu[par_set]
    which_unif <- which(kappa_long < 1e-10)
    n_x_long <- length(x_long)
    den <- rep(0, n_x_long)
    den[which_unif] <- 1/(2*pi)
    if (length(which_unif) < n_x_long)
      den[-which_unif] <- c(duniwnorm_manyx_manypar(x_long[-which_unif],
                                                    kappa_long[-which_unif],
                                                    mu_long[-which_unif],
                                                    omega.2pi.1d))
  }


  if (log) {
    den <- log(den)
  }
  den

}

#' The univariate Wrapped Normal mixtures
#' @inheritParams rvmmix
#' @inheritParams rwnorm
#' @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 wrapped normal 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) wrapped normal distribution with mean parameter \eqn{\mu} and concentration parameter \eqn{\kappa}.
#' @return \code{dwnormmix} computes the density  and \code{rwnormmix} 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
#' dwnormmix(x, kappa, mu, pmix)
#'
#' # number of observations generated from the mixture distribution is n
#' rwnormmix(n, kappa, mu, pmix)
#'
#' @export

rwnormmix <- function(n=1, 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, 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(allpar$kappa <= 0)) stop("kappa must be positive in wnorm")
  if(any(allpar$mu < 0 | allpar$mu >= 2*pi)) allpar$mu <- prncp_reg(allpar$mu)

  clus_label <- cID(tcrossprod(rep(1, n), allpar$pmix), length(allpar$pmix), runif(n))
  samp <- rnorm(length(clus_label), mean = allpar$mu[clus_label], sd = 1/sqrt(allpar$kappa[clus_label]))

  prncp_reg(samp)
}

#' @rdname rwnormmix
#' @export
dwnormmix <- function(x, kappa, mu, pmix, int.displ=3, 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, 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(allpar$kappa <= 0)) stop("kappa must be positive")
  if(any(allpar$mu < 0 | allpar$mu >= 2*pi)) allpar$mu <- prncp_reg(allpar$mu)


  ncomp <- length(kappa)

  allcompden <- vapply(1:ncomp,
                       function(j) dwnorm(x, kappa[j], mu[j],
                                          int.displ, FALSE),
                       rep(0, length(x)))

  mixden <- c(allcompden %*% pmix)

  if (log) {
    log(mixden)
  } else {
    mixden
  }

}


#' Fitting univariate wrapped normal mixtures using MCMC
#' @inheritParams fit_vmsinmix
#'
#' @details
#' Wrapper for \link{fit_angmix} with \code{model = "wnorm"}.
#'
#'
#' @examples
#' # illustration only - more iterations needed for convergence
#' fit.wnorm.20 <- fit_wnormmix(wind$angle, ncomp = 3, n.iter =  20,
#'                              n.chains = 1)
#' fit.wnorm.20
#'
#' @export

fit_wnormmix <- function(...)
{
  fit_angmix(model="wnorm", ...)
}