fit_CAM: Fit the Common Atoms Mixture Model

In sanba: Fitting Shared Atoms Nested Models via MCMC or Variational Bayes

Fit the Common Atoms Mixture Model

Description

fit_CAM fits the common atoms mixture model (CAM) of Denti et al. (2023) with Gaussian kernels and normal-inverse gamma priors on the unknown means and variances. The function returns an object of class SANmcmc or SANvi depending on the chosen computational approach (MCMC or VI).

Usage

fit_CAM(y, group, est_method = c("VI", "MCMC"),
         prior_param = list(),
         vi_param = list(),
         mcmc_param = list())

Arguments

y	Numerical vector of observations (required).
group	Numerical vector of the same length of y, indicating the group membership (required).
est_method	Character, specifying the preferred estimation method. It can be either "VI" or "MCMC".
prior_param	A list containing m0, tau0, lambda0, gamma0 Hyperparameters on (\mu, \sigma^2) \sim NIG(m_0, \tau_0, \lambda_0,\gamma_0). The default is (0, 0.01, 3, 2). hyp_alpha1, hyp_alpha2 If a random \alpha is used, (hyp_alpha1, hyp_alpha2) specify the hyperparameters. The default is (1,1). The prior is \alpha ~ Gamma(hyp_alpha1, hyp_alpha2). alpha Distributional DP parameter if fixed (optional). The distribution is \pi\sim GEM (\alpha). hyp_beta1, hyp_beta2 If a random \beta is used, (hyp_beta1, hyp_beta2) specify the hyperparameters. The default is (1,1). The prior is \beta ~ Gamma(hyp_beta1, hyp_beta2). beta Observational DP parameter if fixed (optional). The distribution is \omega_k \sim GEM (\beta).
vi_param	A list of variational inference-specific settings containing maxL, maxK Integers, the upper bounds for the observational and distributional clusters to fit, respectively. The default is (50, 20). epsilon The threshold controlling the convergence criterion. n_runs Number of starting points considered for the estimation. seed Random seed to control the initialization. maxSIM The maximum number of CAVI iterations to perform. warmstart Logical, if TRUE, the observational means of the cluster atoms are initialized with a k-means algorithm. verbose Logical, if TRUE the iterations are printed.
mcmc_param	A list of MCMC inference-specific settings containing nrep, burn Integers, the number of total MCMC iterations, and the number of discarded iterations, respectively. maxL, maxK Integers, the upper bounds for the observational and distributional clusters to fit, respectively. The default is (50, 20). seed Random seed to control the initialization. warmstart Logical, if TRUE, the observational means of the cluster atoms are initialized with a k-means algorithm. If FALSE, the starting points can be passed through the parameters nclus_start, mu_start, sigma2_start, M_start, S_start, alpha_start, beta_start verbose Logical, if TRUE the iterations are printed.

Details

The common atoms mixture model is used to perform inference in nested settings, where the data are organized into J groups. The data should be continuous observations (Y_1,\dots,Y_J), where each Y_j = (y_{1,j},\dots,y_{n_j,j}) contains the n_j observations from group j, for j=1,\dots,J. The function takes as input the data as a numeric vector y in this concatenated form. Hence, y should be a vector of length n_1+\dots+n_J. The group parameter is a numeric vector of the same size as y, indicating the group membership for each individual observation. Notice that with this specification, the observations in the same group need not be contiguous as long as the correspondence between the variables y and group is maintained.

Model

The data are modeled using a Gaussian likelihood, where both the mean and the variance are observational cluster-specific, i.e.,

y_{i,j}\mid M_{i,j} = l \sim N(\mu_l,\sigma^2_l)

where M_{i,j} \in \{1,2,\dots\} is the observational cluster indicator of observation i in group j. The prior on the model parameters is a normal-inverse gamma distribution (\mu_l,\sigma^2_l)\sim NIG (m_0,\tau_0,\lambda_0,\gamma_0), i.e., \mu_l\mid\sigma^2_l \sim N(m_0, \sigma^2_l / \tau_0), 1/\sigma^2_l \sim Gamma(\lambda_0, \gamma_0) (shape, rate).

Clustering

The model clusters both observations and groups. The clustering of groups (distributional clustering) is provided by the allocation variables S_j \in \{1,2,\dots\}, with

Pr(S_j = k \mid \dots ) = \pi_k \qquad \text{for } \: k = 1,2,\dots

The distribution of the probabilities is \{\pi_k\}_{k=1}^{\infty} \sim GEM(\alpha), where GEM is the Griffiths-Engen-McCloskey distribution of parameter \alpha, which characterizes the stick-breaking construction of the DP (Sethuraman, 1994).

The clustering of observations (observational clustering) is provided by the allocation variables M_{i,j} \in \{1,2,\dots\}, with

Pr(M_{i,j} = l \mid S_j = k, \dots ) = \omega_{l,k} \qquad \text{for } \: k = 1,2,\dots \, ; \: l = 1,2,\dots

The distribution of the probabilities is \{\omega_{l,k}\}_{l=1}^{\infty} \sim GEM(\beta) for all k = 1,2,\dots

Value

fit_CAM returns a list of class SANvi, if est_method = "VI", or SANmcmc, if est_method = "MCMC". The list contains the following elements:

model

Name of the fitted model.

params

List containing the data and the parameters used in the simulation. Details below.

sim

List containing the optimized variational parameters or the simulated values. Details below.

time

Total computation time.

Data and parameters: params is a list with the following components:

• y, group, Nj, J: Data, group labels, group frequencies, and number of groups.

• K, L: Number of distributional and observational mixture components.

• m0, tau0, lambda0, gamma0: Model hyperparameters.

• (hyp_alpha1, hyp_alpha2) or alpha: hyperparameters on \alpha (if \alpha random); or provided value for \alpha (if fixed).

• (hyp_beta1, hyp_beta2) or beta: hyperparameters on \beta (if \beta random); or provided value for \beta (if fixed).

• seed: The random seed adopted to replicate the run.

• epsilon, n_runs: The threshold controlling the convergence criterion and the number of starting points considered.

• nrep, burnin: If est_method = "MCMC", the number of total MCMC iterations, and the number of discarded ones.

Simulated values: depending on the algorithm, it returns a list with the optimized variational parameters or a list with the chains of the simulated values.

Variational inference: sim is a list with the following components:

• theta_l: Matrix of size (maxL, 4). Each row is a posterior variational estimate of the four normal-inverse gamma hyperparameters.

• XI: A list of length J. Each element is a matrix of size (Nj, maxL) containing the posterior variational probability of assignment of the i-th observation in the j-th group to the l-th OC, i.e., \hat{\xi}_{i,j,l} = \hat{\mathbb{Q}}(M_{i,j}=l).

• RHO: Matrix of size (J, maxK). Each row is a posterior variational probability of assignment of the j-th group to the k-th DC, i.e., \hat{\rho}_{j,k} = \hat{\mathbb{Q}}(S_j=k).

• a_tilde_k, b_tilde_k: Vector of updated variational parameters of the beta distributions governing the distributional stick-breaking process.

• a_bar_lk, b_bar_lk: Matrix of updated variational parameters of the beta distributions governing the observational stick-breaking process (arranged by column).

• conc_hyper: If the concentration parameters are chosen to be random, a vector with the four updated hyperparameters.

• Elbo_val: Vector containing the values of the ELBO.

MCMC inference: sim is a list with the following components:

• mu: Matrix of size (nrep, maxL). Each row is a posterior sample of the mean parameter for each observational cluster (\mu_1,\dots\mu_L).

• sigma2: Matrix of size (nrep, maxL). Each row is a posterior sample of the variance parameter for each observational cluster (\sigma^2_1,\dots\sigma^2_L).

• obs_cluster: Matrix of size (nrep, n), with n = length(y). Each row is a posterior sample of the observational cluster allocation variables (M_{1,1},\dots,M_{n_J,J}).

• distr_cluster: Matrix of size (nrep, J), with J = length(unique(group)) Each row is a posterior sample of the distributional cluster allocation variables (S_1,\dots,S_J).

• pi: Matrix of size (nrep, maxK). Each row is a posterior sample of the distributional cluster probabilities (\pi_1,\dots,\pi_{maxK}).

• omega: 3-d array of size (maxL, maxK, nrep). Each slice is a posterior sample of the observational cluster probabilities. In each slice, each column k is a vector (of length maxL) observational cluster probabilities (\omega_{1,k},\dots,\omega_{maxL,k}) for distributional cluster k.

• alpha: Vector of length nrep of posterior samples of the parameter \alpha.

• beta: Vector of length nrep of posterior samples of the parameter \beta.

• maxK: Vector of length nrep of the number of distributional DP components used by the slice sampler.

• maxL: Vector of length nrep of the number of observational DP components used by the slice sampler.

References

Denti, F., Camerlenghi, F., Guindani, M., and Mira, A. (2023). A Common Atoms Model for the Bayesian Nonparametric Analysis of Nested Data. Journal of the American Statistical Association, 118(541), 405-416. DOI: 10.1080/01621459.2021.1933499

Sethuraman, A.J. (1994). A Constructive Definition of Dirichlet Priors, Statistica Sinica, 4, 639–650.

Examples

set.seed(123)
y <- c(rnorm(60), rnorm(40, 5))
g <- rep(1:2, rep(50, 2))

out_vi <- fit_CAM(y, group = g, est_method = "VI", vi_param = list(n_runs = 1,
                  epsilon = .01 ))
out_vi

out_mcmc <- fit_CAM(y = y, group = g, est_method = "MCMC",
                    mcmc_param = list(nrep = 50, burn = 20))
out_mcmc

sanba documentation built on Aug. 8, 2025, 6:15 p.m.