| piv_rel {pivmet} | 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 |
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,\ldots,i_k, to identify the groups and to
relabel the chains:
for each MCMC iteration h=1,\ldots, H (H corresponds to
the argument nMC) and group
j=1,\ldots,k, set
[\mu_j]_h=[\mu_{[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 \leq 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 |
rel_sd |
The relabelled chains of the sd's: a |
rel_weight |
The relabelled chains of the weights: a |
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)