sample_fiSAN_sparsemix {SANple}R Documentation

Sample fiSAN with sparse mixtures

Description

sample_fiSAN_sparsemix is used to perform posterior inference under the finite-infinite shared atoms nested (fiSAN) model with Gaussian likelihood. The model uses a Dirichlet process mixture prior at the distributional level, and a sparse (overfitted) Dirichlet mixture (Malsiner-Walli et al., 2016) at the observational level. The algorithm for the nonparametric component is based on the slice sampler for DPM of Kalli, Griffin and Walker (2011).

Usage

sample_fiSAN_sparsemix(nrep, burn, y, group, 
                 maxK = 50, maxL = 50, 
                 m0 = 0, tau0 = 0.1, lambda0 = 3, gamma0 = 2,
                 hyp_alpha1 = 1, hyp_alpha2 = 1,
                 hyp_beta = 10, 
                 eps_beta = NULL, 
                 alpha = NULL, beta = NULL,
                 warmstart = TRUE, nclus_start = NULL,
                 mu_start = NULL, sigma2_start = NULL,
                 M_start = NULL, S_start = NULL,
                 alpha_start = NULL, beta_start = NULL,
                 progress = TRUE, seed = NULL)

Arguments

nrep

Number of MCMC iterations.

burn

Number of discarded iterations.

y

Vector of observations.

group

Vector of the same length of y indicating the group membership (numeric).

maxK

Maximum number of distributional clusters K (default = 50).

maxL

Maximum number of observational clusters L (default = 50).

m0, tau0, lambda0, gamma0

Hyperparameters on (\mu, \sigma^2) \sim NIG(m_0, \tau_0, \lambda_0,\gamma_0). Default is (0, 0.1, 3, 2).

hyp_alpha1, hyp_alpha2

If a random \alpha is used, (hyp_alpha1,hyp_alpha2) specify the hyperparameters (default = (1,1)). The prior is \alpha ~ Gamma(hyp_alpha1, hyp_alpha2).

hyp_beta, eps_beta

If a random \beta is used, hyp_beta specifies the hyperparameter (default = 10). The prior is \beta ~ Gamma(hyp_beta, hyp_beta*maxL), following Malsiner-Walli et al. (2016). In this case, eps_beta is the tuning parameter of the MH step.

alpha

Distributional DP parameter if fixed (optional). The distribution is \pi\sim GEM (\alpha).

beta

Observational Dirichlet parameter if fixed (optional). The distribution is Dirichlet( rep(beta, maxL) ). Notice that beta should be small to ensure sparsity (e.g. beta = 0.01)

warmstart, nclus_start

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)

mu_start, sigma2_start, M_start, S_start, alpha_start, beta_start

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. alpha_start, alpha_start are numeric.

progress

show a progress bar? (logical, default TRUE).

seed

set a fixed seed.

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 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,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 Dirichlet_L(\beta,\dots,\beta) for all k = 1,2,\dots.

Value

sample_fiSAN_sparsemix returns four objects:

Data and parameters: params is a list with the following components:

nrep

Number of MCMC iterations.

y, group

Data and group vectors.

maxK, maxL

Maximum number of distributional and observational clusters.

m0, tau0, lambda0, gamma0

Model hyperparameters.

(hyp_alpha1,hyp_alpha2) or alpha

Either the hyperparameters on \alpha (if \alpha random), or the value for \alpha (if fixed).

(hyp_beta,eps_beta) or beta

Either the hyperparameter on \beta and MH step size (if \beta random), or the value for \beta (if fixed).

Simulated values: sim is a list with the following components:

mu

Matrix of size (nrep, maxL). Each row is a posterior sample of the mean parameter for each observational cluster (\mu_1,\dots\mu_L).

sigma2

Matrix 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_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_{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_{L,k}) for distributional cluster k.

alpha

Vector of length nrep of posterior samples of the parameter \alpha.

beta

Vector of length nrep of posterior samples of the parameter \beta.

maxK

Vector of length nrep of the number of distributional DP components used by the slice sampler.

References

Kalli, M., Griffin, J.E., and Walker, S.G. (2011). Slice Sampling Mixture Models, Statistics and Computing, 21, 93–105. <doi:10.1007/s11222-009-9150-y>

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>

Sethuraman, A.J. (1994). A Constructive Definition of Dirichlet Priors, Statistica Sinica, 4, 639–650.

Examples

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_fiSAN_sparsemix(nrep = 500, burn = 200, y = y, group = g, 
                              nclus_start = 2,
                              maxK = 20, maxL = 20,
                              beta = 0.01)
out
[Package SANple version 0.1.1 Index]