gibbs {bayess}

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

k

number of components in the mixture (superfluous as it is invariant over the execution of the R code)

mu

matrix of the Gibbs samples on the \mu_i parameters

sig

matrix of the Gibbs samples on the \sigma_i parameters

prog

matrix of the Gibbs samples on the mixture weights

lolik

vector of the observed log-likelihoods along iterations

chibdeno

denominator 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.

See Also

gibbsnorm

Examples

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)