gibbs: Gibbs sampler and Chib's evidence approximation for a generic univariate mixture of normal distributions

Description

This function implements a regular Gibbs sampling algorithm on the posterior distribution associated with a mixture of normal distributions, taking advantage of the missing data structure. It then runs an averaging of the simulations over all permutations of the component indices in order to avoid incomplete label switching and to validate Chib's representation of the evidence. This function reproduces gibbsnorm as its first stage, however it may be much slower because of its second stage.

Usage

gibbs(niter, datha, mix)

Arguments

niter

number of Gibbs iterations

datha

sample vector

mix

list made of k, number of components, p, mu, and sig, starting values of the paramerers, all of size k (see example below)

Value

knumber of components in the mixture (superfluous as it is invariant over the execution of the R code)
mumatrix of the Gibbs samples on the $\mu_i$ parameters
sigmatrix of the Gibbs samples on the $\sigma_i$ parameters
progmatrix of the Gibbs samples on the mixture weights
lolikvector of the observed log-likelihoods along iterations
chibdenodenominator of Chib's approximation to the evidence (see example below)

References

Chib, S. (1995) Marginal likelihood from the Gibbs output.J. American Statist. Associ. 90, 1313-1321.

Examples

Run this code

faithdata=faithful[,1]
mu=rnorm(3,mean=mean(faithdata),sd=sd(faithdata)/10)
sig=1/rgamma(3,shape=10,scale=var(faithdata))
mix=list(k=3,p=rdirichlet(par=rep(1,3)),mu=mu,sig=sig)
resim3=gibbs(100,faithdata,mix)
lulu=order(resim3$lolik)[100]
lnum1=resim3$lolik[lulu]
lnum2=sum(dnorm(resim3$mu[lulu,],mean=mean(faithdata),sd=resim3$sig[lulu,],log=TRUE)+
dgamma(resim3$sig[lulu,],10,var(faithdata),log=TRUE)-2*log(resim3$sig[lulu,]))+
sum((rep(0.5,mix$k)-1)*log(resim3$p[lulu,]))+
lgamma(sum(rep(0.5,mix$k)))-sum(lgamma(rep(0.5,mix$k)))
lchibapprox3=lnum1+lnum2-log(resim3$deno)