label.switching (version 1.8)

sjw: Probabilistic relabelling algorithm

Description

Function to apply the probabilistic relabelling strategy of Sperrin et al (2010). The concept here is to treat the MCMC output as observed data, while the unknown permutations need to be applied to each mcmc data point is treated as unobserved data with associated uncertainty. Then, an EM-type algorithm estimates the weights for each permutation per MCMC data point.

Usage

sjw(mcmc.pars, z, complete, x, init, threshold, maxiter)

Arguments

mcmc.pars

\(m\times K\times J\) array of simulated MCMC parameters.

z

\(m\times n\) integer array of the latent allocation vectors generated from an MCMC algorithm.

complete

function that returns the complete log-likelihood of the mixture model.

x

\(n\)-dimensional data vector/array

init

An (optional) index pointing at the MCMC iteration whose parameters will initialize the algorithm. If it is less or equal to zero, the overall MCMC mean will be used for initialization.

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.

Value

permutations

\(m\times K\) dimensional array of permutations

iterations

integer denoting the number of iterations until convergence

status

returns the exit status

Details

Let \(x=(x_1,\ldots,x_n)\) denote the observed data and \(\boldsymbol{w},\boldsymbol{\theta}\) denote the mixture weights and component specific parameters, respectively. Assume that \(K\) is the the number of components. Then, $$L(\boldsymbol{w},\boldsymbol{\theta}|\boldsymbol{x})=\prod_{i=1}^{n}\sum_{k=1}^{K}w_k f_k(x_i|\theta_k),$$ \(i=1,\ldots,n\) is the observed likelihood of the mixture model. Given the latent allocation variables \(\boldsymbol{z}=(z_1,\ldots,z_n)\), the complete likelihood of the model is defined as: $$L_c(\boldsymbol{w},\boldsymbol{\theta}|\boldsymbol{x},\boldsymbol{z})=\prod_{i=1}^{n}w_{z_{i}}f_{z_{i}}(x_i|\theta_{z_{i}}).$$ Then, complete corresponds to the log of \(L_c\) and should take as input the following: a vector of \(n\) allocations, the observed data and the parameters of the model as a \(K\times J\) array where \(J\) corresponds to the different parameter types of the model. See the example for an implementation at a univariate normal mixture.

References

Sperrin M, Jaki T and Wit E (2010). Probabilistic relabelling strategies for the label switching problem in Bayesian mixture models. Statistics and Computing, 20(3), 357-366.

See Also

permute.mcmc, label.switching

Examples

# NOT RUN {
#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}.  
data("mcmc_output")
mcmc.pars<-data_list$"mcmc.pars"
z<-data_list$"z"
K<-data_list$"K"
x<-data_list$"x"

# mcmc parameters are stored to array \code{mcmc.pars}
# mcmc.pars[,,1]: simulated means of the two components
# mcmc.pars[,,2]: simulated variances
# mcmc.pars[,,3]: simulated weights 
# The number of different parameters for the univariate
# normal mixture is equal to J = 3: means, variances 
# and weights. The generated allocations variables are 
# stored to \code{z}. The observed data is stored to \code{x}.  
# The complete data log-likelihood is defined as follows:
complete.normal.loglikelihood<-function(x,z,pars){
#	x: data (size = n)
#	z: allocation vector (size = n)
#	pars: K\times J vector of normal mixture parameters:
#		pars[k,1] = mean of the k-normal component
#		pars[k,2] = variance of the k-normal component
#		pars[k,3] = weight of the k-normal component
#			k = 1,...,K
	g <- dim(pars)[1] #K (number of mixture components)
	n <- length(x)	#this denotes the sample size
	logl<- rep(0, n)	
 	logpi <- log(pars[,3])
	mean <- pars[,1]
	sigma <- sqrt(pars[,2])
	logl<-logpi[z] + dnorm(x,mean = mean[z],sd = sigma[z],log = TRUE)
	return(sum(logl))
}

#run the algorithm:
run<-sjw(mcmc = mcmc.pars,z = z, 
complete = complete.normal.loglikelihood,x = x, init=0,threshold = 1e-4)
# apply the permutations returned by typing:
reordered.mcmc<-permute.mcmc(mcmc.pars,run$permutations)
# reordered.mcmc[,,1]: reordered means of the two components
# reordered.mcmc[,,2]: reordered variances 
# reordered.mcmc[,,3]: reordered weights
# }