bayesm (version 3.0-2)

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,nprint,Loglike) (R required)

Value

nmix

a list containing: probdraw,zdraw,compdraw

ll

vector of log-likelihood values

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)

• nprint print the estimated time remaining for every nprint'th draw (def: 100)

• LogLike logical flag for compute log-likelihood (def: FALSE)

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 3. http://www.perossi.org/home/bsm-1

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)

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

cat("Summary of Normal Mixture Distribution",fill=TRUE)
summary(out)
tmom=momMix(matrix(pvec,nrow=1),list(comps))
mat=rbind(tmom$mu,tmom$sd)
cat(" True Mean/Std Dev",fill=TRUE)
print(mat)

if(0){
##
## plotting examples
##
plot(out$nmix,Data=dm$x)
}