BAMBI: An R package for Bivariate Angular Mixture Models

Description

BAMBI is an R package that provides functions for fitting (using Bayesian methods) and simulating mixtures of univariate and bivariate angular distributions. Please see the reference for a detailed description of the functionalities of BAMBI.

References

Chakraborty, S., & Wong, S. W. (2021). BAMBI: An R package for fitting bivariate angular mixture models. Journal of Statistical Software, 99 (11), 1-69. doi:10.18637/jss.v099.i11

Deviance Information Criterion (DIC) for angmcmc objects

Description

Deviance Information Criterion (DIC) for angmcmc objects

Usage

DIC(object, form = 2, ...)

Arguments

Details

Given a deviance function D(\theta) = -2 log(p(y|\theta)), and an estimate \theta* = (\sum \theta_i) / N of the posterior mean E(\theta|y), where y denote the data, \theta are the unknown parameters of the model, \theta_1, ..., \theta_N are MCMC samples from the posterior distribution of \theta given y and p(y|\theta) is the likelihood function, the (form 1 of) Deviance Infomation Criterion (DIC) is defined as

DIC = 2 ( (\sum_{s=1}^N D(\theta_s)) / N - D(\theta*) )

The second form for DIC is given by

DIC = D(\theta*) - 4 \hat{var} \log p(y|\theta_s)

where for i = 1, ..., n, \hat{var} \log p(y|\theta) denotes the estimated variance of the log likelihood based on the realizations \theta_1, ..., \theta_N.

Like AIC and BIC, DIC is an asymptotic approximation for large samples, and is only valid when the posterior distribution is approximately normal.

Value

Computes the DIC for a given angmcmc object

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
DIC(fit.vmsin.20)

Add (extra) burnin and thin to angmcmc object after original run

Description

Add (extra) burnin and thin to angmcmc object after original run

Usage

add_burnin_thin(object, burnin.prop = 0, thin = 1)

Arguments

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
lpdtrace(fit.vmsin.20)
# Now add extra burn-in
fit.vmsin.20.burn <- add_burnin_thin(fit.vmsin.20, 0.3)
lpdtrace(fit.vmsin.20.burn)

Create an mcmc.list object from an angmcmc object

Description

Create an mcmc.list object from an angmcmc object

Usage

## S3 method for class 'angmcmc'
as.mcmc.list(x, ...)

Arguments

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)

# now convert it to mcmc.list
library(coda)
fit.vmsin.20.mcmc <- as.mcmc.list(fit.vmsin.20)

Convenience function for extracting angmcmc object, and the value of the model selection criterion corresponding to the best fitted model in stepwise fits

Description

Convenience function for extracting angmcmc object, and the value of the model selection criterion corresponding to the best fitted model in stepwise fits

Usage

bestmodel(step_object)

bestcriterion(step_object)

Arguments

Details

These are convenience functions; the best fitted model and the corresponding value of model selection criterion can also be directly obtained by extracting the elements "fit.best" and "crit.best" from step_object respectively. Note that ⁠bestcriterion} returns: (a) a scalar number (class = ⁠numeric⁠) if ⁠crit⁠used in original⁠fit_incremental_angmix⁠call is⁠'AIC'⁠, ⁠'BIC'or'DIC'⁠, (b) an element of class ⁠bridge⁠from package⁠bridgesamplingifcritisLOGML⁠, (c) an element of class ⁠c("waic", "loo")ifcrit = 'WAIC'⁠, and (d) an element of class ⁠c("psis_loo", "loo")ifcrit = "LOOIC"'. See documentations of these model selection criteria for more details.

Value

bestmodel returns an angmcmc object, and bestcriterion returns the corresponding value of model selection criterion for the best fitted model in step_object.

Examples

# illustration only - more iterations needed for convergence
set.seed(1)
fit.vmsin.step.15 <- fit_incremental_angmix("vmsin", tim8, start_ncomp = 1,
                                            max_ncomp = 3, n.iter = 15,
                                            n.chains = 1,
                                            crit = "WAIC")
fit.vmsin.best.15 <- bestmodel(fit.vmsin.step.15)
fit.vmsin.best.15

crit.best <- bestcriterion(fit.vmsin.step.15)
crit.best

Log Marginal Likelihood via Bridge Sampling for angmcmc objects

Description

Log Marginal Likelihood via Bridge Sampling for angmcmc objects

Usage

## S3 method for class 'angmcmc'
bridge_sampler(samples, ..., ave_over_chains = TRUE)

Arguments

Details

Marginal likelihood is calculated by first converting the angmcmc object samples to an mcmc.list object, and then by passing the resulting mcmc.list object to bridge_sampler. If variablity across multiple chains (if any) are very different, then calling bridge_sampler separately for each chain usually provides more stable results; the final log ML is computed by averaging over chain specific MLs.

Examples

library(future)
library(parallel)
# plan(multisession, gc = TRUE) # parallelize chains

set.seed(100)
MC.fit <- fit_angmix("vmsin", tim8, ncomp=3, n.iter=5000,
                     n.chains = 3)


library(bridgesampling)
bridge_sampler(MC.fit)

Sample circular correlation coefficients

Description

Sample circular correlation coefficients

Usage

circ_cor(
  x,
  type = "js",
  alternative = "two.sided",
  jackknife = FALSE,
  bootse = FALSE,
  n.boot = 100
)

Arguments

Details

circ_cor calculates the (sample) circular correlation between the columns of x. Two parametric (the Jammalamadaka-Sarma (1988, equation 2.6) form "js", and the Fisher-Lee (1983, Section 3) form "fl") and two non-parametric (two versions of Kendall's tau) correlation coefficients are considered. The first version of Kendall's tau ("tau1") is based on equation 2.1 in Fisher and Lee (1982), whereas the second version ("tau2") is computed using equations 6.7-6.8 in Zhan et al (2017).

The cost-complexity for "js", "fl", "tau2" and "tau1" are O(n), O(n^2), O(n^2) and O(n^3) respectively, where n denotes the number of rows in x. As such, for large n evaluation of "tau1" will be slow.

References

Fisher, N. I. and Lee, A. J. (1982). Nonparametric measures of angular-angular association. Biometrika, 69(2), 315-321.

Fisher, N. I. and Lee, A. J. (1983). A correlation coefficient for circular data. Biometrika, 70(2):327-332.

Jammalamadaka, S. R. and Sarma, Y. (1988). A correlation coefficient for angular variables. Statistical theory and data analysis II, pages 349-364.

Zhan, X., Ma, T., Liu, S., & Shimizu, K. (2017). On circular correlation for data on the torus. Statistical Papers, 1-21.

Examples

# generate data from vmsin model
set.seed(1)
dat <- rvmsin(100, 2,3,-0.8,0,0)

# now calculate circular correlation(s) between the 2 columns of dat
circ_cor(dat, type="js")
circ_cor(dat, type="fl")
circ_cor(dat, type="tau1")
circ_cor(dat, type="tau2")

Analytic circular variances and correlations for bivariate angular models

Description

Analytic circular variances and correlations for bivariate angular models

Usage

circ_varcor_model(
  model = "vmsin",
  kappa1 = 1,
  kappa2 = 1,
  kappa3 = 0,
  mu1 = 0,
  mu2 = 0,
  nsim = 10000,
  ...
)

Arguments

Details

The function computes the analytic circular variances and correlations (both Jammalamadaka-Sarma (JS) and Fisher-Lee (FL) forms) for von Mises sine, von Mises cosine and bivariate wrapped normal distributions.

For wnorm2, expressions for the circular variances, JS and FL correlation coefficients can be found in Mardia and Jupp (2009), Jammalamadaka and Sarma (1988) and Fisher and Lee (1983) respectively. For vmsin and vmcos these expressions are provided in Chakraborty and Wong (2018).

Because the analytic expressions in vmsin and vmcos models involve infinite sums of product of Bessel functions, if any of kappa1, kappa2 and abs(kappa3) is larger than or equal to 150, IID Monte Carlo with sample size nsim is used to approximate rho_js for numerical stability. From rho_js, rho_fl is computed using Corollary 2.2 in Chakraborty and Wong (2018), which makes cost-complexity for the rho_fl evaluation to be of order O(nsim) for vmsin and vmcos models. (In general, rho_fl evaluation is of order O(nsim^2)).

In addition, for the vmcos model, when -150 < kappa3 < -1 or 50 < max(kappa1, kappa2, abs(kappa3)) <= 150, the analytic formulas in Chakraborty and Wong (2018) are used; however, the reciprocal of the normalizing constant and its partial derivatives are all calculated numerically via (quasi) Monte carlo method for numerical stability. These (quasi) random numbers can be provided through the argument qrnd, which must be a two column matrix, with each element being a (quasi) random number between 0 and 1. Alternatively, if n_qrnd is provided (and qrnd is missing), a two dimensional sobol sequence of size n_qrnd is generated via the function sobol from the R package qrng. If none of qrnd or n_qrnd is available, a two dimensional sobol sequence of size 1e4 is used.

Value

Returns a list with elements var1, var2 (circular variances for the first and second coordinates), rho_fl and rho_js (circular correlations). See details.

References

Fisher, N. I. and Lee, A. (1983). A correlation coefficient for circular data. Biometrika, 70(2):327-332.

Jammalamadaka, S. R. and Sarma, Y. (1988). A correlation coefficient for angular variables. Statistical theory and data analysis II, pages 349-364.

Mardia, K. and Jupp, P. (2009). Directional Statistics. Wiley Series in Probability and Statistics. Wiley.

Chakraborty, S. and Wong, S, W.K. (2018). On the circular correlation coefficients for bivariate von Mises distributions on a torus. arXiv e-print.

Examples

circ_varcor_model("vmsin", kappa1= 1, kappa2 = 2, kappa3 = 3)

# Monte Carlo approximation
set.seed(1)
dat <- rvmsin(1000, 1, 2, 3)
# sample circular variance
circ_var <- function(x)
  1 - mean(cos(x - atan2(mean(sin(x)), mean(cos(x))) ))
circ_var(dat[, 1])
circ_var(dat[, 2])
circ_cor(dat, "fl")
circ_cor(dat, "js")

Contour plot for angmcmc objects with bivariate data

Description

Contour plot for angmcmc objects with bivariate data

Usage

## S3 method for class 'angmcmc'
contour(
  x,
  fn = "MAP",
  type = "point-est",
  show.data = TRUE,
  xpoints = seq(0, 2 * pi, length.out = 100),
  ypoints = seq(0, 2 * pi, length.out = 100),
  levels,
  nlevels = 20,
  cex = 1,
  col = "red",
  alpha = 0.4,
  pch = 19,
  ...
)

Arguments

Details

contour.angmcmc is an S3 function for angmcmc objects that calls contour from graphics.

To estimate the mixture density required to construct the contour plot, first the parameter vector \eta is estimated by applying fn on the MCMC samples, yielding the (consistent) Bayes estimate \hat{\eta}. Then the mixture density f(x|\eta) at any point x is (consistently) estimated by f(x|\hat{\eta}).

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# now create a contour plot
contour(fit.vmsin.20)

Contourplot for bivariate angular mixture model densities

Description

Contourplot for bivariate angular mixture model densities

Usage

contour_model(
  model = "vmsin",
  kappa1,
  kappa2,
  kappa3,
  mu1,
  mu2,
  pmix = rep(1/length(kappa1), length(kappa1)),
  xpoints = seq(0, 2 * pi, length.out = 100),
  ypoints = seq(0, 2 * pi, length.out = 100),
  levels,
  nlevels = 20,
  xlab = "x",
  ylab = "y",
  col = "black",
  lty = 1,
  main,
  ...
)

Arguments

Examples

contour_model('vmsin', 1, 1, 1.5, pi, pi)
contour_model('vmcos', 1, 1, 1.5, pi, pi)

Density and random deviates from an angmcmc object

Description

Density and random deviates from an angmcmc object

Usage

d_fitted(x, object, type = "point-est", fn = mean, log = FALSE, chain.no, ...)

r_fitted(n = 1, object, type = "point-est", fn = mean, chain.no, ...)

Arguments

Details

If type = 'point-est', density is evaluated/random samples are generated at a point estimate of the parameter values. To estimate the mixture density, first the parameter vector \eta is estimated by applying fn on the MCMC samples (using the function pointest), yielding the (consistent) Bayes estimate \hat{\eta}. Then the mixture density f(x|\eta) at any point x is (consistently) estimated by f(x|\hat{\eta}). The random deviates are generated from the estimated mixture density f(x|\hat{\eta}).

If type == 'post-pred', posterior predictive samples and densities are returned. That is, the average density S^{-1} \sum_{s = 1}^S f(x | \eta_s) is returned in d_fitted, where \eta_1, \dots, \eta_S is the set posterior MCMC samples obtained from object. In r_fitted, first a random sub-sample \eta_{(1)}, \dots, \eta_{(n)} of size n from the set of posterior samples \eta_1, \dots, \eta_S is drawn (with replacement if n > S). Then the i-th posterior predictive data point is generated from the mixture density f(x|\eta_{(i)}) for i = 1,..., n.

Value

d_fitted gives a vector the densities computed at the given points and r_fitted creates a vector (if univariate) or a matrix (if bivariate) with each row being a 2-D point, of random deviates.

Examples

set.seed(1)
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
d_fitted(c(0,0), fit.vmsin.20, type = "post-pred")
d_fitted(c(0,0), fit.vmsin.20, type = "point-est")

r_fitted(10, fit.vmsin.20, type = "post-pred")
r_fitted(10, fit.vmsin.20, type = "point-est")

Density plots for angmcmc objects

Description

Plot fitted angular mixture model density surfaces or curves.

Usage

## S3 method for class 'angmcmc'
densityplot(
  x,
  data = NULL,
  fn = mean,
  type = "point-est",
  log.density = FALSE,
  xpoints = seq(0, 2 * pi, length.out = 35),
  ypoints = seq(0, 2 * pi, length.out = 35),
  plot = TRUE,
  show.hist = ifelse(log.density, FALSE, TRUE),
  xlab,
  ylab,
  zlab = ifelse(log.density, "Log Density", "Density"),
  main,
  ...
)

Arguments

Details

When plot==TRUE, densityplot.angmcmc calls lattice::wireframe or plot from graphics to draw the surface or curve.

To estimate the mixture density, first the parameter vector \eta is estimated by applying fn on the MCMC samples, yielding the (consistent) Bayes estimate \hat{\eta}. Then the mixture density f(x|\eta) at any point x is (consistently) estimated by f(x|\hat{\eta}).

Note that densityplot.angmcmc does not plot the kernel densitie estimates of the MCMC parameters. (These plots can be obtained by first converting an angmcmc object to an mcmc object via as.mcmc.list, and then by using densplot from package coda on the resulting mcmc.list object. Instead, densityplot.angmcmc returns the surface (if 2-D) or the curve (if 1-D) of the fitted model density evaluated at the estimated parameter vector (obtain through pointest).

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# now create density surface with the default first 1/3 as burn-in and thin = 1
library(lattice)
densityplot(fit.vmsin.20)
# the viewing angles can be changed through the argument 'screen'
# (passed to lattice::wireframe)
densityplot(fit.vmsin.20, screen = list(z=-30, x=-60))
densityplot(fit.vmsin.20, screen = list(z=30, x=-60))
# the colors can be changed through 'col.regions'
cols <- grDevices::colorRampPalette(c("blue", "green",
                                      "yellow", "orange", "red"))(100)
densityplot(fit.vmsin.20, col.regions = cols)

# Now fit a vm mixture model
# illustration only - more iterations needed for convergence
fit.vm.20 <- fit_vmmix(wind$angle, ncomp = 3, n.iter =  20,
                             n.chains = 1)
densityplot(fit.vm.20)

Extract MCMC samples for parameters from an angmcmc object

Description

Extract MCMC samples for parameters from an angmcmc object

Usage

extractsamples(object, par.name, comp.label, chain.no, drop = TRUE, ...)

Arguments

Details

The default for both par.name and comp.label are the all possible choices available in object.

Value

Returns a four dimensional array with

dimension 1 - model parameters and mixing proportions dimention 2 - components dimension 3 - MCMC iterations dimension 4 - chain number

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# extract Markov chain realizations for kappa1 from component 1
extr_kappa1_1 <- extractsamples(fit.vmsin.20, "kappa1", 1)
# for kappa1 from component from all components
extr_kappa1 <- extractsamples(fit.vmsin.20, "kappa1")
# for all parameters in component 1
extr_1 <- extractsamples(fit.vmsin.20, comp.label = 1)

Fitting Bivariate and univariate angular mixture models

Description

Fitting Bivariate and univariate angular mixture models

Usage

fit_angmix(
  model = "vmsin",
  data,
  ncomp,
  cov.restrict = "NONE",
  unimodal.component = FALSE,
  start_par = NULL,
  rand_start = rep(FALSE, n.chains),
  method = "hmc",
  perm_sampling = FALSE,
  n.chains = 3,
  chains_parallel = TRUE,
  return_llik_contri = FALSE,
  int.displ = 3,
  epsilon = 0.1,
  L = 10,
  epsilon.random = TRUE,
  L.random = FALSE,
  burnin.prop = 0.5,
  tune.prop = 1,
  thin = 1,
  propscale = 0.05,
  n.iter = 500,
  pmix.alpha = NULL,
  norm.var = 1000,
  autotune = TRUE,
  show.progress = TRUE,
  accpt.prob.upper,
  accpt.prob.lower,
  epsilon.incr = 0.05,
  L.incr = 0.075,
  tune.incr = 0.05,
  tune_ave_size = 100,
  kappa_upper = 150,
  kappa_lower = 1e-04,
  return_tune_param = FALSE,
  qrnd = NULL,
  n_qrnd = NULL,
  ...
)

Arguments

Note

Sampling is done in log scale for the concentration parameters (kappa, kappa1 and kappa2).

Parallelization is done by default when more than one chain is used, but the chains can be run sequentially as well by setting chains_parallel = FALSE. To retain reproducibility while running multiple chains in parallel, the same RNG state is passed at the beginning of each chain. This is done by specifying future.seed = TRUE in future.apply::future_lapply call. Then at the beginning of the i-th chain, before drawing any parameters, i-many Uniform(0, 1) random numbers are generated using runif(i) (and then thrown away). This ensures that the RNG states across chains prior to random generation of the parameters are different, and hence, no two chains can become identical, even if they have the same starting and tuning parameters. This, however creates a difference between a fit_angmix call with multiple chains which is run sequentially by setting chains_parallel = FALSE, and another which is run sequentially because of a sequential plan() (or no plan()), with chains_parallel = TRUE. In the former, different RNG states are passed at the initiation of each chain.

References

Fruhwirth-Schnatter, S. (2011). Label switching under model uncertainty. Mixtures: Estimation and Application, 213-239.

Fruhwirth-Schnatter, S. (2001). Markov chain Monte Carlo estimation of classical and dynamic switching and mixture models. Journal of the American Statistical Association, 96(453), 194-209.

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_angmix("vmsin", tim8,
  ncomp = 3, n.iter = 20,
  n.chains = 1
)
fit.vmsin.20


# Parallelization is implemented via future_lapply from the
# package future.apply. To parallelize, first provide a parallel
# plan(); otherwise the chains will run sequentially.
# Note that not all plan() might work on every OS, as they execute
# functions defined internally in fit_mixmodel. We suggest
# plan(multisession) which works on every OS.

library(future)
library(parallel)
# plan(multisession, gc = TRUE) # parallelize chains

set.seed(1)
MC.fit <- fit_angmix("vmsin", tim8,
  ncomp = 3, n.iter = 5000,
  n.chains = 3
)


pointest(MC.fit)

MC.fix <- fix_label(MC.fit)

contour(MC.fit)
contour(MC.fix)
lpdtrace(MC.fit)

Stepwise fitting of angular mixture models with incremental component sizes and optimum model selection

Description

Stepwise fitting of angular mixture models with incremental component sizes and optimum model selection

Usage

fit_incremental_angmix(
  model,
  data,
  crit = "LOOIC",
  start_ncomp = 1,
  max_ncomp = 10,
  L = NULL,
  fn = mean,
  fix_label = NULL,
  form = 2,
  start_par = NULL,
  prev_par = TRUE,
  logml_maxiter = 10000,
  return_all = FALSE,
  save_fits = FALSE,
  save_file = NULL,
  save_dir = "",
  silent = FALSE,
  return_llik_contri = (crit %in% c("LOOIC", "WAIC")),
  use_best_chain = TRUE,
  alpha = 0.05,
  bonferroni_alpha = TRUE,
  bonferroni_adj_type = "decreasing",
  ...
)

Arguments

Details

The goal is to fit an angular mixture model with an optimally chosen component size K. To obtain an optimum K, mixture models with incremental component sizes between start_ncomp and max_ncomp are fitted incrementally using fit_angmix, starting from K = 1. If the model selection criterion crit is "LOOIC" or "WAIC", then a test of hypothesis H_{0K}: expected log predictive density (elpd) for the fitted model with K components >= elpd for the fitted model with K + 1 components, is performed at every K >= 1. The test-statistic used for the test is an approximate z-score based on the normalized estimated elpd difference between the two models obtained from compare, which provides estimated elpd difference along with its standard error estimate. Because the computed standard error of elpd difference can be overly optimistic when the elpd difference is small (in particular < 4), a conservative worst-case estimate (equal to twice of the computed standard error) is used in such cases. To account for multiplicity among the M = (max_ncomp - start_ncomp) possible sequential tests performed, by default a Bonferroni adjustment to the test level alpha is made. Set ⁠bonferroni_alpha = FALSE} to remove the adjustment. To encourage parsimony in the final model, by default (⁠bonferroni_adj_type = "decreasing"⁠) a decreasing sequence of adjusted alphas of the form ⁠alpha * (0.5)^(1:M) / sum((0.5)^(1:M))⁠is used. Set⁠bonferroni_adj_type = "equal"⁠to use equal sequence of adjusted alphas (i.e.,⁠alpha/M') instead.

The incremental fitting stops if H_{0K} cannot be rejected (at level alpha) for some K >= 1; this K is then regarded as the optimum number of components. If crit is not "LOOIC" or "WAIC" then mixture model with the first minimum value of the model selection criterion crit is taken as the best model.

Note that in each intermediate fitted model, the total number of components (instead of the number of "non-empty components") in the model is used to estimate of the true component size, and then the fitted model is penalized for model complexity (via the model selection criterion used). This approach of selecting an optimal K follows the perspective "let two component specific parameters be identical" for overfitting mixtures, and as such the Dirichlet prior hyper-parameters pmix.alpha (passed to fit_angmix) should be large. See Fruhwirth-Schnatter (2011) for more deltails.

Note that the stability of bridge_sampler used in marginal likelihood estimation heavily depends on stationarity of the chains. As such, while using this criterion, we recommending running the chain long enough, and setting fix_label = TRUE for optimal performance.

Value

Returns a named list (with class = stepfit) with the following seven elements:

fit.all (if return_all = TRUE) - a list all angmcmc objects created at each component size;

fit.best - angmcmc object corresponding to the optimum component size;

ncomp.best - optimum component size (integer);

crit - which model comparison criterion used (one of ⁠"LOOIC", "WAIC", "AIC", "BIC", "DIC"⁠ or "LOGML");

crit.all - all crit values calculated (for all component sizes);

crit.best - crit value for the optimum component size; and

maxllik.all - maximum (obtained from MCMC iterations) log likelihood for all fitted models

maxllik.best - maximum log likelihodd for the optimal model; and

check_min - logical; is the optimum component size less than max_ncomp?

References

Fruhwirth-Schnatter, S.: Label switching under model uncertainty. In: Mengerson, K., Robert, C., Titterington, D. (eds.) Mixtures: Estimation and Application, pp. 213-239. Wiley, New York (2011).

Examples

# illustration only - more iterations needed for convergence
set.seed(1)
fit.vmsin.step.15 <- fit_incremental_angmix("vmsin", tim8, "BIC", start_ncomp = 1,
                                          max_ncomp = 3, n.iter = 15,
                                          n.chains = 1, save_fits=FALSE)
(fit.vmsin.best.15 <- bestmodel(fit.vmsin.step.15))
lattice::densityplot(fit.vmsin.best.15)

Fitting bivariate von Mises cosine model mixtures using MCMC

Description

Fitting bivariate von Mises cosine model mixtures using MCMC

Usage

fit_vmcosmix(...)

Arguments

Details

Wrapper for fit_angmix with model = "vmcos".

Examples

# illustration only - more iterations needed for convergence
fit.vmcos.10 <- fit_vmcosmix(tim8, ncomp = 3, n.iter =  10,
                             n.chains = 1)
fit.vmcos.10

Fitting univariate von Mises mixtures using MCMC

Description

Fitting univariate von Mises mixtures using MCMC

Usage

fit_vmmix(...)

Arguments

Details

Wrapper for fit_angmix with 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

Fitting bivariate von Mises sine model mixtures using MCMC

Description

Fitting bivariate von Mises sine model mixtures using MCMC

Usage

fit_vmsinmix(...)

Arguments

Details

Wrapper for fit_angmix with 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

Fitting bivariate wrapped normal model mixtures using MCMC

Description

Fitting bivariate wrapped normal model mixtures using MCMC

Usage

fit_wnorm2mix(...)

Arguments

Details

Wrapper for fit_angmix with model = "wnorm2".

Examples

# illustration only - more iterations needed for convergence
fit.wnorm2.10 <- fit_wnorm2mix(tim8, ncomp = 3, n.iter =  10,
                             n.chains = 1)
fit.wnorm2.10

Fitting univariate wrapped normal mixtures using MCMC

Description

Fitting univariate wrapped normal mixtures using MCMC

Usage

fit_wnormmix(...)

Arguments

Details

Wrapper for fit_angmix with 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

Fix label switching in angmcmc objects

Description

Fix label switching in angmcmc objects

Usage

fix_label(object, ...)

Arguments

Details

fix_label is a wrapper for label.switching from package label.switching for angmcmc objects. The arguments z, K, complete, mcmc, p and data are appropriately filled in from object. The label.switching argument method can be a scalar or vector; for this wrapper it defaults to "STEPHENS" if the angmcmc was created with permutation sampling (by setting perm_sampling = TRUE in fit_angmix), and to "DATA-BASED" otherwise.

Value

Returns a single angmcmc object or a list of angmcmc objects (according as whether the argument method is a scalar or vector) with label switchings corrected (after burn-in and thin) according to the resulting permutation from label.switching.

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# now apply fix_label
fit.vmsin.20.fix <- fix_label(fit.vmsin.20)

Angular MCMC (angmcmc) Object

Description

Checking for and creating an angmcmc object

Usage

is.angmcmc(object)

angmcmc(...)

Arguments

Details

angmcmc objects are classified lists that are created when any of the five mixture model fitting functions, viz., fit_vmmix, fit_wnormmix, fit_vmsinmix, fit_vmcosmix and fit_wnorm2mix is used. An angmcmc object contains a number of elements, including the dataset, the model being fitted on the dataset and dimension of the model (univariate or bivariate), the tuning parameters used, MCMC samples for the mixture model parameters, the (hidden) component or cluster indicators for data points in each iteration and the (iteration-wise) log likelihood and log posterior density values (both calculated upto some normalizing constants). When printed, an angmcmc object returns a brief summary of the function arguments used to produce the object and the average acceptance rate of the proposals (in HMC and RWMH) used over iterations. An angmcmc object can be used as an argument for the diagnostic and post-processing functions available in BAMBI for making further inferences.

Value

logical. Is the input an angmcmc object?

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
is.angmcmc(fit.vmsin.20)

Finding latent allocation (component indicators) from an angmcmc object

Description

Finding latent allocation (component indicators) from an angmcmc object

Usage

latent_allocation(object, ...)

Arguments

Details

In order to find the latent component indicators, estimates of mixing proportions and model parameters are first computed via pointest. Then, a data point is assigned label j, if the j-th component gives highest density for that point.

Value

Returns a vector of length n, where n is the length (if univariate) or number of rows (if bivariate) of the data used in original fit. i-th entry of the output vector provides component label for the i-th data point.

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# now find latent allocation
latent_allocation(fit.vmsin.20)

Extract Log-Likelihood from angmcmc objects

Description

Extract Log-Likelihood from angmcmc objects

Usage

## S3 method for class 'angmcmc'
logLik(object, method = 1, fn, ...)

Arguments

Details

There are two ways to estimate the log likelihood from the model. If method = 1, then log likelihood is estimated by applying fn (defaults to max, if method = 1) direclty on the log likelihood values from observed during the MCMC run. On the other hand, if method == 2, then parameter estimates are first computed using pointest with fn (defaults to "MODE", if method == 2) applied on the MCMC samples, and then then log likelihood is evaluated at the parameter estimates.

The degrees of the likelihood function is the total number of free parameters estimated in the mixture models, which is equal to 6K - 1 for bivariate models (vmsin, vmcos and wnorm2), or 3K - 1 for univariate models (vm and wnorm), where K denotes the number of components in the mixture model.

Value

Returns an object of class logLik. This is a number (the estimated log likelihood) with attributes "df" (degrees of freedom) and "nobs" (number of observations).

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
logLik(fit.vmsin.20)

Leave-one-out cross-validation (LOO) for angmcmc objects

Description

Leave-one-out cross-validation (LOO) for angmcmc objects

Usage

## S3 method for class 'angmcmc'
loo(x, ...)

Arguments

Details

Note that loo.angmcmc calls loo for computation. If the likelihood contribution of each data point for each MCMC iteration is available in object (can be returned by setting return_llik_contri = TRUE) during fit_angmix call), loo.array is used; otherwise loo.function is called. Computation is much faster if the likelihood contributions are available - however, they are very memory intensive, and by default not returned in fit_angmix.

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1, return_llik_contri = TRUE)
library(loo)
loo(fit.vmsin.20)

Trace and autocorrelation plots of log posterior density or log likelihood from an angmcmc object

Description

Trace and autocorrelation plots of log posterior density or log likelihood from an angmcmc object

Usage

lpdtrace(
  object,
  chain.no,
  use.llik = FALSE,
  plot.autocor = FALSE,
  lag.max = NULL,
  ...
)

Arguments

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# log posterior density trace
lpdtrace(fit.vmsin.20)
# log likelihood trace
lpdtrace(fit.vmsin.20, use.llik = TRUE)

Trace plot for parameters from an angmcmc object

Description

Trace plot for parameters from an angmcmc object

Usage

paramtrace(object, par.name, comp.label, chain.no, ...)

Arguments

Value

Returns a single plot if a single par and a single comp.label is supplied. Otherwise, a series of plots is produced.

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# trace plot for kappa1 in component 1
paramtrace(fit.vmsin.20, "kappa1", 1)
# for kappa1 in all components
paramtrace(fit.vmsin.20, "kappa1")
# for all parameters in component 1
paramtrace(fit.vmsin.20, comp.label = 1)

Summary plots for angmcmc objects

Description

Summary plots for angmcmc objects

Usage

## S3 method for class 'angmcmc'
plot(
  x,
  par.name,
  comp.label,
  chain.no,
  do.paramtrace = TRUE,
  do.lpdtrace = TRUE,
  use.llik = FALSE,
  ...
)

Arguments

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
plot(fit.vmsin.20)

Point estimates for parameters from an angmcmc object

Description

Point estimates for parameters from an angmcmc object

Usage

pointest(object, fn = mean, par.name, comp.label, chain.no, ...)

Arguments

Value

Returns a matrix of point estimates, or vector of point estimates if length(par.name)==1 or length(comp.label)==1.

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# estimate parameters by sample mean
(est_mean <- pointest(fit.vmsin.20))
# estimate parameters by sample median
(est_median <- pointest(fit.vmsin.20, fn = median))
# estimate parameters by MAP
(est_median <- pointest(fit.vmsin.20, fn = "MODE"))

Quantile estimates for parameters from an angmcmc object

Description

Quantile estimates for parameters from an angmcmc object

Usage

## S3 method for class 'angmcmc'
quantile(x, par.name, comp.label, chain.no, probs = seq(0, 1, 0.25), ...)

Arguments

Value

Returns a three dimensional array of quantiles, or a matrix (vector) of quantiles if one (or two) among par.name, comp.label, probs has length 1.

Examples

# first fit a vmsin mixture model
# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
# 0.025th quantiles
(quant_025 <- quantile(fit.vmsin.20, prob = 0.025))
# 0.975th quantiles
(quant_975 <- quantile(fit.vmsin.20, prob = 0.975))
# default quantiles
(quant_def <- quantile(fit.vmsin.20))

The univariate von Mises distribution

Description

The univariate von Mises distribution

Usage

rvm(n, kappa = 1, mu = 0)

dvm(x, kappa = 1, mu = 0, log = FALSE)

Arguments

Details

If mu and kappa are not specified they assume the default values of 0 and 1 respectively.

The univariate von Mises distribution has density

f(x) = 1/(2\pi I_0 (\kappa)) \exp(\kappa \cos(x - mu))

where I_0 (\kappa) denotes the modified Bessel function of the first kind with order 0 evaluated at the point \kappa.

Value

dvm gives the density and 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)

The bivariate von Mises cosine model

Description

The bivariate von Mises cosine model

Usage

rvmcos(
  n,
  kappa1 = 1,
  kappa2 = 1,
  kappa3 = 0,
  mu1 = 0,
  mu2 = 0,
  method = "naive"
)

dvmcos(
  x,
  kappa1 = 1,
  kappa2 = 1,
  kappa3 = 0,
  mu1 = 0,
  mu2 = 0,
  log = FALSE,
  ...
)

Arguments

Details

The bivariate von Mises cosine model density at the point x = (x_1, x_2) is given by

f(x) = C_c (\kappa_1, \kappa_2, \kappa_3) \exp(\kappa_1 \cos(T_1) + \kappa_2 \cos(T_2) + \kappa_3 \cos(T_1 - T_2))

where

T_1 = x_1 - \mu_1; T_2 = x_2 - \mu_2

and C_c (\kappa_1, \kappa_2, \kappa_3) denotes the normalizing constant for the cosine model.

Because C_c involves an infinite alternating series with product of Bessel functions, if kappa3 < -5 or max(kappa1, kappa2, abs(kappa3)) > 50, C_c is evaluated numerically via (quasi) Monte carlo method for numerical stability. These (quasi) random numbers can be provided through the argument qrnd, which must be a two column matrix, with each element being a (quasi) random number between 0 and 1. Alternatively, if n_qrnd is provided (and qrnd is missing), a two dimensional sobol sequence of size n_qrnd is generated via the function sobol from the R package qrng. If none of qrnd or n_qrnd is available, a two dimensional sobol sequence of size 1e4 is used. By default Monte Carlo approximation is used only if kappa3 < -5 or max(kappa1, kappa2, abs(kappa3)) > 50. However, a forced Monte Carlo approximation can be made (irrespective of the choice of kappa1, kappa2 and kappa3) by setting force_approx_const = TRUE. See examples.

Value

dvmcos gives the density and rvmcos 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,
# dvmcos returns single density
dvmcos(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
dvmcos(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
# dvmcos returns a vector with ith element being the density
# evaluated at x with parameter values kappa1[i], kappa2[i],
# kappa3[i], mu1[i] and mu2[i]
dvmcos(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 dvmcos returns a vector 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]
dvmcos(x, kappa1, kappa2, kappa3, mu1, mu2)

# when parameters are all scalars, number of observations generated
# by rvmcos is n
rvmcos(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 rvmcos is the same as the length of the
# recycled vectors
rvmcos(n, kappa1, kappa2, kappa3, mu1, mu2)



## Visualizing (quasi) Monte Carlo based approximations of
## the normalizing constant through density evaluations.

# "good" setup, where the analytic formula for C_c can be
# calculated without numerical issues
# kappa1 = 1, kappa2 = 1, kappa3 = -2, mu1 = pi, mu2 = pi

n_qrnd <-  (1:500)*20
# analytic
good.a <- dvmcos(c(3,3), 1, 1, -2, pi, pi, log=TRUE)
# using quasi Monte Carlo
good.q <- sapply(n_qrnd,
                 function(j)
                   dvmcos(c(3,3), 1, 1, -2, pi, pi,
                          log=TRUE, n_qrnd = j,
                          force_approx_const = TRUE))
# using ordinary Monte Carlo
set.seed(1)
good.r <- sapply(n_qrnd,
                 function(j)
                   dvmcos(c(3,3), 1, 1, -2, pi, pi,
                          log=TRUE,
                          qrnd = matrix(runif(2*j), ncol = 2),
                          force_approx_const = TRUE))


plot(n_qrnd, good.q, ylim = range(good.a, good.q, good.r),
     col = "orange", type = "l",
     ylab = "",
     main = "dvmcos(c(3,3), 1, 1, -2, pi, pi, log = TRUE)")
points(n_qrnd, good.r, col = "skyblue", type = "l")
abline(h = good.a, lty = 2, col = "grey")
legend("topright",
       legend = c("Sobol", "Random", "Analytic"),
       col = c("orange", "skyblue", "grey"),
       lty = c(1, 1, 2))


# "bad" setup, where the calculating C_c
# numerically using the analytic formula is problematic
# kappa1 = 100, kappa2 = 100, kappa3 = -200, mu1 = pi, mu2 = pi

n_qrnd <-  (1:500)*20

# using quasi Monte Carlo
bad.q <- sapply(n_qrnd,
                function(j)
                  dvmcos(c(3,3), 100, 100, -200, pi, pi,
                         log=TRUE, n_qrnd = j,
                         force_approx_const = TRUE))
# using ordinary Monte Carlo
set.seed(1)
bad.r <- sapply(n_qrnd,
                function(j)
                  dvmcos(c(3,3), 100, 100, -200, pi, pi,
                         log=TRUE,
                         qrnd = matrix(runif(2*j), ncol = 2),
                         force_approx_const = TRUE))


plot(n_qrnd, bad.q, ylim = range(bad.q, bad.r),
     col = "orange", type = "l",
     ylab = "",
     main = "dvmcos(c(3,3), 100, 100, -200, pi, pi, log = TRUE)")
points(n_qrnd, bad.r, col = "skyblue", type = "l")
legend("topright",
       legend = c("Sobol", "Random"),
       col = c("orange", "skyblue"), lty = 1)

The bivariate von Mises cosine model mixtures

Description

The bivariate von Mises cosine model mixtures

Usage

rvmcosmix(n, kappa1, kappa2, kappa3, mu1, mu2, pmix, method = "naive", ...)

dvmcosmix(x, kappa1, kappa2, kappa3, mu1, mu2, pmix, log = FALSE, ...)

Arguments

Details

All the argument vectors pmix, kappa1, kappa2, kappa3, mu1 and mu2 must be of the same length ( = component size of the mixture model), with j-th element corresponding to the j-th component of the mixture distribution.

The bivariate von Mises cosine model mixture distribution with component size K = length(pmix) has density

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 j; p[j]; \kappa_1[j], \kappa_2[j], \kappa_3[j]; and \mu_1[j], \mu_2[j] respectively denote the mixing proportion, the three concentration parameters and the two mean parameter for the j-th cluster, j = 1, ..., K, and f(. ; \kappa_1, \kappa_2, \kappa_3, \mu_1, \mu_2) denotes the density function of the von Mises cosine model with concentration parameters \kappa_1, \kappa_2, \kappa_3 and mean parameters \mu_1, \mu_2.

Value

dvmcosmix computes the density and rvmcosmix 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
dvmcosmix(x, kappa1, kappa2, kappa3, mu1, mu2, pmix)

# number of observations generated from the mixture distribution is n
rvmcosmix(n, kappa1, kappa2, kappa3, mu1, mu2, pmix)

The univariate von Mises mixtures

Description

The univariate von Mises mixtures

Usage

rvmmix(n, kappa, mu, pmix)

dvmmix(x, kappa, mu, pmix, log = FALSE)

Arguments

Details

pmix, mu and kappa must be of the same length, with j-th element corresponding to the j-th component of the mixture distribution.

The univariate von Mises mixture distribution with component size K = length(pmix) has density

g(x) = p[1] * f(x; \kappa[1], \mu[1]) + ... + p[K] * f(x; \kappa[K], \mu[K])

where p[j], \kappa[j], \mu[j] respectively denote the mixing proportion, concentration parameter and the mean parameter for the j-th component and f(. ; \kappa, \mu) denotes the density function of the (univariate) von Mises distribution with mean parameter \mu and concentration parameter \kappa.

Value

dvmmix computes the density and 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)

The bivariate von Mises sine model

Description

The bivariate von Mises sine model

Usage

rvmsin(
  n,
  kappa1 = 1,
  kappa2 = 1,
  kappa3 = 0,
  mu1 = 0,
  mu2 = 0,
  method = "naive"
)

dvmsin(x, kappa1 = 1, kappa2 = 1, kappa3 = 0, mu1 = 0, mu2 = 0, log = FALSE)

Arguments

Details

The bivariate von Mises sine model density at the point x = (x_1, x_2) is given by

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

T_1 = x_1 - \mu_1; T_2 = x_2 - \mu_2

and 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 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 method = "naive" (default) a (naive) bivariate rejection sampling scheme with (bivariate) uniform propsoal is used.

Note that although method = "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 n is not large.

Value

dvmsin gives the density and 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)

The bivariate von Mises sine model mixtures

Description

The bivariate von Mises sine model mixtures

Usage

rvmsinmix(n, kappa1, kappa2, kappa3, mu1, mu2, pmix, method = "naive")

dvmsinmix(x, kappa1, kappa2, kappa3, mu1, mu2, pmix, log = FALSE)

Arguments

Details

All the argument vectors pmix, kappa1, kappa2, kappa3, mu1 and mu2 must be of the same length ( = component size of the mixture model), with j-th element corresponding to the j-th component of the mixture distribution.

The bivariate von Mises sine model mixture distribution with component size K = length(p.mix) has density

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 j; p[j]; \kappa_1[j], \kappa_2[j], \kappa_3[j]; and \mu_1[j], \mu_2[j] respectively denote the mixing proportion, the three concentration parameters and the two mean parameter for the j-th component, j = 1, ..., K, and f(. ; \kappa_1, \kappa_2, \kappa_3, \mu_1, \mu_2) denotes the density function of the von Mises sine model with concentration parameters \kappa_1, \kappa_2, \kappa_3 and mean parameters \mu_1, \mu_2.

Value

dvmsinmix computes the density (vector if x is a two column matrix with more than one row) and 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)

The univariate Wrapped Normal distribution

Description

The univariate Wrapped Normal distribution

Usage

rwnorm(n = 1, kappa = 1, mu = 0)

dwnorm(x, kappa = 1, mu = 0, int.displ, log = FALSE)

Arguments

Details

If mu and kappa are not specified they assume the default values of 0 and 1 respectively.

The univariate wrapped normal distribution has density

f(x) = \sqrt(\kappa/(2\pi)) \sum \exp(-\kappa/2 (x - \mu(2\pi\omega))^2)

where the sum extends over all integers \omega, and is approximated by a sum over \omega in \{-M, -M+1, ..., M-1, M \} if int.displ = M.

Value

dwnorm gives the density and 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)

The bivariate Wrapped Normal distribution

Description

The bivariate Wrapped Normal distribution

Usage

rwnorm2(n, kappa1 = 1, kappa2 = 1, kappa3 = 0, mu1 = 0, mu2 = 0, ...)

dwnorm2(
  x,
  kappa1 = 1,
  kappa2 = 1,
  kappa3 = 0,
  mu1 = 0,
  mu2 = 0,
  int.displ,
  log = FALSE
)

Arguments

Details

The bivariate wrapped normal density at the point x = (x_1, x_2) is given by,

f(x) = \sqrt((\kappa_1 \kappa_2 - (\kappa_3)^2)) / (2\pi) \sum \exp(-1/2 * (\kappa_1 (T_1)^2 + \kappa_2 (T_2)^2 + 2 \kappa_3 (T_1) (T_2)) )

where

T_1 = T_1(x, \mu, \omega) = (x_1 - \mu_1(2\pi\omega_1))

T_2 = T_2(x, \mu, \omega) = (x_2 - \mu_1(2\pi\omega_2))

the sum extends over all pairs of integers \omega = (\omega_1, \omega_2), and is approximated by a sum over (\omega_1, \omega_2) in \{-M, -M+1, ..., M-1, M \}^2 if int.displ = M.

Note that above density is essentially the "wrapped" version of a bivariate normal density with mean

\mu = (\mu_1, \mu_2)

and dispersion matrix \Sigma = \Delta^{-1}, with \Delta being a 2 \times 2 symmetric, positive definite matrix with diagonal entries \kappa_1 and \kappa_2 and the off-diagonal entries \kappa_3.

Value

dwnorm2 gives the density and rwnorm2 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,
# dwnorm2 returns single density
dwnorm2(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
dwnorm2(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
# dwnorm2 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]
dwnorm2(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 dwnorm2 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]
dwnorm2(x, kappa1, kappa2, kappa3, mu1, mu2)

# when parameters are all scalars, number of observations generated
# by rwnorm2 is n
rwnorm2(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 rwnorm2 is the same as the length of the
# recycled vectors
rwnorm2(n, kappa1, kappa2, kappa3, mu1, mu2)

The bivariate Wrapped Normal mixtures

Description

The bivariate Wrapped Normal mixtures

Usage

rwnorm2mix(n, kappa1, kappa2, kappa3, mu1, mu2, pmix, ...)

dwnorm2mix(x, kappa1, kappa2, kappa3, mu1, mu2, pmix, int.displ, log = FALSE)

Arguments

Details

All the argument vectors pmix, kappa1, kappa2, kappa3, mu1 and mu2 must be of the same length, with j-th element corresponding to the j-th component of the mixture distribution.

The bivariate wrapped normal mixture distribution with component size K = length(pmix) has density

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 j; p[j]; \kappa_1[j], \kappa_2[j], \kappa_3[j]; and \mu_1[j], \mu_2[j] respectively denote the mixing proportion, the three concentration parameters and the two mean parameter for the j-th component, j = 1, ..., K, and f(. ; \kappa_1, \kappa_2, \kappa_3, \mu_1, \mu_2) denotes the density function of the wrapped normal distribution with concentration parameters \kappa_1, \kappa_2, \kappa_3 and mean parameters \mu_1, \mu_2.

Value

dwnorm2mix computes the density and rwnorm2mix 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
dwnorm2mix(x, kappa1, kappa2, kappa3, mu1, mu2, pmix)

# number of observations generated from the mixture distribution is n
rwnorm2mix(n, kappa1, kappa2, kappa3, mu1, mu2, pmix)

The univariate Wrapped Normal mixtures

Description

The univariate Wrapped Normal mixtures

Usage

rwnormmix(n = 1, kappa, mu, pmix)

dwnormmix(x, kappa, mu, pmix, int.displ = 3, log = FALSE)

Arguments

Details

pmix, mu and kappa must be of the same length, with j-th element corresponding to the j-th component of the mixture distribution.

The univariate wrapped normal mixture distribution with component size K = length(pmix) has density

g(x) = p[1] * f(x; \kappa[1], \mu[1]) + ... + p[K] * f(x; \kappa[K], \mu[K])

where p[j], \kappa[j], \mu[j] respectively denote the mixing proportion, concentration parameter and the mean parameter for the j-th component and f(. ; \kappa, \mu) denotes the density function of the (univariate) wrapped normal distribution with mean parameter \mu and concentration parameter \kappa.

Value

dwnormmix computes the density and 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)

Select chains from angmcmc objects

Description

Select chains from angmcmc objects

Usage

select_chains(object, chain.no, ...)

Arguments

Value

Returns another angmcmc object with only selected chains passed through chain.no

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             L = c(10, 12), chains_parallel = FALSE,
                             n.chains = 2)
fit.vmsin.20
fit.vmsin.20.1 <- select_chains(fit.vmsin.20, 1)
fit.vmsin.20.1

Summary statistics for parameters from an angmcmc object

Description

Summary statistics for parameters from an angmcmc object

Usage

## S3 method for class 'angmcmc'
summary(object, par.name, comp.label, chain.no, ...)

Arguments

Details

Computes (after thinning and discarding burn-in) point estimates with 95% posterior credible sets for all components and all parameters, together with the sample averages of log likelihood and log posterior density.

Value

Returns a list with elements estimate, lower, upper, llik and lpd. estimate is itself a list with three elements: mean, median and mode providing the sample mean, sample median and (sample) MAP estimates.

Note that summary.angmcmc has its own print method, providing a table the estimated mean and 95% credible intervals for each parameter

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1)
summary(fit.vmsin.20)

Surface for bivariate angular mixture model densities

Description

Surface for bivariate angular mixture model densities

Usage

surface_model(
  model = "vmsin",
  kappa1,
  kappa2,
  kappa3,
  mu1,
  mu2,
  pmix = rep(1/length(kappa1), length(kappa1)),
  xpoints = seq(0, 2 * pi, length.out = 30),
  ypoints = seq(0, 2 * pi, length.out = 30),
  log.density = FALSE,
  xlab = "x",
  ylab = "y",
  zlab = ifelse(log.density, "Log Density", "Density"),
  main,
  ...
)

Arguments

Examples

surface_model('vmsin', 1, 1, 1.5, pi, pi)
surface_model('vmcos', 1, 1, 1.5, pi, pi)

Backbone Dihedral Angles of Triose Phosphate Isomerase (8TIM)

Description

A dataset consisting of 490 pairs of backbone dihedral angles (in radian scale [0, 2\pi) ) (\phi, \psi) for the protein Triose Phosphate Isomerase (8TIM). The angles were obtained first by using the DSSP software on the PDB file for 8TIM to get the backbone angles (in degrees), and then by converting all angles into radians. Due to the presence of different secondary structures (helices, sheets and loops) in the protein, the angular data show considerable variability, and is multimodal with noticeably distinct clusters.

Usage

data(tim8)

Format

A data frame with 490 rows and 2 variables (backbone dihedral angles) phi and psi.

Source

8TIM PDB file: http://www.rcsb.org/pdb/explore.do?structureId=8tim.

DSSP software: https://swift.cmbi.umcn.nl/gv/dssp/.

Maximum likelihood estimation of bivariate von Mises parameters

Description

Maximum likelihood estimation of bivariate von Mises parameters

Usage

vm2_mle(data, model = c("vmsin", "vmcos"), ...)

Arguments

Details

The parameters kappa1 and kappa2 are optimized in log scales. The method of optimization used (passed to optim) can be specified through method in ... (defaults to "L-BFGS-B"). Note, however, that lower (0) and upper (2*pi) bounds for mu1 and mu2 are specified; so not all methods implemented in optim will work.

Value

An object of class mle-class.

Examples

pars <- list(kappa1 = 3, kappa2 = 2, kappa3 = 1.5, mu1 = 0.5, mu2 = 1.5)
nsamp <- 2000
model <- "vmsin"
set.seed(100)
dat_gen <- do.call(paste0("r", model), c(list(n = nsamp), pars))

est <- vm2_mle(dat_gen, model = model)
library(stats4)
coef(est)
vcov(est)

Watanabe-Akaike Information Criterion (WAIC) for angmcmc objects

Description

Watanabe-Akaike Information Criterion (WAIC) for angmcmc objects

Usage

## S3 method for class 'angmcmc'
waic(x, ...)

Arguments

Details

Given a deviance function D(\eta) = -2 \log(p(y|\eta)), and an estimate \eta* = (\sum \eta_i) / n of the posterior mean E(\eta|y), where y = (y_1, ..., y_n) denote the data, \eta is the unknown parameter vector of the model, \eta_1, ..., \eta_N are MCMC samples from the posterior distribution of \eta given y and p(y|\eta) is the likelihood function, the Watanabe-Akaike Information Criterion (WAIC) is defined as

WAIC = LPPD - p_W

where

LPPD = \sum_{i=1}^n \log (N^{-1} \sum_{s=1}^N p(y_i|\eta_s) )

and (form 1 of)

p_W = 2 \sum_{i=1}^n [ \log (N^{-1} \sum_{s=1}^N p(y_i|\eta_s) ) - N^{-1} \sum_{s=1}^N \log \:p(y_i|\eta_s) ].

An alternative form (form 2) for p_W is given by

p_W = \sum_{i=1}^n \hat{var} \log p(y_i|\eta)

where for i = 1, ..., n, \hat{var} \log p(y_i|\eta) denotes the estimated variance of \log p(y_i|\eta) based on the realizations \eta_1, ..., \eta_N.

Note that waic.angmcmc calls waic for computation. If the likelihood contribution of each data point for each MCMC iteration is available in object (can be returned by setting return_llik_contri = TRUE) during fit_angmix call), waic.array is used; otherwise waic.function is called. Computation is much faster if the likelihood contributions are available - however, they are very memory intensive, and by default not returned in fit_angmix.

Value

Computes the WAIC for a given angmcmc object.

Examples

# illustration only - more iterations needed for convergence
fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
                             n.chains = 1, return_llik_contri = TRUE)
library(loo)
waic(fit.vmsin.20)

Saturna Island wind directions

Description

A dataset consisting of 239 observations on wind direction in radians (original measurements were in 10s of degrees), measured at Saturna Island, British Columbia, Canada during October 1-10, 2016 (obtained from Environment Canada website). There was a severe storm during October 4-7, which caused significant fluctuations among the wind directions. As a result the angular data show a clear multimodality.

Usage

data(wind)

Format

A data frame with 239 rows and 2 columns; the column "angle" provides the angular direction (in radian) and the column day provides the days on which the data points were collected (ranges between 1-10, corresponding to October 1-10, 2016).

Source

Environment Canada: https://climate.weather.gc.ca/climate_data/data_quality_e.html.

CBC news on the storm: https://www.cbc.ca/news/canada/british-columbia/storm-bc-1.3795204.

Wrap angles into [-pi, pi] or [0, 2*pi]

Description

Wrap angles into [-pi, pi] or [0, 2*pi]

Usage

zero_to_2pi(x)

minuspi_to_pi(x)

Arguments

Details

minuspi_to_pi wraps x into [-pi, pi], while zero_to_pi wraps x into [0, 2*pi].

Examples

dat <- matrix(runif(100, -pi, pi), ncol=2)
dat1 <- zero_to_2pi(dat)
dat2 <- minuspi_to_pi(dat1)
all.equal(dat, dat2)