gibbs: Gibbs sampler and Chib's evidence approximation for a generic...

In bayess: Bayesian Essentials with R

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)

bayess documentation built on May 29, 2024, 3:58 a.m.