bayesm (version 2.0-9)

rnmixGibbs: Gibbs Sampler for Normal Mixtures

Description

rnmixGibbs implements a Gibbs Sampler for normal mixtures.

Usage

rnmixGibbs(Data, Prior, Mcmc)

Arguments

Data

list(y)

Prior

list(Mubar,A,nu,V,a,ncomp) (only ncomp required)

Mcmc

list(R,keep) (R required)

Value

• a list containing:
• probdraw
  R/keep x ncomp array of mixture prob draws
• zdraw
  R/keep x nobs array of indicators of mixture comp identity for each obs
• compdraw
  R/keep lists of lists of comp parm draws

concept

• bayes
• MCMC
• normal mixtures
• Gibbs Sampling

Details

Model: $y_i$ $\sim$ $N(mu_{ind_i},Sigma_{ind_i})$. ind $\sim$ iid multinomial(p). p is a ncomp x 1 vector of probs. Priors: $mu_j$ $\sim$ $N(mubar,Sigma_j (x) A^{-1})$. $mubar=vec(Mubar)$. $Sigma_j$ $\sim$ IW(nu,V). note: this is the natural conjugate prior -- a special case of multivariate regression. $p$ $\sim$ Dirchlet(a). Output of the components is in the form of a list of lists. compsdraw[[i]] is ith draw -- list of ncomp lists. compsdraw[[i]][[j]] is list of parms for jth normal component. jcomp=compsdraw[[i]][j]]. Then jth comp $\sim$ $N(jcomp[[1]],Sigma)$, $Sigma$ = t(R)%*%R, $R^{-1}$ = jcomp[[2]]. List arguments contain:

• y

{ n x k array of data (rows are obs) } Mubar{ 1 x k array with prior mean of normal comp means (def: 0)} A{ 1 x 1 precision parameter for prior on mean of normal comp (def: .01)} nu{ d.f. parameter for prior on Sigma (normal comp cov matrix) (def: k+3)} V{ k x k location matrix of IW prior on Sigma (def: nuI)} a{ ncomp x 1 vector of Dirichlet prior parms (def: rep(5,ncomp))} ncomp{ number of normal components to be included } R{ number of MCMC draws } keep{ MCMC thinning parm: keep every keepth draw (def: 1)}

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 3. http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html

See Also

rmixture, rmixGibbs ,eMixMargDen, momMix, mixDen, mixDenBi

Examples

##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

set.seed(66)
dim=5;  k=3   # dimension of simulated data and number of "true" components
sigma = matrix(rep(0.5,dim^2),nrow=dim);diag(sigma)=1
sigfac = c(1,1,1);mufac=c(1,2,3); compsmv=list()
for(i in 1:k) compsmv[[i]] = list(mu=mufac[i]*1:dim,sigma=sigfac[i]*sigma)
comps = list() # change to "rooti" scale
for(i in 1:k) comps[[i]] = list(mu=compsmv[[i]][[1]],rooti=solve(chol(compsmv[[i]][[2]])))
pvec=(1:k)/sum(1:k)

nobs=500
dm = rmixture(nobs,pvec,comps)

Data=list(y=dm$x)
ncomp=9
Prior=list(ncomp=ncomp)
Mcmc=list(R=R,keep=1)
out=rnmixGibbs(Data=Data,Prior=Prior,Mcmc=Mcmc)

tmom=momMix(matrix(pvec,nrow=1),list(comps))
if(R < 1000) {begin=1} else {begin=500}
pmom=momMix(out$probdraw[begin:R,],out$compdraw[begin:R])
mat=rbind(tmom$mu,pmom$mu)
rownames(mat)=c("true","post expect")
cat("mu and posterior expectation of mu",fill=TRUE)
print(mat)
mat=rbind(tmom$sd,pmom$sd)
rownames(mat)=c("true","post expect")
cat("std dev and posterior expectation of sd",fill=TRUE)
print(mat)
mat=rbind(as.vector(tmom$corr),as.vector(pmom$corr))
rownames(mat)=c("true","post expect")
cat("corr and posterior expectation of corr",fill=TRUE)
print(t(mat))

if(0){
##
## plotting examples
##
## check true and estimated marginal densities
end=R
grid=NULL
for (i in 1:dim){
  rgi=range(dm$x[,i])
  gr=seq(from=rgi[1],to=rgi[2],length.out=50)
  grid=cbind(grid,gr)
}

tmden=mixDen(grid,pvec,comps)
pmden=eMixMargDen(grid,out$probdraw[begin:end,],out$compdraw[begin:end])

## plot the marginal on third variable
plot(range(grid[,3]),c(0,1.1*max(tmden[,3],pmden[,3])),type="n",xlab="",ylab="density")
lines(grid[,3],tmden[,3],col="blue",lwd=2)
lines(grid[,3],pmden[,3],col="red",lwd=2)

## compute implied bivariate marginal distributions
i=1
j=2
rxi=range(dm$x[,1])
rxj=range(dm$x[,2])
xi=seq(from=rxi[1],to=rxi[2],length.out=50)
xj=seq(from=rxj[1],to=rxj[2],length.out=50)
den=matrix(0,ncol=length(xi),nrow=length(xj))
for(ind in as.integer(seq(from=begin,to=end,length.out=100))){
  den=den+mixDenBi(i,j,xi,xj,out$probdraw[ind,],out$compdraw[[ind]])
}
tden=matrix(0,ncol=length(xi),nrow=length(xj))
tden=mixDenBi(i,j,xi,xj,pvec,comps)

par(mfrow=c(2,1))
image(xi,xj,tden,col=terrain.colors(100),xlab="",ylab="")
contour(xi,xj,tden,add=TRUE,drawlabels=FALSE)
title("True Bivariate Marginal")
image(xi,xj,den,col=terrain.colors(100),xlab="",ylab="")
contour(xi,xj,den,add=TRUE,drawlabels=FALSE)
title("Posterior Mean of Bivariate Marginal")
}