variational_fSAN {SANvi} | R Documentation |
Mean Field Variational Bayes estimation of fSAN
Description
variational_fSAN
is used to perform posterior inference under the finite shared atoms nested (fSAN) model with Gaussian likelihood (originally proposed in D'Angelo et al., 2023).
The model uses finite Dirichlet mixtures for both the distributional and observational levels of the model.
Usage
variational_fSAN(y, group, maxL = 30, maxK = 20,
m0 = 0, tau0 = .01, lambda0 = 3, gamma0 = 2,
alpha_bar = .005, beta_bar = .005,
epsilon = 1e-6, seed = NULL, maxSIM = 1e5,
warmstart = TRUE,verbose = FALSE)
Arguments
y |
Numerical vector of observations (required). |
group |
Numerical vector of the same length of |
maxL , maxK |
integers, the upper bounds for the observational and distributional clusters to fit, respectively |
m0 , tau0 , lambda0 , gamma0 |
Hyperparameters on |
alpha_bar |
the hyperparameter of the symmetric distributional Dirichlet distribution. |
beta_bar |
the hyperparameter of the symmetric observational Dirichlet distribution. |
epsilon |
the tolerance that drives the convergence criterion adopted as stopping rule |
seed |
random seed to control the initialization. |
maxSIM |
the maximum number of CAVI iteration to perform. |
warmstart |
logical, if |
verbose |
logical, if |
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(\alpha/K,\dots,\alpha/K)
.
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(\beta/L,\dots,\beta/L)
for all k = 1,\dots,K
.
Here, the dimension L
is fixed.
Value
variational_fSAN
returns a list of class SANvb
containing four objects:
-
model
: name of the fitted model. -
params
: list containing the data and the parameters used in the simulation. Details below. -
sim
: list containing the simulated values (optimized variational parameters). 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 fitted distributional and observational clusters.
m0, tau0, lambda0, gamma0
Model hyperparameters.
epsilon, seed
The threshold controlling the convergence criterion and the random seed adopted to replicate the run.
alpha_bar, beta_bar
the hyperparameters governing all the finite Dirichlet distributions at the distributional and observational level.
Simulated values:
sim
is a list with the following components:
theta_l
Matrix of size (L,4). Each row is a posterior variational estimate of the four normal-inverse gamma hyperparameters.
Elbo_val
Vector containing the values of the ELBO.
XI
A list of length J. Each element is a matrix of size (N, L) 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, K). 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.
alpha_bar_k
Vector of updated variational parameters of the Dirichlet distributions governing the distributional clustering.
beta_bar_lk
Matrix of updated variational parameters of the Dirichlet distributions governing the observational clustering (arranged by column).
References
D’Angelo, L., Canale, A., Yu, Z., and Guindani, M. (2023). Bayesian nonparametric analysis for the detection of spikes in noisy calcium imaging data. Biometrics, 79(2), 1370–1382. DOI: 10.1111/biom.13626
Examples
set.seed(123)
y <- c(rnorm(50),rnorm(50,5))
g <- rep(1:2,rep(50,2))
est <- variational_fSAN(y, g, verbose = FALSE,
epsilon = 1e-2, maxL=15, maxK=10)