Nothing
# 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", ...)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.