Source-linked AI summary
dirichletprocess: An R Package for Fitting Complex Bayesian Nonparametric Models
Gordon J. Ross, Dean Markwick, Priyanshu Tiwari
TL;DR
Parametric assumptions can be unjustified when prior information and data are insufficient, while DPMM inference is difficult to implement. The dirichletprocess package provides flexible R objects, built-in models, and sampling algorithms so DPMMs can serve as components in user-specified models and broader statistical workflows.
Problem
Parametric distributional assumptions may lack sufficient prior information and data, while DPMM posterior inference typically requires substantial statistical and programming expertise.
Method
The package represents DPMMs as extensible R objects, provides several mixture models and posterior sampling algorithms, and handles Markov chain Monte Carlo inference.
Results
The package enables DPMMs to function as building blocks inside user-specified models rather than restricting users to a fixed collection of tasks.
Takeaways & Limitations
Users can apply DPMMs to density estimation, clustering, and hierarchical models while retaining flexibility to construct broader statistical workflows.
Takeaways & Limitations
Posterior-draw functions conditional on the current fitted state should not be interpreted as MCMC-averaged summaries.
Abstract
from arXiv · showhide
The dirichletprocess package provides software for creating flexible Dirichlet process objects. Users can perform nonparametric Bayesian analysis using Dirichlet processes without the need to program their own inference algorithms. Instead, the user can utilise our pre-built models or specify their own models whilst allowing the dirichletprocess package to handle the Markov chain Monte Carlo sampling. Our Dirichlet process objects can act as building blocks for a variety of statistical models including: density estimation, clustering and prior distributions in hierarchical models.
1. Introduction
Dirichlet process mixture models provide flexible alternatives to restrictive parametric assumptions, but their computational complexity has limited adoption. The dirichletprocess package addresses this by providing extensible R objects, pre-built mixture models, and posterior sampling tools for user-specified statistical workflows.
- Motivation: Nonparametric models let distributions learn essential features from data when parametric assumptions lack sufficient justification.
- Dirichlet process mixtures: DPMMs represent unknown densities with an infinite-dimensional mixture whose mixing distribution has a Dirichlet process prior.The mixture kernel may be Gaussian, and the model can approximate any continuous distribution arbitrarily accurately.
- Hierarchical models: DPMMs can serve as nonparametric prior specifications within hierarchical models, including models that share information across group-level parameters.Using a DPMM for group means avoids committing to a potentially inappropriate parametric prior when group means are measured imprecisely.
- Computational challenge: Analytically intractable DPMM posteriors require simulation algorithms, creating statistical and programming demands that have hindered wider adoption.Gibbs sampling and Metropolis-Hastings are among the methods used to draw posterior samples.
- Package design: The package represents DPMMs as R objects that can be used as building blocks inside user-specified models rather than automating only a fixed task range.Its emphasis is flexibility and extensibility for models that use DPMMs within broader estimation schemes.
- Package features: The package includes Gaussian, Beta, multivariate Normal, and Weibull mixture kernels, conjugate and nonconjugate posterior sampling, and an object-based R interface.Core sampling routines for the main mixture models are implemented in compiled code while R handles object representation and user-facing workflows.
2. Background information
Dirichlet process models provide flexible priors that avoid committing to a fixed parametric family, while DPMMs represent continuous densities through mixtures and support clustering. The package implements posterior sampling and lets users choose kernels and base measures for broader model construction.
- Motivation: Nonparametric priors let distributions adapt to data when available information does not justify a fixed parametric family.This avoids restricting inference to families such as the Normal or Exponential.
- Dirichlet process foundations: A Dirichlet process prior has base measure G0 and concentration parameter α, with α controlling the balance between prior information and observed data.The posterior combines G0 with the empirical distribution, weighted according to α.
- Dirichlet process mixtures: Dirichlet process draws are discrete, so DPMMs place the DP prior on mixture parameters that generate observations through a typically continuous kernel.Repeated parameter values induce mixture components and allow the model to represent clusters.
- Package scope: The package returns posterior samples of both mixture parameters and F, supporting clustering, density-focused inference, and other downstream uses.Users can specify their own kernel k and base measure G0, while implemented choices also provide selected hyperparameter updates.
- Inference: Posterior inference repeatedly samples mixture parameters with the Chinese Restaurant Process sampler, integrating out the distribution F.Conjugate models permit analytic updates, whereas non-conjugate cases require numerical techniques such as Algorithm 8.
3. Package Overview
The package represents Dirichlet process mixture models as flexible R objects that support standard fitting, posterior summaries, clustering, and hierarchical-model building. It provides pre-built kernels while allowing users to control sampling or extend models with new kernels.
- Core design: The package represents Dirichlet process mixture models as abstract R objects that serve as building blocks within hierarchical models.Its emphasis is on flexibility and extensibility rather than automating only fixed tasks such as density estimation or regression.
- Inference: Fit provides out-of-the-box MCMC fitting by updating cluster allocations, cluster parameters, and the concentration parameter α after the mixture kernel is specified.Users can adjust sampling behavior, including how often α is updated, when default updating may delay convergence.
- Core design: It includes Gaussian, multivariate Gaussian, Beta, and Weibull mixture-kernel implementations for nonparametric Bayesian analysis.The package also supports marginal and truncated stick-breaking representations of Dirichlet process samples.
- Extensibility: The package supports custom Dirichlet process models by reusing its resampling functions, illustrated with Poisson and Gamma kernels not currently implemented in the package.It also provides plotting for fitted densities, posterior summaries, and pointwise credible intervals.
- Examples: Examples demonstrate multimodal density estimation, cluster recovery, hierarchical shared components, and recovery of an intensity shape when the full dataset is unobserved.The examples include Old Faithful density estimation, faithful-data clustering, hierarchical group densities, common clusters across datasets, and a point-process application.
- Inference: Users can manually control resampling and inspect the current Dirichlet process state, enabling DP objects to participate inside user-specified hierarchical models.Manual update functions provide control over the MCMC routine, while retained samples must be stored explicitly rather than being appended automatically.
4. Advanced Features
Advanced users can customize Dirichlet process mixture objects by defining kernel, prior, posterior, predictive, or Metropolis–Hastings components. The package then fits these models through its conjugate or nonconjugate Chinese Restaurant Process samplers.
- When customization is needed: The advanced controls support changing base-measure hyperparameters, addressing convergence concerns, or using mixture kernels beyond the defaults.These controls are primarily needed when default functions and specifications are inadequate for a particular application.
- Custom mixture specifications: The package represents mixture models through S3 objects whose kernel specification includes likelihood, prior draws, and base-measure parameters.Conjugate models additionally provide posterior-draw and predictive functions.
- Conjugate models: Conjugate mixture models use direct posterior sampling and analytically evaluable marginal distributions within Neal’s conjugate Algorithm 4.The Fit function performs the Chinese Restaurant Process sampler after these functions are specified.
- Nonconjugate models: Nonconjugate models use Metropolis–Hastings proposals, requiring prior-density, parameter-proposal, and step-size specifications.The resulting model is fitted with Neal’s Algorithm 8 through Fit.
- Object construction: Users can create a dirichletprocess object from the data, mixing distribution, and concentration parameter α, after which the package’s class functions become available.Default Gaussian, Beta, multivariate normal, and Weibull constructors provide non-informative base-measure parameters.
4.2. Creating New Dirichlet process Mixture Types
The package is designed to be extended with new Dirichlet process mixture types rather than limiting users to its built-in kernels. New S3 objects that follow the package specification can reuse existing resampling and prediction functions.
- Motivation: Built-in kernels cannot cover every desired application, so the package supports user-defined mixture models beyond Gaussian, Beta, and Weibull kernels.The extension mechanism is intended for complex models outside the package’s implemented scope.
- S3 extension mechanism: A new mixture type is created as an S3 object that encapsulates the model and conforms to the package’s required specifications.Conforming objects allow the package’s resampling and prediction functions to continue operating on the new type.
- Worked example: The package documents construction of a new S3 type through an extended example so users can implement their own mixture models.The example explains the object construction process in detail.
Conjugate Mixture
A conjugate Poisson mixture can be added by defining its likelihood, Gamma prior draws, posterior sampling, and predictive distribution. The resulting custom model is fitted with the package’s standard tools and closely approximates its generating distribution.
- Model specification: A Poisson mixture uses a Gamma conjugate prior, with a one-element parameter vector because the Poisson distribution has one parameter.The mixing distribution stores the prior parameters α0 and β0.
- Model specification: The custom Poisson type defines a likelihood function that evaluates the Poisson density for the single parameter θ.Parameters remain internally represented as lists, making θ[[1]] necessary even for one-element parameters.
- Conjugate inference: PriorDraw samples Poisson parameters from the Gamma base measure, while PosteriorDraw directly samples from the conjugate posterior.The posterior draw uses updated prior parameters based on the observed data.
- Conjugate inference: The Predictive function evaluates the marginal data distribution analytically for the conjugate Poisson mixture.This predictive function is translated into the package’s expected R interface.
- Fitting and validation: After construction with MixingDistribution, the Poisson model is used with DirichletProcessCreate, Initialise, Fit, and PosteriorSummary.The demonstration generates observations from a two-component Poisson mixture and retains MCMC samples for posterior summaries.
- Fitting and validation: The fitted Poisson mixture closely approximates the true generating distribution, showing that custom mixtures retain the package’s fitting and summary tools.Figure 8 compares the true distribution with the estimated posterior mean and pointwise credible intervals.
Nonconjugate Mixture
A nonconjugate Gamma mixture requires Metropolis–Hastings components because its posterior is not directly tractable. The package accommodates these additions and produces a fitted mixture that closely approximates the generating distribution.
- Model specification: A Gamma mixture is nonconjugate, so its required functions must support Metropolis–Hastings posterior sampling.The Gamma distribution has two parameters, α and β.
- Model specification: The custom model defines a Gamma likelihood and draws both parameters from exponential prior distributions.The likelihood uses the two components of θ as the Gamma parameters.
- Metropolis–Hastings inference: Metropolis–Hastings additionally requires a prior-density function and constrained parameter proposals for the positive parameters α and β.The proposal uses absolute values after normal perturbations to maintain positivity.
- Metropolis–Hastings inference: MixingDistribution stores the prior parameters and proposal scale h, with the model marked as nonconjugate before fitting.The Gamma example uses prior parameters c(0.1, 0.1) and proposal steps c(0.1, 0.1).
- Fitting and validation: The nonconjugate object is initialized and fitted through the same general Dirichlet process workflow used for custom mixtures.The example applies DirichletProcessCreate, Initialise, Fit, and PosteriorSummary to simulated Gamma data.
- Fitting and validation: The fitted Gamma mixture closely approximates the true generating distribution.Figure 9 reports this comparison for the simulated example.
4.3. Extended Example - Working with Censored Observations
The package can be extended to analyze censored survival data by combining a custom censored-Weibull likelihood with existing Dirichlet process mixture functionality. Fitting two treatment datasets produces chain-averaged density and survival estimates consistent with a prior qualitative comparison.
- Data and model: A Dirichlet process mixture of Weibull distributions provides a nonparametric approach for modeling censored survival data.The package does not directly support censoring, but its extensible design allows users to add this functionality.
- Data and model: Censored survival observations record event times with an indicator for right censoring, requiring adapted likelihoods.The example concerns leukaemia remission times from two treatment groups.
- Custom likelihood: The custom likelihood translates uncensored observations into Weibull densities and censored observations into survival-probability contributions.The implementation validates positive Weibull parameters and returns the resulting likelihood values for sampling.
- Custom likelihood: The custom weibullcens class reuses Weibull and non-conjugate mixture behavior, so the full fitting machinery does not need reimplementation.Method dispatch applies the censored-Weibull likelihood while retaining general package functionality.
- Fitting and results: 500 MCMC iterations were fitted for each treatment, after rescaling remission times and initializing separate censored-Weibull mixture objects.The fitted objects yield chain-averaged density and survival estimates for comparison.
- Fitting and results: The resulting density and survival estimates were consistent with the qualitative treatment comparison reported by Kottas (2006).Figure 10 presents the chain-averaged point estimates for both treatments.
4.4. Resampling Component Indexes and Parameters
When fitting a Dirichlet process object, the package resamples component assignments and cluster parameters using different procedures for conjugate and non-conjugate mixtures.
- Resampling procedures: Fit resamples component indexes and parameters using Neal’s Algorithm 4 for conjugate models and Algorithm 8 for non-conjugate models.ClusterComponentUpdate and ClusterParameterUpdate perform most of the computation.
- Conjugate mixtures: Conjugate mixtures draw new component indexes and cluster parameters directly from predictive and posterior distributions, making fitting efficient.Users can modify the base-distribution prior parameters to support the underlying data.
- Non-conjugate mixtures: Non-conjugate mixtures propose new components from the prior and sample cluster parameters with Metropolis-Hastings under the parameters’ support constraints.For Weibull parameters, the random walk is restricted to positive values.
Overriding Default Behaviour
Users can override the package’s default resampling behavior while retaining its other S3 methods and object structures.
- Custom sampling: Custom ClusterComponentUpdate and ClusterParameterUpdate functions can replace the default sampling schemes for conjugate and non-conjugate mixtures.The remaining package infrastructure continues to be available to the customized mixture model.
- Auxiliary variables: For non-conjugate mixtures, users can change m, the number of auxiliary variables drawn for each datapoint during initialization.The default is m = 3; increasing the number can lead to more changes in component indexes.
4.5. Resampling the Base Measure, G0
The package can update the base measure G0 during fitting when its parameters are uncertain or difficult to set appropriately in advance.
- Updating G0: Hyper-prior parameters on G0 allow the base measure to be updated at each fitting iteration from the current cluster parameters.This provides an alternative when manually selecting suitable base-measure parameters is impractical.
- Updating G0: Users enable base-measure resampling with Fit(dp,...,updatePrior=TRUE) for mixture distributions that support it.The package specifies hyper-prior distributions for the implemented mixture kernels and permits changing their default hyperparameter priors.
4.6. Component Prediction
ClusterLabelPredict estimates cluster allocations for new data from a fitted Dirichlet process, returning predicted labels and associated cluster information. In the faithful dataset example, the withheld test points were correctly assigned to their appropriate clusters.
- Prediction method: The faithful example withholds the last five observations, fits a multivariate normal Dirichlet process to the training data, and predicts their labels.The fitted object uses 1,000 iterations before prediction.
- Prediction method: ClusterLabelPredict estimates new-data cluster allocations using probabilities for existing clusters or a new cluster, then samples a label.The calculation uses an algorithm appropriate to conjugate or non-conjugate mixture models.
- Returned outputs: The function returns predicted labels, assigned cluster parameters, points per cluster, and the total number of clusters.The number of clusters can change with each prediction.
- Example result: The returned prediction results can be assembled into a dataframe and plotted against the training data.The example plots eruption duration against waiting time and distinguishes predicted test points with larger triangle markers.
- Example result: The test data were correctly identified with the appropriate cluster in the faithful dataset visualization.Training points and predicted test points are plotted with cluster colours, while test points use triangle markers.