R: Sample CAM

sample_CAM {SANple}

R Documentation

Sample CAM

Description

sample_CAM is used to perform posterior inference under the common atoms model (CAM) of Denti et al. (2023) with Gaussian likelihood. The model uses Dirichlet process mixtures (DPM) at both the observational and distributional levels. The implemented algorithm is based on the nested slice sampler of Denti et al. (2023), based on the algorithm of Kalli, Griffin and Walker (2011).

Usage

sample_CAM(nrep, y, group, 
           maxK = 50, maxL = 50, 
           m0 = 0, tau0 = 0.1, lambda0 = 3, gamma0 = 2,
           hyp_alpha1 = 1, hyp_alpha2 = 1,
           hyp_beta1 = 1, hyp_beta2 = 1,
           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.
`y`	Vector of observations.
`group`	Vector of the same length of y indicating the group membership (numeric).
`maxK`	Maximum number of distributional clusters (default = 50).
`maxL`	Maximum number of observational clusters (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 is (1,1). The prior is `\alpha` ~ Gamma(`hyp_alpha1`, `hyp_alpha2`).
`hyp_beta1`, `hyp_beta2`	If a random `\beta` is used, (`hyp_beta1`, `hyp_beta2`) specify the hyperparameters. Default is (1,1). The prior is `\beta` ~ Gamma(`hyp_beta1`, `hyp_beta2`).
`alpha`	Distributional DP parameter if fixed (optional). The distribution is `\pi\sim GEM (\alpha)`.
`beta`	Observational DP parameter if fixed (optional). The distribution is `\omega_k \sim GEM (\beta)`.
`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). Default is `nclus_start = min(c(maxL, 30))`.
`mu_start`, `sigma2_start`, `M_start`, `S_start`, `alpha_start`, `beta_start`	Starting points of the MCMC chains (optional). `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 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,2,\dots\} 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,2,\dots\}, with

Pr(M_{i,j} = l \mid S_j = k, \dots ) = \omega_{l,k} \qquad \text{for } \: k = 1,2,\dots \, ; \: l = 1,2,\dots

The distribution of the probabilities is \{\omega_{l,k}\}_{l=1}^{\infty} \sim GEM(\beta) for all k = 1,2,\dots

Value

sample_CAM 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:

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_beta1,hyp_beta2) or beta: Either the hyperparameters on \beta (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_{maxL,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.
maxL: Vector of length nrep of the number of observational DP components used by the slice sampler.

References

Denti, F., Camerlenghi, F., Guindani, M., and Mira, A. (2023). A Common Atoms Model for the Bayesian Nonparametric Analysis of Nested Data. Journal of the American Statistical Association, 118(541), 405-416. <doi:10.1080/01621459.2021.1933499>

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>

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_CAM(nrep = 500, y = y, group = g,
                  nclus_start = 2,
                  maxL = 20, maxK = 20)
out

[Package SANple version 0.1.0 Index]