piv_rel: Performing the pivotal relabelling step and computing the...
In pivmet: Pivotal Methods for Bayesian Relabelling and k-Means Clustering

piv_rel

R Documentation

Performing the pivotal relabelling step and computing the relabelled posterior estimates

Description

This function allows to perform the pivotal relabelling procedure described in Egidi et al. (2018) and to obtain the relabelled posterior estimates.

Usage

piv_rel(mcmc)

Arguments

mcmc

The output of the MCMC sampling from piv_MCMC.

Details

Prototypical models in which the label switching problem arises are mixture models, as explained in the Details section of the piv_MCMC function.

These models are unidentified with respect to an arbitrary permutation of the labels 1,...,k. Relabelling means permuting the labels at each iteration of the Markov chain in such a way that the relabelled chain can be used to draw inferences on component-specific parameters.

We assume here that a MCMC sample is obtained for the posterior distribution of a Gaussian mixture model–for instance via piv_MCMC function–with a prior distribution which is labelling invariant. Furthermore, suppose that we can find k units, one for each group, which are (pairwise) separated with (posterior) probability one (that is, the posterior probability of any two of them being in the same group is zero). It is then straightforward to use the k units, called pivots in what follows and denoted by the indexes i_1,…,i_k, to identify the groups and to relabel the chains: for each MCMC iteration h=1,…, H (H corresponds to the argument nMC) and group j=1,…,k, set

[μ_j]_h=[μ_{[Z_{i_{j}}]_h}]_h;

[Z_{i}]_h=j \mbox{ for } i:[Z_i]_h=[Z_{i_{j}}]_h.

The applicability of this strategy is limited by the existence of the pivots, which is not guaranteed. The existence of the pivots is a requirement of the method, meaning that its use is restricted to those chains—or those parts of a chain—for which the pivots are present. First, although the model is based on a mixture of k components, each iteration of the chain may imply a different number of non-empty groups. Let then [k]_h ≤q k be the number of non-empty groups at iteration h,

[k]_h = \#\{j: [Z_i]_h=j\mbox{ for some }i\},

where \#A is the cardinality of the set A. Hence, the relabelling procedure outlined above can be used only for the subset of the chain for which [k]_h=k; let it be

\mathcal{H}_k=\{h:[k]_h= k\},

which correspond to the argument true.iter given by piv_MCMC. This means that the resulting relabelled chain is not a sample (of size H) from the posterior distribution, but a sample (of size \#\mathcal{H}_k) from the posterior distribution conditional on there being (exactly) k non-empty groups. Even if k non-empty groups are available, however, there may not be k perfectly separated units. Let us define

\mathcal{H}^{*}_k=\{ h\in\mathcal{H}_k : \exists r,s \mbox{ s.t. } [Z_{i_r}]_h=[Z_{i_s}]_h \}

that is, the set of iterations where (at least) two pivots are in the same group. In order for the pivot method to be applicable, we need to exclude iterations \mathcal{H}^{*}_k; that is, we can perform the pivot relabelling on \mathcal{H}_k- \mathcal{H}^{*}_{k}, corresponding to the argument final_it.

Value

This function gives the relabelled posterior estimates–both mean and medians–obtained from the Markov chains of the MCMC sampling.

`final_it`	The final number of valid MCMC iterations, as explained in Details.
`final_it_p`	The proportion of final valid MCMC iterations.
`rel_mean`	The relabelled chains of the means: a `final_it` \times k matrix for univariate data, or a `final_it` \times D \times k array for multivariate data.
`rel_sd`	The relabelled chains of the sd's: a `final_it` \times k matrix for univariate data, or a `final_it` \times D matrix for multivariate data.
`rel_weight`	The relabelled chains of the weights: a `final_it` \times k matrix.

Author(s)

Leonardo Egidi legidi@units.it

References

Egidi, L., Pappadà, R., Pauli, F. and Torelli, N. (2018). Relabelling in Bayesian Mixture Models by Pivotal Units. Statistics and Computing, 28(4), 957-969.

Examples


#Univariate simulation
## Not run: 
N   <- 250
nMC <- 2500
k   <- 3
p   <- rep(1/k,k)
x   <- 3
stdev <- cbind(rep(1,k), rep(20,k))
Mu    <- seq(-trunc(k/2)*x,trunc(k/2)*x,length=k)
W     <- c(0.2,0.8)
sim   <- piv_sim(N = N, k = k, Mu = Mu,
                 stdev = stdev, W=W)
res   <- piv_MCMC(y = sim$y, k =k, nMC = nMC)
rel   <- piv_rel(mcmc=res)

## End(Not run)

#Bivariate simulation
## Not run: 
N <- 200
k <- 3
D <- 2
nMC <- 5000
M1  <- c(-.5,8)
M2  <- c(25.5,.1)
M3  <- c(49.5,8)
Mu  <- matrix(rbind(M1,M2,M3),c(k,2))
Sigma.p1 <- diag(D)
Sigma.p2 <- 20*diag(D)
W <- c(0.2,0.8)
sim <- piv_sim(N = N, k = k, Mu = Mu,
               Sigma.p1 = Sigma.p1,
               Sigma.p2 = Sigma.p2, W = W)
res <- piv_MCMC(y = sim$y, k = k, nMC = nMC)
rel <- piv_rel(mcmc = res)
piv_plot(y=sim$y, mcmc=res, rel_est = rel, type="chains")
piv_plot(y=sim$y, mcmc=res, rel_est = rel,
         type="hist")

## End(Not run)

pivmet documentation built on March 7, 2023, 6:34 p.m.