DPMM 的 MCMC 采样非常具有挑战性，原因有很多，主要原因是模型是无限的并且分布不是那么容易使用。使用诸如 Metropolis 之类的算法并非易事，因为在您指定提案分布和接受率方面实际上有相当多的自由度。

您可以在此处追溯 Radford Neal（多伦多大学）的工作中 DPMM MCMC 算法的部分演变。

您可以从 1998 年的这篇论文开始，然后继续他在 2000 年和 2005 年对共轭和非共轭模型的拆分和合并。这些将帮助您了解 DPMM 的 MCMC 采样背后的主要思想。

当您感到勇敢时，您可以了解最新的进展。其中之一是 Chang 和 Fisher 的 (MIT) 2013 年论文，您可以在此处找到。您还可以在 Jason Chang 的网站上找到代码和演示。

如果您需要对 DP 进行更一般的数学探索，您可以在网上找到其中的许多。Yee Whye Teh (UCL) 是一个很好的作品。我会发布更多链接，但此时我的代表还不够:-)

它不是一篇论文，但我发现了一些Matlab 代码，它为无限高斯混合模型实现了 DP 先验。该代码使用 Gibbs 采样来推断一些输入数据的 GMM（以及混合物中的成分数量）。代码非常易读，它帮助我在一个具体的例子中看到了 DP 的运行情况。

我也有同样的感觉。这是我来过的最接近的：

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

从第 37 页开始解释的算法是可以理解的，但是我仍然希望看到一个非常简单的示例逐步写出，因此我对每个术语的含义更有信心。

我试图在 R 中实现该算法，下面是我的代码，效率不高，不确定是否全部正确！

#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)

我认为几乎不可能有这样的论文，因为概念和实现细节并不直截了当。因此，即使是最好的论文也会有所涉及，你需要经历这些。