Simple introduction to MCMC with Dirichlet process prior?

Question

I'm looking for a simple and easy to read introduction to using MCMC with a dirichlet process prior. Or perhaps using MCMC in any machine learning scenario, eg Gaussian Process.

I've been circling around various papers and tutorials (Neal, Teh, Sahu, Ghahramani, Ferguson, Escobar and West, ... ), but I can't find anything that explains in straightforward procedural terms, specifying explicitly each equation I would need to implement MCMC for DP or GP.

Edit: I got this far: my mcmc dp progress

Answer 1

Its not a paper, but I found some Matlab code that implements a DP prior for an infinite Gaussian mixture model. The code uses Gibbs sampling to infer a GMM (and the number of components in the mixture) over some input data. The code is pretty readable and it has helped me quite a bit to see the DP in action in a concrete example.

Answer 2

MCMC sampling for DPMM is quite challenging and that's for many reasons, the main one being that the model is infinite and the distribution is not that easy to work with. Employing algorithms such as metropolis is non-trivial since there are actually quite a few degrees of freedom in how you specify your proposal distributions and acceptance ratios.

You can trace part of the evolution of DPMM MCMC algorithms in the work of Radford Neal (University of Toronto) here.

You could start with this 1998 paper and move on to his split and merge for conjugate and non-conjugate models from 2000 and 2005. These will help you understand the main ideas behind MCMC sampling for DPMMs.

When you feel brave you can hop on to the most recent advances. One of them is Chang's and Fisher's (MIT) 2013 paper which you can find here. You can also find code and a demo on Jason Chang's website.

If you need a more general mathematical exploration of DPs you can find many of them online. A good one is by Yee Whye Teh (UCL). I would post more links but my reps are not enough at this point :-)

Answer 3

I have the same feeling. This is as closest as I've come:

http://www.cs.cmu.edu/~kbe/dp_tutorial.pdf

The algorithm explained starting at page 37 is understandable, however I still wish to see a very simple example written out step by step so I'm more confident of which each term means.

I attempted to implement the algorithm in R, below is my code, not efficient and not sure if it is all correct!

#Dirichlet mixture of normals
#flat prior on mu - posterior normal
#flat prior on precision - posterior gamma


###generate data
library(mixtools)
n=100
x=rnormmix(n,lambda=c(.5,.5),mu=c(-4,2),sigma=c(1,1))


###get probs for each component
cProbs=function(x,c,etamu,etavar,alpha,i,n){
  c1=c[-i]
  cProb=rep(NA,length(unique(c1)))

  for(j in 1:length(cProb))
    cProb[j]=length(which(c1==unique(c1)[j]))/(n-1+alpha)*dnorm(x[i],etamu[j],sqrt(etavar[j]))

  return(cProb)  
}

###get probs for new component
newProb=function(x,alpha,i,n){
  s=sqrt(sum((x-mean(x))^2)/(n-1))

  #predictive distribution is t
  newProb1=alpha/(n-1+alpha)*dt((x[i]-mean(x))/(s*sqrt(1+1/n)),n-1)

  return(newProb1)  
}


#parms
alpha=.01
nsim=100


#initialize
etamu=mean(x)
etavar=var(x)
c=rep(1,n)
idx=0

#loop
repeat{
  idx=idx+1
  for(i in 1:n){

    ####if c_i singleton remove
    if(sum(c==c[i])==1){
      etamu=etamu[-c[i]]
      etavar=etavar[-c[i]]
      c[c>c[i]]=c[c>c[i]]-1
      c[i]=1
    }

    ####draw new c_i
    #get probabilities for components
    probsExisting=cProbs(x,c,etamu,etavar,alpha,i,n)
    newP=newProb(x,alpha,i,n)

    #sample
    temp=sample(1:(length(unique(c))+1),size=1,prob=c(probsExisting,newP))
    newC=ifelse(temp==length(unique(c))+1,1,0)
    c[i]=temp

    ####if new c_i draw new eta
    if(newC==1){
      #var: posterior is inverse gamma
      newVar1=1/rgamma(1,(n-1)/2,(n-1)/2*sum((x-mean(x))^2)/(n-1))
      etavar=c(etavar,newVar1)

      #mean: posterior is normal
      newMean1=rnorm(1,x[i],sqrt(newVar1/n))
      etamu=c(etamu,newMean1)

    }
  } #end loop over observations

  #updata etas
  for(i in 1:length(unique(c))){
    n1=sum(c==i)
    if(n1>1){
      temp=x[which(c==i)]
      etavar[i]=1/rgamma(1,(n1-1)/2,(n1-1)/2*sum((temp-mean(temp))^2)/(n1-1)) 
      etamu[i]=rnorm(1,mean(temp),sqrt(etavar[i]/n1))
    }
  } #end update etas

  if(idx==nsim) break
} #end repeat



###plot results
probs=rep(NA,length(etamu))
for(i in 1:length(probs))
  probs[i]=sum(c==i)/n

grid=seq(min(x)-1,max(x)+1,length=500)
dens=rep(NA,length=length(grid))
for(i in 1:length(grid))
  dens[i]=sum(probs*dnorm(grid[i],etamu,sqrt(etavar)))

hist(x,freq=FALSE)
lines(grid,dens,col='red',lwd=2)

Answer 4

My belief is that there can hardly be such a paper, because the concepts and implementation details are not straighforward. Hence, even the best papers will be somewhat involved and you need to go through that.

Answer 5

(Edited) Since, no one had added this. Let me contribute here.

Sethuraman's stick breaking representation for Dirichlet Process is incredibly useful for understanding DP model and also for simulating DP. Basically, Sethuraman tells us that one can think of the weights which appear in Dirichlet Process as a product of stick breaking process (beta). Here's the link for the paper. http://www.jstor.org/stable/24305538?seq=1#page_scan_tab_contents

One of the main reason, why DP is hard to handle because it's literally an infinite mixture model. An excellent paper by Ishwaran tells us that it's alright to approximate DP as a finite mixture model. Here's the paper https://people.eecs.berkeley.edu/~jordan/sail/readings/archive/ishwaran-Mixture.pdf.

Once, you have internalized above methods, it's fairly straightforward to do Gibbs Sampling with DP with say Normal Base Distribution. Let's say that your data y also comes from the Normal Distribution. let's also that the cluster mean of your data (y) comes from DP. Then, here are the steps to do a simple Gibbs Sampling with DP

1. Choose a large K (total number of cluster).
2. Initialize cluster mean and variance for each cluster.
3. Initialize the stick breaking prior. 2#. Choose a data point
4. Compute the probability that data point lies in a particular cluster .
5. Do the Step 3 for all of the cluster.
6. Normalize the above probability, so that it adds up to 1.
7. Using the above probability allocate the data point to a particular cluster.
8. Do step 3-6 for all your data points.
9. Now update the number of points in all the cluster.
10. Using basic Bayesian Methods, update the cluster means and cluster variance.
11. Repeat the above process for however number of iteration you want.
12. Discard the firs half of your sample.
13. Viola, you did a MCMC for your DP model.

Here's an excellent R code shared by Brian Neelon, which uses the above two concept to make an elegant DP code. http://people.musc.edu/~brn200/r/DP03.r