fit_fSAN: Fit the Finite Shared Atoms Mixture Model

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

Fit the Finite Shared Atoms Mixture Model

Description

fit_fSAN fits the finite shared atoms nested (fSAN) 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_fSAN(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). a_dirichlet The hyperparameter of the symmetric distributional Dirichlet distribution. The default is 1/maxK. 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 iteration 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 verbose Logical, if TRUE the iterations are printed.

Details

Data structure

The finite 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,\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 Gamma(\lambda_0, \gamma_0) (shape, rate).

Clustering

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

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

The distribution of the probabilities is (\pi_1,\dots,\pi_{K})\sim Dirichlet_K(a,\dots,a). Here, the dimension K is fixed.

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,\dots,K \, ; \: l = 1,\dots,L.

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

Value

fit_fSAN 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.

• a_dirichlet: Provided value for a.

• b_dirichlet: Provided value for b.

• 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) posterior variational probability of assignment 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_dirichlet_k: Vector of updated variational parameters of the Dirichlet distribution governing the distributional clustering.

• b_dirichlet_lk: Matrix of updated variational parameters of the Dirichlet distributions governing the observational clustering (arranged by column).

• 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 of each observational cluster (\mu_1,\dots\mu_L).

• sigma2: Matrix of size (nrep, maxL). Each row is a posterior sample of the variance parameter of 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_{\code{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_{\code{maxL},k}) for distributional cluster k.

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_fSAN(y, group = g, est_method = "VI", vi_param = list(n_runs = 1))
out_vi

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

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