Description Usage Arguments Details Value Author(s) References See Also Examples
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).
1 2 |
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. |
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)
simulation output as a list of draws containing:
krec |
cluster indicator variables |
wrec |
cluster weights |
phirec |
theta cluster parameters |
varrec |
sigma cluster parameters |
Adelino Ferreira da Silva, Universidade Nova de Lisboa, Faculdade de Ciencias e Tecnologia, Portugal, afs@fct.unl.pt.
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.
readsliceimg
,
postdataseg
,
postdpmixciz
,
postimgclgrp
,
postimgcomps
,
postkcluster
,
premask
,
readsliceimg
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | ## 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.