fit_fiSAN: Fit the Finite-Infinite Shared Atoms Mixture Model
In sanba: Fitting Shared Atoms Nested Models via MCMC or Variational Bayes

fit_fiSAN

R Documentation

Fit the Finite-Infinite Shared Atoms Mixture Model

Description

fit_fiSAN fits the finite-infinite shared atoms nested (fiSAN) mixture model 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_fiSAN(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 \text{GEM} (\alpha)`. `b_dirichlet` The hyperparameter of the symmetric observational Dirichlet distribution. The default is 1/`maxL`.
`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`. `verbose` Logical, if `TRUE` the iterations are printed.

Details

Data structure

The finite-infinite 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:

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

where M_{i,j} \in \{1,\dots,L \} 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 \text{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,\dots,L\}, with:

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

The distribution of the probabilities is (\omega_{1,k},\dots,\omega_{L,k})\sim \text{Dirichlet}_L(b,\dots,b) for all k = 1,2,\dots. Here, the dimension L is fixed.

Value

fit_fiSAN 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).
b_dirichlet: Provided value for b.
seed: The random seed adopted to replicate the run.
epsilon, n_runs: If est_method = "VI", 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), the posterior variational assignment probabilities \hat{\mathbb{Q}}(M_{i,j}=l).
RHO: Matrix of size (J, maxK), with the posterior variational assignment probabilities \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.
conc_hyper: If the concentration parameter is random, this contains its updated hyperparameters.
b_dirichlet_lk: Matrix of updated variational parameters of the Dirichlet distributions governing observational clustering.
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) with samples of the observational cluster means.
sigma2: Matrix of size (nrep, maxL) with samples of the observational cluster variances.
obs_cluster: Matrix of size (nrep, n) with posterior samples of the observational cluster allocations.
distr_cluster: Matrix of size (nrep, J) with posterior samples of the distributional cluster allocations.
pi: Matrix of size (nrep, maxK) with posterior samples of the distributional cluster weights.
omega: Array of size (maxL, maxK, nrep) with observational cluster weights.
alpha: Vector of length nrep with posterior samples of \alpha.
maxK: Vector of length nrep with number of active distributional components.

Examples

set.seed(123)
y <- c(rnorm(60), rnorm(40, 5))
g <- rep(1:2, rep(50, 2))
plot(density(y[g==1]), xlim = c(-5,10), main = "Group-specific density")
lines(density(y[g==2]), col = 2)

out_vi <- fit_fiSAN(y, group = g, est_method = "VI",
                    vi_param = list(n_runs = 1))
out_vi

out_mcmc <- fit_fiSAN(y = y, group = g, est_method = "MCMC")
out_mcmc

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