bayesm (version 3.1-6)

rnmixGibbs: Gibbs Sampler for Normal Mixtures

Description

rnmixGibbs implements a Gibbs Sampler for normal mixtures.

Usage

rnmixGibbs(Data, Prior, Mcmc)

Value

A list containing:

ll

\(R/keep x 1\) vector of log-likelihood values

nmix

a list containing: probdraw, zdraw, compdraw (see “nmix Details” section)

Arguments

Data

list(y)

Prior

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

Mcmc

list(R, keep, nprint, Loglike)

Author

Peter Rossi, Anderson School, UCLA, .

Details

Model and Priors

\(y_i\) \(\sim\) \(N(\mu_{ind_i}, \Sigma_{ind_i})\)
ind \(\sim\) iid multinomial(p) with \(p\) an \(ncomp x 1\) vector of probs

\(\mu_j\) \(\sim\) \(N(mubar, \Sigma_j (x) A^{-1})\) with \(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)

Argument Details

Data = list(y)

y:	\(n x k\) matrix of data (rows are obs)

Prior = list(Mubar, A, nu, V, a, ncomp) [only ncomp required]

Mubar:	\(1 x k\) vector with prior mean of normal component means (def: 0)
A:	\(1 x 1\) precision parameter for prior on mean of normal component (def: 0.01)
nu:	d.f. parameter for prior on Sigma (normal component cov matrix) (def: k+3)
V:	\(k x k\) location matrix of IW prior on Sigma (def: nu*I)
a:	\(ncomp x 1\) vector of Dirichlet prior parameters (def: rep(5,ncomp))
ncomp:	number of normal components to be included

Mcmc = list(R, keep, nprint, Loglike) [only R required]

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keepth draw (def: 1)
nprint:	print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)
LogLike:	logical flag for whether to compute the log-likelihood (def: FALSE)

nmix Details

nmix is a list with 3 components. Several functions in the bayesm package that involve a Dirichlet Process or mixture-of-normals return nmix. Across these functions, a common structure is used for nmix in order to utilize generic summary and plotting functions.

probdraw:	\(ncomp x R/keep\) matrix that reports the probability that each draw came from a particular component
zdraw:	\(R/keep x nobs\) matrix that indicates which component each draw is assigned to
compdraw:	A list of \(R/keep\) lists of \(ncomp\) lists. Each of the inner-most lists has 2 elemens: a vector of draws for mu and a matrix of draws for the Cholesky root of Sigma.

References

For further discussion, see Chapter 3, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

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)

# change to "rooti" scale
comps = list() 
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$nmix)

tmom = momMix(matrix(pvec,nrow=1), list(comps))
mat = rbind(tmom$mu, tmom$sd)
cat(" True Mean/Std Dev", fill=TRUE)
print(mat)

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