--- title: "Introduction to the nrMCmix package" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to the nrMCmix package} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: REFERENCES.bib header-includes: - \usepackage{amsmath} - \usepackage{mathtools} - \usepackage{amsfonts} - \usepackage{amssymb} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Introduction The `nrMCmix` package provides an interface for Bayesian inference in mixture models via _marginal_ schemes. More precisely, it implements different versions of the _non-reversible_ sampler described by @ascolani2026fastnonreversiblesamplerbayesian. # How the package works The model to be fitted is specified using the `specify_model` function, while the sampling algorithm is defined using `specify_sampler`. The `nrMCmix` function then runs a Markov chain invariant with respect to the _collapsed_ posterior distribution of the previously specified mixture model. The `print`, `summary`, and `plot` methods can then be used to summarize the information contained in the posterior distribution. ## 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 [@pitman1997two], 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 @dahl2022sequentially, $\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*} ### Data generation with the `generate_data` function Using the output of the `specify_model` function, it is possible to simulate data from the specified model by passing it to the model argument. The argument `K` specifies the number of mixture components; setting `K = 0` corresponds to the infinite mixture case. The argument `D` specifies the number of dimensions and is used when kernel is either `"MGaussian"` or `"MBernoulli"`. Finally, `n` specifies the desired sample size, while `xi` specifies the Poisson rate used to generate the number of trials for each observation when simulating data under the binomial kernel. ## 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 [@diebolt1994estimation], 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 @ascolani2026fastnonreversiblesamplerbayesian, 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 @ascolani2026fastnonreversiblesamplerbayesian. Setting `informed = TRUE` employs a random-neighborhood informed proposal [@zanella2020informed], 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 [@dahl2022sequentially;@jain2002split]. 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 @neal2003slice. 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 [@plummer2006coda]. 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_ [@meilua2007comparing], implemented in the `salso` `R` package [@dahl2026salso]. - `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`. # Practical examples This section provides practical examples illustrating the use of the `nrMCmix` package for Bayesian mixture modelling with non-reversible samplers. ```{r, echo = FALSE, eval = TRUE} set.seed(123) ``` ## Galaxy dataset ```{r setup, eval = TRUE, echo = FALSE} library(nrMCmix) ``` 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. ```{r load, echo = TRUE, eval = TRUE} 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. ```{r galaxy, echo = TRUE, eval = TRUE} 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) ``` Printing the resulting object provides a summary of the model and sampling algorithm used. ```{r print, echo = TRUE, eval = TRUE} fit print(summary(fit)) ``` We can further investigate the output of the model using the `plot` method. ```{r plot, echo = TRUE, eval = TRUE, fig.align='center',fig.width=5, fig.height=4} plot(fit, n_breaks = 33) ``` We first assess the convergence of the algorithm by examining the trace plot of the logarithm of the (unnormalized) posterior density. ```{r log_p, echo = TRUE, eval = TRUE, fig.align='center',fig.width=5, fig.height=4} 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. ```{r K, fig.align='center',fig.width=5, fig.height=4} plot(fit,4, par_name = "K") ``` Finally, we examine the posterior distributions of the concentration and discount parameters. ```{r densities, fig.align='center',fig.width=5, fig.height=3} 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. ```{r iris} 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. ```{r iris_plot, fig.align='center',fig.width=5, fig.height=5} { 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. ```{r iris_plot2, fig.align='center',fig.width=5, fig.height=5} plot(fit,3) ``` ## 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. ```{r binomial} 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. ```{r, PY_binomial} 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. ```{r binomial_plot, fig.align='center',fig.width=5, fig.height=5} plot(fit) ``` 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 @ascolani2026fastnonreversiblesamplerbayesian. - The composition of the previous non-reversible sampler with global updates using the split and merge method of @dahl2022sequentially. - The composition of the split and merge method with an informed version [@zanella2020informed] 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. ```{r binomial_comparison} #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) ``` ## Bibliography