Stephens' algorithm
Stephens (2000) developed a relabelling algorithm that makes the permuted sample points to agree as much as possible on the n\times K matrix of classification probabilities, using the Kullback-Leibler divergence. The algorithm's input is the matrix of allocation probabilities for each MCMC iteration.
stephens(p, threshold, maxiter)
p |
m\times n \times K dimensional array of allocation probabilities of the n observations among the K mixture components, for each iteration t = 1,…,m of the MCMC algorithm. |
threshold |
An (optional) positive number controlling the convergence criterion. Default value: 1e-6. |
maxiter |
An (optional) integer controlling the max number of iterations. Default value: 100. |
For a given MCMC iteration t=1,…,m, let w_k^{(t)} and θ_k^{(t)}, k=1,…,K denote the simulated mixture weights and component specific parameters respectively. Then, the (t,i,k) element of p corresponds to the conditional probability that observation i=1,…,n belongs to component k and is proportional to p_{tik} \propto w_k^{(t)} f(x_i|θ_k^{(t)}), k=1,…,K, where f(x_i|θ_k) denotes the density of component k. This means that:
p_{tik} = \frac{w_k^{(t)} f(x_i|θ_k^{(t)})}{w_1^{(t)} f(x_i|θ_1^{(t)})+… + w_K^{(t)} f(x_i|θ_K^{(t)})}.
In case of hidden Markov models, the probabilities w_k should be replaced with the proper left (normalized) eigenvector of the state-transition matrix.
permutations |
m\times K dimensional array of permutations |
iterations |
integer denoting the number of iterations until convergence |
status |
returns the exit status |
Panagiotis Papastamoulis
Stephens, M. (2000). Dealing with label Switching in mixture models. Journal of the Royal Statistical Society Series B, 62, 795-809.
#load a toy example: MCMC output consists of the random beta model
# applied to a normal mixture of \code{K=2} components. The number
# of observations is equal to \code{n=5}. The number of MCMC samples
# is equal to \code{m=300}. The matrix of allocation probabilities
# is stored to matrix \code{p}.
data("mcmc_output")
# mcmc parameters are stored to array \code{mcmc.pars}
mcmc.pars<-data_list$"mcmc.pars"
# mcmc.pars[,,1]: simulated means of the two components
# mcmc.pars[,,2]: simulated variances
# mcmc.pars[,,3]: simulated weights
# the computed allocation matrix is p
p<-data_list$"p"
run<-stephens(p)
# apply the permutations returned by typing:
reordered.mcmc<-permute.mcmc(mcmc.pars,run$permutations)
# reordered.mcmc[,,1]: reordered means of the components
# reordered.mcmc[,,2]: reordered variances
# reordered.mcmc[,,3]: reordered weightsPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.