fit_CAM: 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())

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.

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$

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

Run this code

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

Run the code above in your browser using DataLab