BayesianNonparametrics.jl is a Julia package implementing state-of-the-art Bayesian nonparametric models for medium-sized unsupervised problems. The software package brings Bayesian nonparametrics to non-specialists allowing the widespread use of Bayesian nonparametric models. Emphasis is put on consistency, performance and ease of use allowing easy access to Bayesian nonparametric models inside Julia.
BayesianNonparametrics.jl allows you to
- explain discrete or continous data using: Dirichlet Process Mixtures or Hierarchical Dirichlet Process Mixtures
- analyse variable dependencies using: Variable Clustering Model
- fit multivariate or univariate distributions for discrete or continous data with conjugate priors
- compute point estimtates of Dirichlet Process Mixtures posterior samples
- julia version 0.5
- packages listed in REQUIREMENTS file
You can clone the package into your running julia 0.5 installation using
Pkg.clone("https://github.com/OFAI/BayesianNonparametrics.jl.git")
Documentation is available on readthedocs: http://bnpjl.readthedocs.org/en/latest/
The following example illustrates the use of BayesianNonparametrics.jl for clustering of continuous observations using a Dirichlet Process Mixture of Gaussians.
After loading the package:
using BayesianNonparametrics
we can generate a 2D synthetic dataset (or use a multivariate continuous dataset of interest)
(X, Y) = bloobs(randomize = false)
and construct the parameters of our base distribution:
μ0 = vec(mean(X, 1))
κ0 = 5.0
ν0 = 9.0
Σ0 = cov(X)
H = WishartGaussian(μ0, κ0, ν0, Σ0)
After defining the base distribution we can specify the model:
model = DPM(H)
which is in this case a Dirichlet Process Mixture. Each model has to be initialised, one possible initialisation approach for Dirichlet Process Mixtures is a k-Means initialisation:
modelBuffer = init(X, model, KMeansInitialisation(k = 10))
The resulting buffer object can now be used to apply posterior inference on the model given
models = train(modelBuffer, DPMHyperparam(), Gibbs(maxiter = 500))
You shoud see the progress of the sampling process in the command line. After applying Gibbs sampling, it is possible explore the posterior based on their posterior densities,
densities = Float64[m.energy for m in models]
number of active components
activeComponents = Int[sum(m.weights .> 0) for m in models]
or the groupings of the observations:
assignments = [m.assignments for m in models]
The following animation illustrates posterior samples obtained by a Dirichlet Process Mixture:
Alternatively, one can compute a point estimate based on the posterior similarity matrix:
A = reduce(hcat, assignments)
(N, D) = size(X)
PSM = ones(N, N)
M = size(A, 2)
for i in 1:N
for j in 1:i-1
PSM[i, j] = sum(A[i,:] .== A[j,:]) / M
PSM[j, i] = PSM[i, j]
end
end
and find the optimal partition which minimizes the lower bound of the variation of information:
mink = minimum([length(m.weights) for m in models])
maxk = maximum([length(m.weights) for m in models])
(peassignments, _) = pointestimate(PSM, method = :average, mink = mink, maxk = maxk)
The grouping wich minimizes the lower bound of the variation of information is illustrated in the following image: