Model specification through the specify_model function

Let \(y = (y_1,\dotsc,y_n)\) denote an observed sample of \(n\) data points in a \(D\)-dimensional space, the Bayesian mixture models implemented in package nrMCmix assume that \[ y \,\mid\, P \overset{iid}\sim \int_\Theta f(y\,\mid\,\theta)\text{d}P(\theta), \label{eq:model} \] where \(f\) is the mixture kernel indexed by the parameter \(\theta \in\Theta\), and \(P\) is the random mixing distribution. The latter is assumed to be either a probability mass function over \(K\) mixture components, with a \(\text{Dir}_{K}(\alpha)\) prior, or a realization of a Pitman–Yor process (Pitman and Yor 1997), with concentration parameter \(\alpha\) and discount parameter \(\delta\).

The kernel argument specify the kernel family \(f\). In particular,

• kernel = "MGaussian" specifies a multivariate Gaussian kernel. When dense = TRUE, a Normal-Inverse-Wishart distribution is assigned to the parameters \(\theta\): \[ \theta = (\mu,\Sigma) \sim \mathrm{NIW}(\mu_0,\lambda_0,a_0,B_0). \] When dense=FALSE, , there are two possible specifications. Setting sigma2 to a non-NULL value corresponds to assuming a spherical covariance structure and assigning a multivariate normal prior to \(\theta = \mu\): \[ \mu \sim \mathrm{N}_D(\mu_0,\lambda_0^{-1}\sigma^2I_D). \] On the other hand, if sigma2 = NULL, an independent Normal-Inverse-Gamma prior is assigned to each dimension: \[ (\mu_d,\sigma^2_d)_{d = 1}^D \overset{iid}\sim \mathrm{NIG}(\mu_0,\lambda_0,\alpha_0,\beta_0). \]

• kernel = "Gaussian" specifies a univariate Gaussian kernel. As in the multivariate case, setting sigma2 = NULL corresponds to assuming a Normal-Inverse-Gamma prior, whereas setting sigma2 to a non-NULL value corresponds to the simpler Gaussian specification.

• kernel = "Poisson" specifies a Poisson kernel, with \(\theta = \lambda\) representing the rate parameter. A conjugate Gamma prior is assumed: \[ \lambda \sim \mathrm{Ga}(\alpha_0,\beta_0). \]

• kernel = "Binomial" specifies a binomial kernel with a conjugate Beta prior: \[ \theta = p \sim \mathrm{Beta}(\alpha_0,\beta_0). \]

• kernel = "MBernoulli" can be used to model multivariate dichotomous data. In particular, \[ f(y_1,\dotsc,y_D\,\mid\,\theta) = \prod_{d = 1}^D \theta_d^{y_d}(1-\theta_d)^{1-y_d}, \] with \[ \theta_1,\dotsc,\theta_D \overset{iid}\sim \mathrm{Beta}(\alpha_0,\beta_0). \]

The arguments hyperpar_alpha, hyperpar_delta, and hyperpar_baseline, when set to TRUE, assign hyperpriors to \(\alpha\), \(\delta\), and the parameters of the prior on \(\theta\), respectively.

In the finite mixture case, a hyperprior on \(\alpha\) assumes that it is a symmetric vector with a single unique entry: \[ \alpha \sim \mathrm{Ga}(a_0,b_0). \] For the infinite mixture case, following David B. Dahl and Newcomb (2022), \(\alpha\) is a scalar with prior \[ \alpha \sim \mathrm{Unif}(a_0,b_0), \] with default values \(a_0 = 0\) and \(b_0 = 15\). The discount parameter is assigned the prior \[ \delta \sim \mathrm{Unif}(\text{lwr_bound},\text{upr_bound}) \] with default values \(\text{lwr_bound} = 0\) and \(\text{upr_bound} = 1\).

Depending on the choice of kernel, different hyperpriors are assigned to the parameters of the baseline distribution. These are summarized below: \[\begin{align*} &\mu_0 \sim \mathrm{N}_D(m_0,V_0I_D),\, \lambda_0 \sim \mathrm{Ga}(L_1,L_2),\, \alpha_0 \sim \mathrm{Ga}(A_1,A_2),\\ &\beta_0 \sim \mathrm{Ga}(B_1,B_2),\, B_0^{-1} \sim \mathrm{Wishart}(\nu_0,\Psi_0). \end{align*}\]

Sampler specification throuth the specify_sampler function

The nrMCmix package provides different Markov chain Monte Carlo (MCMC) algorithms targeting the collapsed posterior distribution of the Bayesian mixture models described above. The specify_sampler function allows the user to select the desired sampling algorithm. In particular, setting gibbs = TRUE selects the classical marginal Gibbs sampler (Diebolt and Robert 1994), which updates the cluster allocation of each observation one at a time. Setting reversible to either TRUE or FALSE selects one of the two main samplers described in Ascolani, Manildo, and Zanella (2026), namely, the reversible and the non-reversible one.

In this case, the parameter xi regulates the refreshment rate. At each iteration, with probability \(\texttt{xi}/n\), the directional variables are randomly shuffled. This mechanism is only used to avoid periodicity issues and has no practical impact provided that xi is set to a sufficiently small value. The parameter s regulates the type of non-reversible sampler. Setting s to a value greater than zero corresponds to the modified non-reversible sampler described in Ascolani, Manildo, and Zanella (2026).

Setting informed = TRUE employs a random-neighborhood informed proposal (Zanella 2020), with parameter m providing an upper bound on the neighborhood size and g regulating the type of balancing function. In particular, g = 0 corresponds to \(g(t) = \sqrt{t}\), g = 1 corresponds to the Metropolis function \(g(t) = \min{1,t}\), while g = 2 corresponds to the Barker function \(g(t) = \frac{t}{1+t}\).

When K = 0 is passed to nrMCmix (see below), the infinite mixture case allows for split-and-merge global updates (David B. Dahl and Newcomb 2022; Jain 2002). These updates are enabled by setting SAM = TRUE. The parameter n_restricted_steps regulates the number of restricted Gibbs steps used to improve the acceptance probability of the global proposal, while thin_SAM specifies the number of subsequent global transitions before refining the cluster allocation with local updates. Moreover, setting NUSAMS = TRUE uses non-uniform weights to select the observations involved in split or merge proposals. The sampling weights can either be provided through the NU_weights argument or estimated during the warm-up phase.

Finally, slice_max_iter and scale regulate the maximum number of iterations and the scale, respectively, of the stepping-out univariate slice sampler of Neal (2003). These updates are employed when updating univariate hyperparameters. For multivariate hyperparameters, exact Gibbs updates are obtained by first simulating the collapsed mixture parameters.

MCMC inference using the nrMCmix function

The main function of the nrMCmix package is also named nrMCmix. It requires the object y, which can be either a matrix or a vector, and K, which specifies the number of mixture components. Setting K = 0 corresponds to the infinite mixture case. The model and sampler arguments take as input the outputs of the specify_model and specify_sampler functions, respectively. These determine the type of Bayesian mixture model to be fitted to the data and the sampling algorithm to be employed.

The number of MCMC iterations, discarded samples, and thinned samples are controlled by the arguments N, warm_up, and thin, respectively. The argument refresh_console controls how frequently the user receives updates on the progress of the Markov chain. The argument thin_scan controls the length of each sequence of local updates, with a default value equal to \(n\). Finally, setting save_parameters = TRUE and save_configurations = TRUE returns the posterior draws of the collapsed parameters (i.e., \(\theta\)) and hyperparameters, and of the cluster allocations, respectively.

The function returns an nrMCmix object, which is a named list containing, among others, the following elements:

• log_p: the trace of the logarithm of the (unnormalized) posterior density.

• entropy: the trace of the entropy of the data partition.

• n_pred_calls: the trace of the number of predictive density evaluations performed after each sequence of local updates.

• K: in the infinite mixture case, the trace of the number of clusters.

• hyperpars: the hyperparameters of the model.

• acceptance_rates: the acceptance rates of the Metropolis-Hastings proposals. When SAM = TRUE, the first columns refer to split proposals, the second to merge proposals, and the third to local updates.

• max_excursions: the trace of the maximum excursion, i.e., the maximum number of consecutive iterations without a direction flip, recorded during a sequence of local updates.

• c: the posterior draws of the clustering configurations. If save_configurations = FALSE, only the final configuration is returned.

• atoms: a list of length N. Each element contains the sampled values of \(\theta\) and, when applicable, the hyperparameters.

• y, the data used to fit the model.

Posterior summaries with the functions print, summary and plot

An nrMCmix object can be analyzed using the print, summary, and plot methods. The print method provides a concise summary of the model and sampling algorithm used. The summary method provides summaries of the log_p, entropy, and K traces, including quantiles, means, and effective sample sizes, computed using the coda package (Plummer et al. 2006).

The plot method provides several graphical summaries of the model fit. The type of plot is selected using the which_plot argument.

• which_plot = 1 overlays either univariate or bivariate posterior densities on a histogram or scatter plot, respectively. The argument x1_grid specifies the grid on which the density is computed for the first dimension, while x2_grid specifies the corresponding grid for the second dimension in the bivariate case. The arguments idx1 and idx2 specify the first and second dimensions, respectively. The argument B controls the number of posterior draws used to estimate the density. The argument l_grid specifies the length of the equally spaced grid automatically generated when the corresponding grid is NULL. Finally, n_breaks controls the number of histogram bins.

• which_plot = 2 provides a scatter plot of one or two dimensions of the data, with points colored according to the estimated optimal partition obtained by minimizing the loss specified by the loss argument. The default is the variation of information loss (Meilă 2007), implemented in the salso R package (David B. Dahl et al. 2026).

• which_plot = 3 provides a graphical representation of the posterior similarity matrix computed using the compute_similarity_matrix function.

• which_plot = 4 provides a one-dimensional trace plot for a parameter or hyperparameter of interest. Other posterior summaries, such as log_p, entropy, and K, can also be plotted. The parameter of interest is specified through the par_name argument. When the parameter is multidimensional, which_coordinate specifies the coordinate of interest. For matrix-valued parameters, such as \(\Sigma\) in the multivariate Gaussian case with dense = TRUE, which_coordinate can be a vector of length two specifying the row and column of the desired matrix element. The f argument can be used to apply a transformation to the selected quantity before plotting, allowing the trace of any function of the parameter of interest to be visualized. For example, setting f = exp allows the trace of a parameter stored on the log scale to be plotted on its original scale.

Note: The extraction and transformation of posterior parameters described above is handled internally by the extract_params function. This function can also be used independently to extract a parameter, select a specific coordinate, and optionally apply a transformation to the extracted quantity.

• which_plot = 5 is analogous to which_plot = 4, but displays a smooth posterior density estimate obtained using the density function from the stats package.

Furthermore, all plots can be customized using the ... argument, which allows graphical parameters to be passed to the underlying plotting functions. Finally, setting ggplot2 = TRUE constructs the plot using the grammar of graphics framework, as implemented by ggplot2.

Galaxy dataset

We start by analyzing the classic galaxy dataset, available from the bmixture package. The dataset contains \(82\) observations of the velocities (in \(1000\) km/second) of distant galaxies diverging from our own, observed in six well-separated conic sections of the Corona Borealis.

data("galaxy", package = "bmixture")

We specify the model and sampling algorithm using the specify_model and specify_sampler functions, respectively. In particular, we consider a Pitman-Yor mixture of Gaussian kernels with unknown variance. All hyperparameters are assigned hyperpriors, avoiding the need to specify them manually. We run the MCMC algorithm for \(20,000\) iterations, discarding the first \(10,000\) iterations as warm-up and retaining one sample every \(10\) iterations. For the sampling algorithm, we use non-reversible local updates together with split-and-merge global updates.

model <- specify_model("Gaussian",
                       hyperpar_alpha = TRUE,
                       hyperpar_delta = TRUE,
                       hyperpar_baseline = TRUE)

sampler <- specify_sampler(SAM = TRUE)

fit <- nrMCmix(y = galaxy, model = model, sampler = sampler,
               warm_up = 1000, thin = 10)
#> Warm-up currently at 10%
#> Warm-up currently at 20%
#> Warm-up currently at 30%
#> Warm-up currently at 40%
#> Warm-up currently at 50%
#> Warm-up currently at 60%
#> Warm-up currently at 70%
#> Warm-up currently at 80%
#> Warm-up currently at 90%
#> Warm-up currently at 100%
#> Sampling currently at 10%
#> Sampling currently at 20%
#> Sampling currently at 30%
#> Sampling currently at 40%
#> Sampling currently at 50%
#> Sampling currently at 60%
#> Sampling currently at 70%
#> Sampling currently at 80%
#> Sampling currently at 90%
#> Sampling currently at 100%

Printing the resulting object provides a summary of the model and sampling algorithm used.

fit
#> Pitman-Yor mixture model, using a Gaussian kernel with unknown variance. 
#> The concentration parameter is assigned a uniform prior over (0.010000,15.000000).
#> The discount parameter is assigned a uniform prior over (0.000000,1.000000).
#> A hyperprior is specified for the baseline distribution. 
#> The MCMC algorithm uses a non-reversible allocation sampler (with refreshment rate 0.500 every full data scan). 
#> After 82 local updates, 1 Split and Merge updates using 0 restricted steps are performed.
print(summary(fit))
#>                      Min.      1st Qu.       Median         Mean     3rd Qu.
#> log_p        -231.0428427 -185.0344248 -173.3346085 -174.9915898 -160.829318
#> entropy         0.4452904    0.4452904    0.6729786    0.8215492    1.153795
#> n_components    3.0000000    3.0000000    4.0000000    4.6420000    6.000000
#>                   Max.      ESS
#> log_p        -155.3917 484.8063
#> entropy         1.8091 165.5282
#> n_components   12.0000 138.0556
#> 
#> Split rate:  0.0057
#> Merge rate:  0.0046
#> Local scan acceptance rate:  0.07168
#> Median maximum number of excurions per scan:  3

We can further investigate the output of the model using the plot method.

plot(fit, n_breaks = 33)
#> 
#> Computing the posterior densities... done!

We first assess the convergence of the algorithm by examining the trace plot of the logarithm of the (unnormalized) posterior density.

plot(fit,4,par_name = "log_p")

We can also inspect the trace of the number of clusters. In this example, the chain explores configurations containing at least three clusters.

plot(fit,4, par_name = "K")

Finally, we examine the posterior distributions of the concentration and discount parameters.

par(mfrow = c(1,2))
plot(fit,5,par_name = "alpha", lab_name = "alpha")

plot(fit,5,par_name = "delta", lab_name = "delta")

par(mfrow = c(1,1))

Iris dataset

We next consider the iris dataset, which contains measurements, in centimeters, of sepal length, sepal width, petal length, and petal width for \(50\) flowers from each of three iris species: Iris setosa, Iris versicolor, and Iris virginica. We use a multivariate Gaussian kernel with independent Normal-Inverse-Gamma priors across dimensions.

data(iris)

y <- as.matrix(iris[1:4])
fit <- nrMCmix(y, sampler = specify_sampler(SAM = TRUE),
               model = specify_model("MGaussian", dense = FALSE),
               warm_up = 1000, thin = 10, refresh_console = 0)

We visualize the fit using the plot function. Here, we display all pairwise combinations of the four dimensions. On the diagonal, we show the corresponding univariate posterior density overlaid on the empirical distribution, while off the diagonal we show the corresponding bivariate posterior density.

{
  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){
      plot(fit,idx1 = i, verbose = FALSE, ggplot2 = FALSE)
    }else{
      plot(fit,idx1 = i, idx2 = j, verbose = FALSE, ggplot2 = FALSE)
    }
  }
}
par(op)
}

We can also use which_plot = 3 to visualize the posterior similarity matrix. Since the observations in the iris dataset are ordered by species, the block structure of the matrix provides a visual representation of the inferred clustering.

plot(fit,3)
#> 
#> Computing the similarity matrix... done!

Simulated dataset

We conclude with a simulated example. We choose the binomial kernel and generate \(1,000\) observations from a mixture with \(K = 7\) components.

model <- specify_model("Binomial", alpha0 = 2.1, beta0 = 1.4)

set.seed(123)
y <- generate_data(n = 1000, K = 7, model = model)$y

We then fit a Pitman-Yor mixture with hyperpriors on the concentration parameter, discount parameter, and baseline distribution.

fit <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(SAM = TRUE),
               warm_up = 1000,
               refresh_console = 0)

Once again, we use the plot method to assess the fitted model. In this case, we visualize the data together with the inferred clustering, using the variation of information loss function to obtain a point estimate of the partition.

plot(fit)
#> Warning in plot.nrMCmix(fit): 'which_plot = 1' applies only to Gaussian or
#> Poisson kernels, switched to 'which_plot = 2'.
#> 
#> Estimating the partition... done!

Now, we use the function extract_params to more easily handle the comparison of different samplers according to the effective sample size per number of evaluation of the the partition entropy. We compare four algorithms:

• The marginal Gibbs sampler.

• The non-reversible sampler of Ascolani, Manildo, and Zanella (2026).

• The composition of the previous non-reversible sampler with global updates using the split and merge method of David B. Dahl and Newcomb (2022).

• The composition of the split and merge method with an informed version (Zanella 2020) of the non-reversible local sampler.

We fix the upper bound for the neighborhood size of the latter scheme as m = 10. As can be seen, this provides a better mixing chain but with an higher computational cost.

#set the seed
set.seed(123)

#marginal gibbs
MG <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(gibbs = TRUE),
               warm_up = 100, N = 1000,refresh_console = 0)

#non-reversible
NR <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(),
               warm_up = 100, N = 1000,refresh_console = 0)

#non-reversible + split and merge
SAM_NR <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(SAM = TRUE),
               warm_up = 100, N = 1000,refresh_console = 0)

#non-reversible informed + split and merge
SAM_NR_inf10 <- nrMCmix(y,model = specify_model("Binomial",
                                       hyperpar_alpha = TRUE,
                                       hyperpar_delta = TRUE,
                                       hyperpar_baseline = TRUE),
               sampler = specify_sampler(SAM = TRUE,
                                         informed = TRUE,
                                         m = 10),
               warm_up = 100, N = 1000,refresh_console = 0)

#get the entropy traces
entropy <- cbind(MG = extract_params(MG,"entropy"),
                 NR = extract_params(NR,"entropy"),
                 SAM = extract_params(SAM_NR,"entropy"),
                 SAM_NR_inf10 = extract_params(SAM_NR_inf10,
                                             "entropy"))
      
#get the overall computational cost 
#(number of predictive density evaluations)
cost <- cbind(MG = extract_params(MG,"n_pred_calls"),
              NR = extract_params(NR,"n_pred_calls"),
              SAM = extract_params(SAM_NR,"n_pred_calls"),
              SAM_NR_inf10 = extract_params(SAM_NR_inf10,
                                            "n_pred_calls"))
             

df <- data.frame(log(apply(entropy,2,coda::effectiveSize)),
                 log(apply(cost,2,sum)),
                 log(apply(entropy,2,coda::effectiveSize))-
                   log(apply(cost,2,sum)))
colnames(df) <- c("log(ESS)","log(n° target evaluations)","log(ESS) per target evaluation")
rownames(df) <- c("Marginal Gibbs", "Non-reversible",
                  "Non-reversible + Split and Merge",
                  "Non-reversible-informed (m = 20) + Split and Merge")
knitr::kable(df)