dpmixsim: Dirichlet Process mixture model for clustering and image...

In dpmixsim: Dirichlet Process Mixture model simulation for clustering and image segmentation

Description

dpmixsim implements a Dirichlet Process mixture (DPM) model. The DPM model is a Bayesian nonparametric methodology that relies on MCMC simulations for exploring mixture models with an unknown number of components. The function implements conjugate models with normal structure (conjugate normal-normal DP mixture model).

Usage

dpmixsim(x, M=1, a=1, b=1, upalpha=1, a0=2, b0=2, maxiter=4000, rec=3000,
  fsave=NA, kmax=30, nclinit=NA, minvar=0.001)

Arguments

x	scaled input data as vector in range {0,1}
M	DP precision hyperparameter
a	Gamma prior hyperparameter
b	Gamma prior hyperparameter
upalpha	is a logical variable for simulations with {automatic,fixed} calibration of the precision hyperparameter M (default = TRUE)
a0	Gamma prior hyperparameter for M (default 2)
b0	Gamma prior hyperparameter for M (default 2)
maxiter	maximum number of MCMC iteration steps
rec	record the last rec iteration steps
fsave	filename for saving the MCMC simulation (def: NULL do not save)
kmax	maximum number of clusters in the simulation, (default 30)
nclinit	number of initial clusters to use at the beginning of the simulation. If not specified (NA) the number of initial clusters is equal to the length of x (one element per cluster); (default: NA)
minvar	minimum value admissible for a cluster variance (default=0.001). Decreasing minval may improve resolution (distribution fitness), but increases the maximum number of admissible clusters (kmax). In this case, you may have to increase (kmax) as well.

Details

Consider n observations x_1,...,x_n which we regard as exchangeable. We model the distribution from which the x_i are drawn as a mixture of distributions. Dirichlet process mixture models are based on Dirichlet process priors for the primary parameters θ_i. DP mixture models assume that the prior distribution function G itself is uncertain, drawn from a Dirichlet process G ~ DP(M G_0), with base prior G_0 and precision parameter M. This specification may be expressed by the hierarchical model:

x_i ~ N(.|θ_i, σ^2)
θ_i ~ G
G ~ DP(M N(0,1))
σ^{-2} ~ Gamma(a,b)

Value

simulation output as a list of draws containing:

krec	cluster indicator variables
wrec	cluster weights
phirec	theta cluster parameters
varrec	sigma cluster parameters

Author(s)

Adelino Ferreira da Silva, Universidade Nova de Lisboa, Faculdade de Ciencias e Tecnologia, Portugal, [e-mail].

References

Adelino Ferreira da Silva, A Dirichlet process mixture model for brain MRI tissue classification, Medical Image Analysis 11 (2007) 169-182.

Adelino Ferreira da Silva, Bayesian mixture models of variable dimension for image segmentation, Comput. Methods Programs Biomed. 94 (2009) 1-14.

See Also

readsliceimg, postdataseg, postdpmixciz, postimgclgrp, postimgcomps, postkcluster, premask, readsliceimg

Examples

## Not run: 
## Example 1: simple test using `galaxy' data
  data("galaxy")
  x0 <- galaxy$speed
  x  <- prescale(x0) 
  maxiter <- 4000; rec <- 3000; ngrid <- 100
  res <- dpmixsim(x, M=1, a=1, b=0.1, upalpha=1, maxiter=maxiter, rec=rec,
    nclinit=4)
  z <-  postdpmixciz(x=x, res=res, rec=rec, ngrid=ngrid, plot=T)
  ##
  res <- dpmixsim(x, M=2, a=1, b=0.001, upalpha=0, maxiter=maxiter,
    rec=rec, nclinit=4)
  z <-  postdpmixciz(x, res=res, rec=rec, ngrid=ngrid, plot=T)
##-----------------
## Example 2: 
  demo(testMarronWand)
##-----------------
## Example 3: MRI segmentation
## Testing note: this example should reproduce the equivalent segmented
## images used in the author's references 
  slicedata <- readsliceimg(fbase="t1_pn3_rf0", swap=FALSE)
  image(slicedata$niislice, col=gray((0:255)/256), main="original image")
  x0 <- premask(slicedata, subsamp=TRUE)
  x  <- prescale(x0) 
  rec <- 3000
  res <- dpmixsim(x, M=1, a=1, b=1, upalpha=1, maxiter=4000,
      rec=rec, nclinit=8, minvar=0.002)
  ## post-simulation
  ngrid <- 200
  z <- postdpmixciz(x, res=res, rec=rec, ngrid=ngrid, plot=TRUE)
  x0 <- premask(slicedata, subsamp=FALSE) # use full-sized image after estimation 
  x  <- prescale(x0) 
  cx   <- postdataseg(x, z, ngrid=ngrid)
  cat("*** view grouped segmentations:\n")
  postimgclgrp(slicedata$mask, cx, palcolor=FALSE)
  cat("*** display all clusters:\n")
  postimgcomps(slicedata$mask, cx)
  cat("*** re-cluster with 4 clusters:\n")
  postkcluster(slicedata$mask, cx, clk=4)

## End(Not run)

dpmixsim documentation built on May 2, 2019, 5:45 p.m.