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).
fit_fSAN(y, group, est_method = c("VI", "MCMC"),
prior_param = list(),
vi_param = list(),
mcmc_param = list())fit_fSAN returns a list of class SANvi, if est_method = "VI", or SANmcmc, if est_method = "MCMC". The list contains the following elements:
modelName of the fitted model.
paramsList containing the data and the parameters used in the simulation. Details below.
simList containing the optimized variational parameters or the simulated values. Details below.
timeTotal 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\).
Numerical vector of observations (required).
Numerical vector of the same length of y, indicating the group membership (required).
Character, specifying the preferred estimation method. It can be either "VI" or "MCMC".
A list containing
m0, tau0, lambda0, gamma0Hyperparameters on \((\mu, \sigma^2) \sim NIG(m_0, \tau_0, \lambda_0,\gamma_0)\). The default is (0, 0.01, 3, 2).
a_dirichletThe hyperparameter of the symmetric distributional Dirichlet distribution. The default is 1/maxK.
b_dirichletThe hyperparameter of the symmetric observational Dirichlet distribution. The default is 1/maxL.
A list of variational inference-specific settings, containing
maxL, maxKIntegers, the upper bounds for the observational and distributional clusters to fit, respectively. The default is (50, 20).
epsilonThe threshold controlling the convergence criterion.
n_runsNumber of starting points considered for the estimation.
seedRandom seed to control the initialization.
maxSIMThe maximum number of CAVI iteration to perform.
warmstartLogical, if TRUE, the observational means of the cluster atoms are initialized with a k-means algorithm.
verboseLogical, if TRUE the iterations are printed.
A list of MCMC inference-specific settings, containing
nrep, burnIntegers, the number of total MCMC iterations, and the number of discarded iterations, respectively.
maxL, maxKIntegers, the upper bounds for the observational and distributional clusters to fit, respectively. The default is (50, 20).
seedRandom seed to control the initialization.
warmstartLogical, 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
verboseLogical, if TRUE the iterations are printed.
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.
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
Run the code above in your browser using DataLab