Package {nrMCmix}

Type:	Package
Title:	Non-Reversible Markov Chains for Bayesian Mixture Models
Version:	0.3
Date:	2026-10-13
Description:	Provides Bayesian inference for finite and infinite mixture models using Markov chain Monte Carlo (MCMC) methods. Supports Gaussian, multivariate Gaussian, Poisson, Binomial, and multivariate Bernoulli mixture kernels with different covariance structures and flexible prior specifications. Provides tools for data generation, posterior density evaluation, clustering and similarity analysis, and advanced MCMC sampling strategies, including reversible and non-reversible samplers: Ascolani et al. (2026) <doi:10.48550/arXiv.2510.03226>, informed proposals: Zanella (2020) <doi:10.1080/01621459.2019.1585255>, collapsed Gibbs sampling: Diebolt and Robert (1994) <doi:10.1111/j.2517-6161.1994.tb01985.x>, and split-and-merge samplers: Dahl and Newcomb (2021) <doi:10.1080/00949655.2021.1998502>.
License:	MIT + file LICENSE
Encoding:	UTF-8
Depends:	R (≥ 4.1.0)
Imports:	base, coda, graphics, grDevices, fields, ggplot2, salso, Rcpp (≥ 1.0.12), Rdpack, stats
RdMacros:	Rdpack
LinkingTo:	Rcpp, RcppArmadillo
Suggests:	knitr, rmarkdown, bmixture, datasets
RoxygenNote:	7.3.1
VignetteBuilder:	knitr
NeedsCompilation:	yes
Packaged:	2026-10-13 16:00:28 UTC; paolomanildo
Author:	Filippo Ascolani [aut], Paolo Manildo [aut, cre], Giacomo Zanella [aut]
Maintainer:	Paolo Manildo <paolo.manildo@studenti.unipd.it>
Repository:	CRAN
Date/Publication:	2026-10-22 10:40:02 UTC

Non-Reversible Markov Chains for Bayesian Mixture Models

Description

Provides Bayesian inference for finite and infinite mixture models using Markov chain Monte Carlo (MCMC) methods. Supports Gaussian, multivariate Gaussian, Poisson, Binomial, and multivariate Bernoulli mixture kernels with different covariance structures and flexible prior specifications. Provides tools for data generation, posterior density evaluation, clustering and similarity analysis, and advanced MCMC sampling strategies, including reversible and non-reversible samplers: Ascolani et al. (2026) <doi:10.48550/arXiv.2510.03226>, informed proposals: Zanella (2020) <doi:10.1080/01621459.2019.1585255>, collapsed Gibbs sampling: Diebolt and Robert (1994) <doi:10.1111/j.2517-6161.1994.tb01985.x>, and split-and-merge samplers: Dahl and Newcomb (2021) <doi:10.1080/00949655.2021.1998502>.

Details

The DESCRIPTION file:

Index of help topics:

compute_similarity_matrix
                        Compute the posterior similarity matrix
extract_params          Extract a Parameter from a 'nrMCmix' object.
generate_data           Simulate data from a Bayesian mixture-model
                        generative process.
nrMCmix                 MCMC inference for Bayesian mixture models
nrMCmix-package         Non-Reversible Markov Chains for Bayesian
                        Mixture Models
plot.nrMCmix            Plot diagnostics for an 'nrMCmix' model
posterior_density1      Evaluate the posterior predictive density
posterior_density2      Evaluate the bivariate posterior predictive
                        density
print.nrMCmix           Print an 'nrMCmix' model fit
print.summary.nrMCmix   Print an 'nrMCmix' model summary
priorK                  Prior distribution of the number of clusters
specify_model           Specify the model and prior distributions
specify_sampler         Specify the sampling algorithm
summary.nrMCmix         Summarize an nrMCmix model fit

The package provides functions for inference on Bayesian mixture models with conjugate priors using a variety of marginal MCMC schemes. The number of mixture components can either be fixed a priori or inferred using a Pitman–Yor mixture model. The package primarily implements the non-reversible methods developed by Ascolani et al. (2026), while also providing the split-and-merge methods of Dahl and Newcomb (2022) and novel random-neighborhood informed proposals for mixture models, inspired by Zanella (2020). These three families of algorithms can be combined to construct flexible sampling schemes for Bayesian mixture models.

Author(s)

Filippo Ascolani [aut] (<https://orcid.org/0009-0009-3899-8877>), Paolo Manildo [aut, cre] (<https://orcid.org/0009-0006-7056-7012>), Giacomo Zanella [aut] (<https://orcid.org/0000-0001-6727-9884>)

Maintainer: Paolo Manildo <paolo.manildo@studenti.unipd.it>

References

Ascolani, F. Manildo, P. and Zanella, G. (2026). Fast non-reversible samplers for Bayesian mixture models. doi:10.48550/arXiv.2510.03226.

Dahl, D. B. and Newcomb, S. (2021). Sequentially allocated merge-split samplers for conjugate Bayesian nonparametric models. doi:10.1080/00949655.2021.1998502.

Zanella, G. (2020). Informed proposals for discrete state-space models. doi:10.1080/01621459.2019.1585255.

Examples

#load the iris data
data(iris)

#fit the model
set.seed(132)
fit <- nrMCmix(as.matrix(iris[1:4]),
               K = 0,
               model = specify_model("MGaussian", dense = TRUE,
                                     alpha0 = 4.01,
                                     hyperpar_baseline = TRUE),
               sampler = specify_sampler())

#check the model fit
{
  op <- par(no.readonly = TRUE)
  par(mfrow = c(4,4), mar = c(1.1,1.1,1.1,1.1))
  for(i in 1:4){
    for(j in 1:4){
      if(i == j){
        #univariate plot
        plot(fit,idx1 = i, verbose = FALSE, gg = FALSE)
      }else{
        #bivariate plot
        plot(fit,idx1 = i, idx2 = j, verbose = FALSE, gg = FALSE)
      }
    }
  }
  par(op)
}

#check the summary
summary(fit)

#plot the number of components distribution
plot(fit,5,par_name = "K")

Compute the posterior similarity matrix

Description

Computes the posterior similarity matrix from sampled allocation configurations. Each entry gives the proportion of MCMC iterations in which two units are assigned to the same mixture component.

Usage

compute_similarity_matrix(c)

Arguments

Value

A numeric n by n symmetric matrix, where the $(i, j)$ entry is the proportion of sampled configurations in which units $i$ and $j$ belong to the same component.

Examples

#specify the model
model <- specify_model("Gaussian", sigma2 = 1, lambda0 = 0.1)

#generate some data
set.seed(123)
y <- generate_data(n = 100, K = 7, model = model)$y

#run a non-reversible split and merge MCMC for a DP mixture
fit <- nrMCmix(y, K = 0, model = model,
               sampler = specify_sampler(SAM = TRUE))

#compute the similarity matrix
S <- compute_similarity_matrix(fit$c)               

Extract a Parameter from a nrMCmix object.

Description

Extracts a parameter or a transformation of the parameters, from the output of the model-fitting procedure of function nrMCmix.

Usage

extract_params(
  X,
  par_name = "log_p",
  lab_name = NULL,
  which_coordinate = 1,
  f = function(x) x
)

Arguments

Value

A data frame containing the extracted parameter. The returned column is named according to lab_name. If lab_name is NULL, the name is generated from par_name and which_coordinate.

Examples

#specify the model
model <- specify_model("Gaussian", sigma2 = 1, lambda0 = 0.1)

#generate some data
set.seed(123)
y <- generate_data(n = 100, K = 7, model = model)$y

#run a non-reversible split and merge MCMC for a DP mixture
fit <- nrMCmix(y, K = 0, model = model,
 sampler = specify_sampler(SAM = TRUE))

#extract the posterior draws of the largest cluster proportion
largest_cluster <- extract_params(fit,"c","n1",f = function(x) max(table(x)) / length(x))

Simulate data from a Bayesian mixture-model generative process.

Description

Simulate data from a Bayesian mixture-model generative process.

Usage

generate_data(n = 100L, K = 0L, model = NULL, D = 2L, xi = 100)

Arguments

Value

A list containing y, the simulated observations, and groups, the corresponding mixture-component assignments.

Examples

#specify the model
model <- specify_model("Gaussian", sigma2 = 1, lambda0 = 0.1)

#generate some data
set.seed(123)
y <- generate_data(n = 100, K = 7, model = model)$y

MCMC inference for Bayesian mixture models

Description

Performs Markov chain Monte Carlo inference for Bayesian finite or nonparametric mixture models.

Usage

nrMCmix(
  y = NULL,
  K = 0L,
  model = NULL,
  sampler = NULL,
  init = NULL,
  N = 1000L,
  warm_up = 100L,
  thin = 1L,
  thin_scan = 0L,
  refresh_console = 0.1,
  save_parameters = TRUE,
  save_configurations = TRUE
)

Arguments

Value

An object of class nrMCmix, containing the MCMC output, including:

log_p

Log-posterior values across MCMC iterations.

n_pred_calls

Number of predictive density evaluations at each iteration.

entropy

Entropy of the allocation configuration.

K

Number of occupied components at each iteration for nonparametric mixtures.

hyperpars

Sampled hyperparameters.

acceptance_rates

Acceptance rates of the MCMC moves, when available.

max_excursions

Maximum number of excursions per scan, when available.

c

Sampled allocation configurations, or the final configuration if save_configurations = FALSE.

atoms

Sampled component parameters, when save_parameters = TRUE.

y

The input data.

Examples

#specify the model
model <- specify_model("Gaussian", sigma2 = 1, lambda0 = 0.1)

#generate some data
set.seed(123)
y <- generate_data(n = 100, K = 7, model = model)$y

#run a non-reversible split and merge MCMC for a DP mixture
fit <- nrMCmix(y, K = 0, model = model,
 sampler = specify_sampler(SAM = TRUE))

Plot diagnostics for an nrMCmix model

Description

Produces graphical summaries and diagnostics from a fitted nrMCmix object.

Usage

## S3 method for class 'nrMCmix'
plot(
  x,
  which_plot = 1L,
  line_col1 = NULL,
  line_col2 = NULL,
  x1_grid = NULL,
  x2_grid = NULL,
  idx1 = 1L,
  idx2 = NULL,
  B = 100L,
  l_grid = 101,
  n_breaks = NULL,
  verbose = TRUE,
  loss = salso::VI(),
  method = "complete",
  par_name = "log_p",
  lab_name = NULL,
  which_coordinate = 1L,
  f = function(x) x,
  ggplot2 = TRUE,
  ...
)

Arguments

Details

Five types of plots are available:

1. Data fit: for Gaussian, multivariate Gaussian, and Poisson kernels, displays the observed data together with posterior density draws and their pointwise mean. For multivariate Gaussian models, bivariate density plots can also be produced.

2. Dispersion plot: displays the observations coloured according to the estimated partition obtained using salso::salso(). For multivariate Gaussian models, bivariate dispersion plots can also be produced.

3. Posterior similarity matrix: displays the posterior similarity matrix computed from the sampled configurations.

4. Trace plot: displays the trace plot of a parameter (or function of parameters) across MCMC iterations, extracted using function extract_params.

5. Density plot: works as the previous one but displays the MCMC reconstruction of the posterior density.

For which_plot = 1, the available plots depend on the kernel. For multivariate Gaussian models, specifying idx2 produces a bivariate density contour plot. Otherwise, a univariate histogram is overlaid with posterior density draws and their pointwise posterior mean. The posterior draws of the density are obtained using the functions posterior_density1 for the univariate case and posterior_density2 for the bivariate case.

For which_plot = 2, the function estimates a representative partition using salso::salso() and colours observations according to their estimated cluster assignment. This plot requires multiple configuration draws to be available in X$c.

For which_plot = 3, the posterior similarity matrix is computed from X$c using function compute_similarity_matrix and it is displayed as a heatmap with the rows and columns ordered according to the output of hclust.

For which_plot = 4 and which_plot = 5, the requested parameter is extracted from the nrMCmix object using function extract_params.

Value

The function is called for its graphical side effect and does not return a value of interest.

Examples

#simulate multivariate gaussian data
set.seed(123)
y <- generate_data(n = 300, D = 2, K = 3,
                   model = specify_model("MGaussian", lambda0 = 0.9,
                                         dense = TRUE, alpha0 = 2.1))$y
plot(y)
#fit the model using nrMCmix
fit <- nrMCmix(y,model = specify_model("MGaussian", dense = TRUE,
                                       alpha0 = 2.01))

#model fit
{
  par(mfrow = c(2,2))
  
  #univariate density plot
  plot(fit, idx1 = 1, gg = FALSE)
  
  #bivariate density plots
  plot(fit, idx1 = 1, idx2 = 2, gg = FALSE)
  plot(fit, idx1 = 2, idx2 = 1, gg = FALSE)
  
  #univariate density plot
  plot(fit, idx1 = 2, gg = FALSE)
  par(mfrow = c(1,1))
}

#clustering
plot(fit, 2, idx1 = 1, idx2 = 2)

#posterior uncertainty using the similarity matrix
plot(fit,3)

#trace plot of the log posterior
plot(fit,4,par_name = "log_p")

#density plot of the numeber of cluster
plot(fit,5,par_name = "K")

Evaluate the posterior predictive density

Description

Computes the posterior predictive density at a set of evaluation points for a single variable. The density is obtained by averaging over the mixture components and, for nonparametric models, the remaining prior probability mass.

Usage

posterior_density1(draw, hyperpars, xj, j = 0L)

Arguments

Value

A numeric vector containing the posterior predictive density evaluated at each element of xj.

Examples

#specify the model
model <- specify_model("MGaussian", sigma2 = 1, lambda0 = 0.1)

#generate some data
set.seed(123)
y <- generate_data(n = 100, K = 7, model = model)$y

#run a non-reversible split and merge MCMC for a DP mixture
fit <- nrMCmix(y, K = 0, model = model,
 sampler = specify_sampler(SAM = TRUE))
 
#compute the density for the first draw
dens <- posterior_density1(fit$atoms[[1]],fit$hyperpars,seq(-5,5,l=101))

Evaluate the bivariate posterior predictive density

Description

Computes the posterior predictive density on a two-dimensional grid of evaluation points for a pair of variables under the Gaussian mixture model. For nonparametric models, the contribution of the remaining prior probability mass is also included.

Usage

posterior_density2(draw, hyperpars, xj, xl, j = 0L, l = 1L)

Arguments

Value

A numeric matrix containing the posterior predictive density on the grid defined by xj and xl.

Examples

#specify the model
model <- specify_model("MGaussian", sigma2 = 1, lambda0 = 0.1)

#generate some data
set.seed(123)
y <- generate_data(n = 100, K = 7, model = model)$y

#run a non-reversible split and merge MCMC for a DP mixture
fit <- nrMCmix(y, K = 0, model = model,
 sampler = specify_sampler(SAM = TRUE))
 
#compute the density for the first draw
dens <- posterior_density2(fit$atoms[[1]],fit$hyperpars,
                           seq(-5,5,l=101), seq(-4,4,l = 101),
                           0,1)

Print an nrMCmix model fit

Description

Prints the model specifications stored in an nrMCmix object.

Usage

## S3 method for class 'nrMCmix'
print(x, ...)

Arguments

Value

Invisibly returns NULL. Called primarily for its printing side effect.

Examples

#simulate data
set.seed(123)
y <- generate_data(n = 100, K = 3, 
                   model = specify_model("Poisson", 
                                         alpha0 = 10, 
                                         beta0 = 3))$y[,1]

#fit the model
fit <- nrMCmix(y,K=3,model = specify_model("Poisson"))

#check the specifics
fit

Print an nrMCmix model summary

Description

Prints the posterior summaries and MCMC diagnostics contained in an nrMCmix summary object.

Usage

## S3 method for class 'summary.nrMCmix'
print(x, ...)

Arguments

Value

Invisibly returns NULL. Called primarily for its printing side effect.

Prior distribution of the number of clusters

Description

Computes the prior probability distribution of the number of clusters induced by a two-parameter nonparametric prior.

Usage

priorK(alpha = 1, delta = 0, n = 100L)

Arguments

Value

A numeric vector of length n, where the K-th element is the prior probability that the number of clusters is K.

Examples

  priorK(alpha = 1, delta = 0, n = 100)

Specify the model and prior distributions

Description

Constructs a list specifying the kernel, clustering prior, and kernel-specific prior distributions used by the model.

Usage

specify_model(
  kernel = "Gaussian",
  alpha = NULL,
  delta = 0,
  hyperpar_alpha = FALSE,
  hyperpar_delta = FALSE,
  a0 = 0.01,
  b0 = 15,
  lwr_bound = 0,
  upr_bound = 1,
  mu0 = NULL,
  lambda0 = 0.01,
  sigma2 = NULL,
  alpha0 = 0.01,
  beta0 = 0.01,
  dense = FALSE,
  B0 = NULL,
  hyperpar_baseline = FALSE,
  m0 = NULL,
  V0 = 10,
  L1 = 0.1,
  L2 = 0.1,
  A1 = 0.1,
  A2 = 0.1,
  B1 = 0.1,
  B2 = 0.1,
  nu0 = 0.1,
  Psi0 = NULL
)

Arguments

Value

A list containing the model specification and the corresponding prior and hyperprior parameters.

Examples

model <- specify_model(kernel = "Gaussian")

model <- specify_model(
  kernel = "MGaussian",
  dense = FALSE,
  hyperpar_alpha = TRUE
)

Specify the sampling algorithm

Description

Constructs a list specifying the options and tuning parameters of the sampling algorithm.

Usage

specify_sampler(
  reversible = FALSE,
  informed = FALSE,
  gibbs = FALSE,
  SAM = FALSE,
  NUSAMS = FALSE,
  xi = 0.5,
  s = 0,
  m = 20L,
  g = 0L,
  n_restricted_steps = 0L,
  thin_SAM = 1L,
  NU_weights = NULL,
  return_weights = FALSE,
  slice_max_iter = 10L,
  scale = 0.5
)

Arguments

Value

A list containing the selected sampling options and the corresponding tuning parameters.

Examples

sampler <- specify_sampler()

sampler <- specify_sampler(
  reversible = FALSE,
  informed = TRUE,
  m = 20
)

sampler <- specify_sampler(
  SAM = TRUE,
  NUSAMS = TRUE
)

Summarize an nrMCmix model fit

Description

Computes posterior summary statistics and MCMC diagnostics for an nrMCmix model fit. Summaries are provided for the log-posterior, entropy, and, when available, the number of mixture components. Acceptance rates and excursion statistics are also reported when available.

Usage

## S3 method for class 'nrMCmix'
summary(object, ...)

Arguments

Value

A list containing:

posterior_summaries

A data frame containing the minimum, first quartile, median, mean, third quartile, maximum, and effective sample size (ESS) for each monitored quantity.

split_rate

Mean split-move acceptance rate, when available.

merge_rate

Mean merge-move acceptance rate, when available.

local_acceptance_rate

Mean local-move acceptance rate, when available.

excursions

Median maximum excursion, when available.

Examples

#simulate data
set.seed(123)
y <- generate_data(n = 100, K = 3, 
                   model = specify_model("Poisson", 
                                         alpha0 = 10, 
                                         beta0 = 3))$y[,1]

#fit the model
fit <- nrMCmix(y,K=3,model = specify_model("Poisson"))

#check the summary
summary(fit)