sample_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 overfitted (sparse) Dirichlet mixtures (Malsiner-Walli et al., 2016) at both the observational and distributional level.
sample_fSAN(nrep, burn, y, group,
maxK = 50, maxL = 50,
m0 = 0, tau0 = 0.1, lambda0 = 3, gamma0 = 2,
alpha = 0.01, beta = 0.01,
warmstart = TRUE, nclus_start = NULL,
mu_start = NULL, sigma2_start = NULL,
M_start = NULL, S_start = NULL,
progress = TRUE, seed = NULL)sample_fSAN returns 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 (MCMC chains). Details below.
time: total computation time.
Data and parameters:
params is a list with the following components:
nrepNumber of MCMC iterations.
y, groupData and group vectors.
maxK, maxLMaximum number of distributional and observational clusters.
m0, tau0, lambda0, gamma0Model hyperparameters.
Simulated values:
sim is a list with the following components:
muMatrix of size (nrep, maxL).
Each row is a posterior sample of the mean parameter for each observational cluster \((\mu_1,\dots\mu_L)\).
sigma2Matrix 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_clusterMatrix 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_clusterMatrix 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)\).
piMatrix of size (nrep, maxK).
Each row is a posterior sample of the distributional cluster probabilities \((\pi_1,\dots,\pi_{K})\).
omega3-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_{L,k})\) for distributional cluster \(k\).
alphaVector of length nrep of posterior samples of the parameter \(\alpha\).
betaVector of length nrep of posterior samples of the parameter \(\beta\).
Number of MCMC iterations.
Number of discarded iterations.
Vector of observations.
Vector of the same length of y indicating the group membership (numeric).
Maximum number of distributional clusters \(K\) (default = 50).
Maximum number of observational clusters \(L\) (default = 50).
Hyperparameters on \((\mu, \sigma^2) \sim NIG(m_0, \tau_0, \lambda_0,\gamma_0)\). Default is (0, 0.1, 3, 2).
Distributional Dirichlet parameter (default alpha = 0.01). The distribution is Dirichlet( rep(alpha, maxK) ).
Observational Dirichlet parameter (default beta = 0.01). The distribution is Dirichlet( rep(beta, maxL) ).
Initialization of the observational clustering.
warmstart is logical parameter (default = TRUE) of whether a kmeans clustering should be used to initialize the chains.
An initial guess of the number of observational clusters can be passed via the nclus_start parameter (optional)
Starting points of the MCMC chains (optional). Default is nclus_start = min(c(maxL, 30)).
mu_start, sigma2_start are vectors of length maxL.
M_start is a vector of observational cluster allocation of length N.
S_start is a vector of observational cluster allocation of length J.
show a progress bar? (logical, default TRUE).
set a fixed seed.
Data structure
The overfitted mixture common atoms 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 univariate 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,\dots,\alpha)\).
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,\dots,\beta)\) for all \(k = 1,\dots,K\).
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>
Malsiner-Walli, G., Frühwirth-Schnatter, S. and Grün, B. (2016). Model-based clustering based on sparse finite Gaussian mixtures. Statistics and Computing 26, 303–324. <doi:10.1007/s11222-014-9500-2>
set.seed(123)
y <- c(rnorm(40,0,0.3), rnorm(20,5,0.3))
g <- c(rep(1,30), rep(2, 30))
plot(density(y[g==1]), xlim = c(-5,10))
lines(density(y[g==2]), col = 2)
out <- sample_fSAN(nrep = 500, burn = 200, y = y, group = g,
nclus_start = 2,
maxK = 20, maxL = 20,
alpha = 0.01, beta = 0.01)
out
Run the code above in your browser using DataLab