added readme, DPM code and point estimate code.

Author: trappmartin

README.md

@@ -1,4 +1,76 @@
-# BayesianNonparametrics
+# BayesianNonparametrics.jl
+
+## Examples
+### Dirichlet Process Mixture
+The following example illustrates the use of BayesianNonparametrics.jl for clustering of continuous observations. 
+
+After loading the package:
+```julia
+using BayesianNonparametrics
+```
+we can generate a 2D synthetic dataset (or use a multivariate continuous dataset of interest)
+```julia
+(X, Y) = bloobs(randomize = false)
+```
+and construct the parameters of our base distribution:
+```julia
+μ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:
+```julia
+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:
+```julia
+modelBuffer = init(X, model, KMeansInitialisation(k = 10))
+```
+The resulting buffer object can now be used to apply posterior inference on the model given $X$. In the following we apply Gibbs sampling for 500 iterations without burn in or thining:
+```julia
+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,
+```julia
+densities = Float64[m.energy for m in models]
+```
+number of active components
+```julia
+activeComponents = Int[sum(m.weights .> 0) for m in models]
+```
+or the groupings of the observations:
+```julia
+assignments = [m.assignments for m in models]
+```
+The following animation illustrates posterior samples obtained by a Dirichlet Process Mixture: 
+![alt text](posteriorSamples.gif "Posterior Sample")
+
+Alternatively, one can compute a point estimate based on the posterior similarity matrix:
+```julia
+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:
+```julia
+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:
+![alt text](pointestimate.png "Point Estimate")
 
 [![Build Status](https://travis-ci.org/trappmartin/BayesianNonparametrics.jl.svg?branch=master)](https://travis-ci.org/trappmartin/BayesianNonparametrics.jl)
 

REQUIRE

@@ -1 +1,3 @@
 julia 0.5
+Distributions
+Combinatorics

demo/DirichletProcessMixture.jl

@@ -0,0 +1,37 @@
+using BayesianNonparametrics
+
+(X, Y) = bloobs(randomize = false)
+
+μ0 = vec(mean(X, 1))
+κ0 = 5.0
+ν0 = 9.0
+Σ0 = cov(X)
+
+H = WishartGaussian(μ0, κ0, ν0, Σ0)
+
+model = DPM(H)
+initialisation = KMeansInitialisation(k = 10)
+
+modelBuffer = init(X, model, initialisation)
+models = train(modelBuffer, DPMHyperparam(), Gibbs(maxiter = 500))
+
+densities = Float64[m.energy for m in models]
+activeComponents = Int[sum(m.weights .> 0) for m in models]
+assignments = [m.assignments for m in models]
+
+A = reduce(hcat, assignments)
+(N, D) = size(X)
+
+PSM = zeros(N, N)
+M = size(A, 2)
+for i in 1:N
+  for j in 1:i
+    PSM[i, j] = sum(A[i,:] .== A[j,:]) / M
+    PSM[j, i] = sum(A[i,:] .== A[j,:]) / M
+  end
+end
+
+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)

pointestimate.png

@@ -1,5 +1,17 @@
 module BayesianNonparametrics
 
-# package code goes here
+  using Distributions, Combinatorics, Clustering, ProgressMeter
+
+  include("common.jl")
+  include("math.jl")
+  include("utils.jl")
+  include("datasets.jl")
+  include("distributions.jl")
+  include("distfunctions.jl")
+  include("inits.jl")
+  include("inference.jl")
+  include("models.jl")
+  include("dpmm.jl")
+  include("vipointestimate.jl")
 
 end # module

posteriorSamples.gif

@@ -0,0 +1,44 @@
+export init, train
+
+"abstract Hyperparameters"
+abstract AbstractHyperparam;
+
+"Abstract Model Data Object"
+abstract AbstractModelData;
+
+"Abstract Model Buffer Object"
+abstract AbstractModelBuffer;
+
+abstract ModelType;
+
+abstract InitialisationType;
+
+abstract PosteriorInference;
+
+function init(X, model::ModelType, init::InitialisationType)
+  throw(ErrorException("Initialisation $(init) for $(model) is not available."))
+end
+
+function extractpointestimate(B::AbstractModelBuffer)
+  throw(ErrorException("No point estimate available for $(typeof(B))."))
+end
+
+function sampleparameters!(B::AbstractModelBuffer, P::AbstractHyperparam)
+  throw(ErrorException("Parameter sampling for $(typeof(B)) using $(typeof(P)) is not available."))
+end
+
+function gibbs!(B::AbstractModelBuffer)
+  throw(ErrorException("Gibbs sampling for $(typeof(B)) is not available."))
+end
+
+function slicesample!(B::AbstractModelBuffer)
+  throw(ErrorException("Slice sampling for $(typeof(B)) is not available."))
+end
+
+function variationalbayes!(B::AbstractModelBuffer)
+  throw(ErrorException("Variational inference for $(typeof(B)) is not available."))
+end
+
+function train(B::AbstractModelBuffer, P::AbstractHyperparam, I::PosteriorInference)
+  throw(ErrorException("Posterior inference for $(typeof(B)) is not available."))
+end

src/BayesianNonparametrics.jl

@@ -0,0 +1,26 @@
+export bloobs
+
+function bloobs(;centers = 3, samples = 100, randomize = true)
+
+  μ = reduce(hcat, [rand(Uniform(-10, 10), centers) for i in 1:2])
+  Σ = reduce(hcat, [rand(Uniform(0.5, 4), centers) for i in 1:2])
+
+  samplespcenter = ones(Int, centers) * round(Int, samples / centers)
+
+  for i = 1:(centers % samples)
+      samplespcenter[i] += 1
+  end
+
+  X = reduce(hcat, [rand(MvNormal(ones(2) .* μ[i], eye(2) .* Σ[i]), samplespcenter[i]) for i in 1:centers])'
+  Y = reduce(vcat, [ones(Int, samplespcenter[i]) * i for i in 1:centers])
+
+  ids = collect(1:size(X, 1))
+  if randomize
+    shuffle!(ids)
+  end
+
+  X = X[ids,:]
+  Y = Y[ids]
+
+  return (X, Y)
+end