跳到主要内容
访问密钥 NCBI主页 MyNCBI主页 主要内容 主导航
数学心理学杂志。作者手稿；PMC 2020年8月1日提供。
以最终编辑形式发布为：
数学心理学杂志。2019年8月；91: 128–144.
2019年5月20日在线发布。 数字对象标识：2016年10月10日/j.jmp.2019.04.04
预防性维修识别码：第6583910页
NIHMSID公司：美国国立卫生研究院1528335
PMID：31217637

Dirichlet过程混合建模教程

李月林,a、，b、，* 伊丽莎白·斯科菲尔德,一和米塔·戈宁b条
作者信息 版权和许可信息 PMC免责声明

摘要

贝叶斯非参数（BNP）模型作为认知的理论模型和分析工具，在心理学中越来越重要。然而，现有的教程往往处于非技术人员无法理解的抽象级别。本教程旨在帮助初学者通过仔细而明确地研究重要但经常被忽略的推导来理解关键概念，重点是将数学与Dirichlet过程混合模型（DPMM）的实际计算解决方案联系起来，该模型是最广泛使用的BNP方法之一。通过研究产生抽象概念的理论，抽象概念向非技术性读者明确而具体。用统计语言R编写的一个可公开访问的计算机程序逐行解释，以帮助读者理解计算算法。该算法还与本期刊中一个可访问的教程中称为“中餐厅流程”的构造方法相关联(Gershman&Blei，2012年)。总的目标是帮助读者更全面地理解理论和应用，以便他们可以在自己的工作中应用BNP方法，并利用本教程中的技术细节来开发新的方法。

关键词：贝叶斯非参数，Dirichlet过程，混合模型，Gibbs抽样，中餐厅过程

1 介绍

贝叶斯非参数（BNP）方法以与传统方法根本不同的方式解决模型复杂性。例如，在传统的聚类分析（如k-means）中，分析员拟合了几种不同复杂度的模型(k个= 1, 2, … ,K（K）集群），然后根据模型比较指标选择所需的模型。在传统的分析中，参数的数量总是固定的。然而，BNP方法基于一个统计框架，该框架原则上包含无穷多个参数。一个这样的统计框架是一个称为Dirichlet过程（DP）的随机过程。本教程旨在向好奇的非技术人员解释DP。

我们主要为那些想学习BNP但发现大多数学术论文抽象程度超出他或她的现有能力的学生撰写。文献中通常用抽象概念描述DP，作为一个随机过程，它将Dirichlet分布从固定类别的共轭先验推广到无限多类别的先验。术语非参数的在这种广泛的背景下，基本上意味着模型原则上有无限多的参数。Austerweil等人（2015）描述此定义。在此框架中定义的BNP模型不限于特定的参数类。对允许的模型有最低限度的假设。模型参数化中这种灵活性的一个突出特点是，随着数据的增加，模型的结构复杂性可能会增加。DP的特征还包括分布上的分布（例如。，2017年;Jara等人，2011年;Jara，2017年;Escobar&West，1995年)。阅读这些材料的学生可能会很好地吸收所有这些概念，但他或她可能对这些略显抽象的概念的起源只有模糊的印象。本教程旨在使这些概念更加具体、明确和透明。

以BNP方法的应用为重点的文章提供了有限的技术细节（例如。，Navarro等人，2006年;卡拉巴索斯，2017年;卡拉巴索斯和沃克，2009年)。更多技术论文，例如最近的几个来源(Orbanz&Teh，2011年;Jara，2017年;Müller等人，2015年;Hjort等人，2012年;Ghoshal&van der Vaart，2015年)解释DP by Measure理论，非技术人员可能不熟悉这个概念。经典论文（例如。，弗格森，1973年;Blackwell&MacQueen，1973年;安东尼亚克，1974年)以及最近的概述(Jara，2017年)也是抽象的。为了最小化这些抽象概念，DP通常使用以下几个施工方法，隐喻有助于解释计算是如何在更直观的层面上实现的(Gershman&Blei，2012年)。然而，重要的细节和推导通常不包括在内。我们认为，有必要弥合这些方法之间的差距，以帮助非技术读者更全面地了解DP是如何在自足的操作指南中实际实现的。

涉及BNP方法的应用迅速增长，包括人类认知(安德森，1991年;奥斯特维尔和格里菲斯，2013年;柯林斯和弗兰克出版社，2013年)，建模个人可变性(Navarro等人，2006年;Liew等人，2016年)和心理测量学(卡拉巴索斯和沃克，2009年,2008)。DP是支持BNP的最简单、使用最广泛的方法之一。已经有很多教程和操作指南（例如。，Gershman&Blei，2012年;Houpt，2018年;奥尔班兹，2018;库鲁克利德斯，2018年）。本教程提供了一些不同的内容。我们的目标是让非技术人员更容易获得技术细节。因此，与基本介绍中通常可用的内容相比，对细节的描述更为明确。我们旨在详细描述DP集群分析如何生成新集群。DP是否以某种方式计算了观测值与其簇之间的拟合度？省略的技术细节是本教程的重点，在使用DP混合建模的示例性聚类分析中进行了介绍。当我们继续学习本教程时，我们将插入对省略推导的解释。目标读者是定量心理学或相关领域的高级研究生，熟悉概率论和条件概率，但不熟悉度量理论。需要熟悉R编程和贝叶斯统计（例如，在Lee，2012年)。然而，所需的编程技巧并不比用R编写简单函数和从标准概率分布中进行随机抽样更复杂。

本教程的组织方式如下。为了建立一个基本框架，我们从一个简单的有限高斯混合模型开始K（K）固定且已知。然后我们详细解释了DP如何将其扩展到无限混合模型。接下来，解释称为中餐厅流程（CRP）的施工方法(奥尔德斯，1985年;Blackwell&MacQueen出版社，1973年)。然后，我们介绍了一个用统计语言R编写的计算机程序，以演示DP在数据积累时产生新集群时的“外观”。R程序改编自教程布罗德里克（2017）中吉布斯采样算法的实现尼尔（2000），算法3）。我们的论述也可能对技术专家有用，因为它可以更新尼尔（2000）以及其他学术论文。

2 有限混合模型

2.1. 假设数据

中的假设数据集图1将在本教程中使用，以使解释更加明确和具体。每个观察由两个评估组成，均按z刻度标准化。每个点代表一名参与者的两个测量值，例如焦虑(年1)和抑郁症状(年2)分别是。观察结果是通过以下方式进行模拟的K（K）=4个已知簇，每个簇由具有平均值的双变量高斯分布表示μk个和协方差矩阵σk个2。所有4个集群的假定已知参数如所示图1.DPMM有望在没有明确告知的情况下找到一个4集群解决方案。

保存图片、插图等的外部文件。对象名为nihms-1528335-f0001.jpg

一个假设的心理评估分数示例，来自分为四个已知组的240个人。从每个聚类中，从具有指定平均值和协方差矩阵的二元高斯分布中随机抽取60个个体。还绘制了基本分布的第95个百分位椭球。

2.2. 有限混合模型中的模型似然函数

我们从一个简单的有限混合模型开始，该模型包含簇数K（K）该模型的似然定义为高斯分布的混合。

第页(年我∣μ1,…,μk个,σ12,…,σk个2,π1,…,πk个)=∑k个=1K（K）πk个N个(年我;μk个,σk个2),
(1)

其中有一个观察结果年我根据具有指定平均值的高斯分布进行评估(μk个)和协方差(σk个2)然后根据混合比例进行加权πk个（分配给簇的权重之和必须为1）。

接下来，我们引入一个指示符变量c（c）我表示第i个参与者的集群成员，其中c（c）我取离散值1、，K（K），所以c（c）我=k个表示我将第个观测值分配给k个第个集群。诸如c（c）我混合模型中是潜在变量。它们的具体值是混合概率的随机实现πk个（想想铸造模具）。原因稍后会变得更清楚，c（c）我之所以在这里介绍，是因为它在简化计算方面起着重要作用。如果我们知道每个观测值属于哪个簇，则根据其簇平均值和协方差计算每个人的数据可能性：

第页(年我∣c（c）我=k个)=N个(年我∣μk个,σk个2).

2.3. 基于贝叶斯定理的联合后验分布

然后应用贝叶斯定理获得模型参数的联合后验分布：

第页(μ,σ2,c（c）∣年)=第页(年∣μ,σ2,c（c）)第页(μ,σ2,c（c）)第页(年)∝第页(年∣μ,σ2,c（c）)第页(μ)第页(σ2)第页(c（c）).
(2)

在左边，我们有聚类参数和聚类分配的联合后验分布，我们要对其进行估计。在右边我们有第页(年∣μ,σ2,c（c）)，中的数据似然函数方程式（1）乘以先验值。粗体字体用于表示向量（例如，簇平均值μk个作为μ)。这一推导表明了贝叶斯统计中的一个重要概念——“后验概率是先验概率”(Lee，2012年，第2.1.2节）。从联合后验分布中，我们可以导出每个参数的分布。联合概率第页(μ,σ2,c（c）)变成第页(μ)第页(σ2)第页(c（c）)因为这三个参数假定是相互独立的。在这里第页(μ)第页(σ2)第页(c（c）)是下面描述的三个先前分布的乘积。

2.3.1. 优先级

我们需要指定之前μ,σ2和π。我们对所有簇平均值使用正常先验公共值，第页(μk个)~N个(μ0,σ02)，用协方差指定σ02例如，如果我们设置，第页(μk个)~N个(μ0=0,σ02=1)，然后我们显式地表示聚类平均值的标准正态先验。如果分析员对聚类平均值的先验知识有限，那么较大的误差幅度可能更好地反映不确定性，例如。，N个(μ0=0,σ02=16).在这里μ0和σ02是超参数参数的μk个.

逆Wishart先验分布适用于多元正态数据的协方差。最简单的多元正态数据遵循二元正态，如图1。在中使用符号Gelman等人（2003年，表A.1），Inv-Wishart(S公司−1,v（v）)previor可以设置为v（v）=2个自由度和一个比例矩阵美国。例如，在对二元法线进行聚类时，我们可以使用逆比例矩阵S公司−1=(100010)，表示方差（对角线条目）和零协变量（非对角线项）的高度不确定性，将由数据更新。

在建模单变量正态分布时，逆Wishart分布简化为逆Gamma分布形状参数γ和率参数β.

第页(σk个2∣γ,β)∼G公司−1(γ,β)∝(1∕σk个2)γ−1经验(−β∕σk个2).

中的派生附录A.3将基于单变量情况，以保持简单。注意超参数(μ0,τ0)和(γ,β)可以有自己的前科。超参数的先验被称为超验函数。为了简单起见，我们没有使用超验函数。

混合比例，πk个，在参数为α/K（K）:

第页(π1,…,πk个∣α)∼迪里克莱(α∕K（K）,…,α∕K（K）)=Γ(α)Γ(α∕K（K）)K（K）∏k个=1K（K）πk个α∕K（K）−1,
(3)

哪里α/K（K）为简单起见，=1表示平坦的Dirichlet分布。Γ（·）是伽玛函数，用希腊字母Γ表示。伽马函数由积分定义（例如。，霍格和克雷格，1995年第3.3节）。当gamma函数的输入是正整数时，输出就是我们熟悉的阶乘函数，偏移量为1Γ(n个)=============================================================================================================================================================(n个– 1)!, 它在本教程中使用（例如，Γ（7）=6×5×…×2×1=720）。伽马函数不应与上述定义的伽马分布混淆。读者可能会注意到，分布的形状只受以下因素的影响α因为K（K）是有限混合模型中的一个固定和已知常数。因此，K（K）本身不是参数。但是，α/K（K）后来起到了关键作用，当我们在容纳无限混合模型之前扩展Dirichlet，在该模型中簇数不再固定，它会发生变化，并可能接近无穷大。

不熟悉Dirichlet发行版的读者可能会发现图2有帮助。Dirichlet分布是范畴和/或多项式分布的共轭先验分布。从具有三个簇的Dirichlet分布中，我们可以得出比例样本[π1,π2,π三]=[0.5，0.3，0.2]，这将是方程式（1）.图2（a） ，α=[1，1，1]，表示权重均匀分布的平坦先验（点表示先验的潜在实现）。子图（b），带α=[0.2，0.2，0.2]，示出了在边缘附近具有权重的先验（例如[π1,π2,π三] = [0.50, 0.49, 0.01]). 子图（c），带α=[2,7,17]，表示先前认为第三个簇具有最高权重（例如[π1,π2,π三]=[0.2，0.1，0.6]）。中的值α可以概念化为先前的样本大小（例如。，α=[2,7,17]，之前的26个样本分布在3个亚组中）。的元素α不需要采用相同的值。有关其他视觉解释，请参见Li等人（2018）.

保存图片、插图等的外部文件。对象名为nihms-1528335-f0002.jpg

具有3个成分的假设混合模型的Dirichlet先验示例。Dirichlet分布是分类变量（如混合比例）的共轭先验分布[π1,π2,π三]，绘制在三维单纯形上。每个符号代表一种可能的实现[π1,π2,π三]产生于潜在的Dirichlet previor。的长度αparameter定义簇的数量，特定的参数值影响混合比例的分布。子图（a）表明了先前的观点[π1,π2,π三]可以呈现在三维单纯形上均匀分布的任何三元组。子图（b）显示了一种先验信念，即混合比例主要位于三维单纯形的边缘，其中一个混合比例接近于零（例如[π1,π2,π三] = [0.50, 0.49, 0.01]). 子图（c）显示了π三什么时候α= [2, 7, 17].

2.4. 条件后验分布

从关节后部分布方程式（2）我们可以解析地计算出μk个,σk个2、和πk个吉布斯采样。你所需要做的就是重组方程式（2）这样，您可以将结果识别为已知分布，从中可以提取样本。然而，如果我们使用方差的倒数而不是方差，这些代数操作会更简单。单变量正态分布中方差的倒数称为精度类似地，多元正态分布的精度是逆协方差矩阵。参见Gelman等人（2003年），第2.6节），以讨论贝叶斯统计中精度的作用。这些代数操作以簇精度显示在附录A.

集群意味着。

接下来我们介绍聚类协方差的精度，τk个，作为逆协方差矩阵k个聚类平均值为

第页(μk个∣年,c（c）=k个)∝第页(年,c（c）=k个∣μk个)第页(μk个)至N个(年‒k个n个k个τk个+μ0τ0n个k个τk个+τ0,n个k个τk个+τ0).
(4)

详细的逐步推导见附录A.2。这个等式起初可能看起来令人困惑和恐惧。但我们可以通过正常的魔术中已知的事实来接近它。这里我们试图确定k个具有多个正常分布数据点的第个簇，其中第页(年,c（c）=k个∣μk个)第页(μk个)基于贝叶斯规则。我们有μk个~N个(μ0,τ0)，使用τ0=1∕σ02表示关于先验均值的精度。我们还有一些观察到的数据点年，平均值为年‒k个=1n个k个∑k个[我=1]n个k个年我，其中k个[我]表示求和我嵌套在k个th簇并除以n个k个，该群集中的人数。我们可以将此信息表示为年~N个(年‒k个,τk个)，其中τk个表示集群平均值周围数据的剩余不确定性年‒k个在这种情况下（当我们有几个正常的观测值和一个正常的先验值的正常共轭时），已知以下是真的。样本平均值的后验分布μk个精度是先验精度加上样本精度乘以样本中的观测值数量，精确到n个k个τk个+τ0同上。请注意，如果我们使用精度而不是协方差来处理，操作会更简单。

这一事实可以在入门教科书中找到（例如。，Lee，2012年第2.3节）以及更高级的教科书（例如。，Gelman等人，2003年第2.6节）。我们可以考虑n个k个τk个+τ0作为先验精度中关于聚类均值的组合不确定性加上已知时关于聚类均值周围数据的剩余不确定性，这很直观。在引用的教科书中，你还可以找到后验分布平均值的推导。它是由先验精度加权的先验平均值加上由观测次数和样本精度加权的样本平均值，该总和除以权重之和，精确到年-k个n个k个τk个+μ0τ0n个k个τk个+τ0.

群集精度。

后验精度为：

第页(τk个∣年)∝τn个k个∕2经验(−n个k个∑[我]c（c）=k个(年我−年‒k个)2)×τα−1经验(−βτ)=τα−1+n个k个∕2经验(−τ(β+12∑[我]c（c）=k个K（K）(年我−年-k个)2)),
(5)

其中，先前的伽马分布G公司(γ,β)与数据相结合，生成具有更新形状参数的后验伽马分布α– 1 +n个k个/2和速率参数β+12∑[我]c（c）=k个K（K）(年我−年‒k个)2。有关详细信息，请参阅附录A.3.详细推导方程式（4）和(5)可以帮助初学者更好地理解为什么共轭模型采用特定形式。

潜在群集成员身份。

指标的联合分布c（c）1,c（c）2, … ,c（c）k个是混合比例的多项式实现πk个提高到n个k个，每个K（K）集群。

第页(c（c）1,…,c（c）k个∣π1,…,πk个)=∏k个=1K（K）πk个n个k个.
(6)

回想一下πk个遵循Dirichlet分布，密度函数由参数控制α(方程式（3）)，而且c（c）k个指标是实现πk个因此，后验Dirichlet分布具有更新的参数n个1+α/K（K）–1，n个k个+α/K（K）– 1.

什么时候？K（K）有限且已知，采样模型相对简单。我们可以转向方程式（3）–(6)吉布斯采样算法：

πk个(t吨)∼迪里克莱(n个1+α∕K（K）−1,…,n个k个+α∕K（K）−1),μk个(t吨)∣c（c）(t吨−1),年,τk个(t吨−1)∼N个(年‒k个n个k个τk个(t吨−1)+μ0τ0n个k个τk个(t吨−1)+τ0,n个k个τk个(t吨−1)+τ0),τk个(t吨)∣c（c）(t吨−1),年,μk个(t吨−1)∼伽马射线(α−1+n个k个∕2,β+12∑[我]c（c）=k个(年我−年‒k个)2),c（c）我(t吨)∼多项式(2,π1(t吨)N个(年我∣μ1(t吨),τ1(t吨)),…,πk个(t吨)N个(年我∣μk个(t吨),τk个(t吨))).

我们使用上标(t吨)和(t吨–1）分别对最近采样的参数和之前的参数进行索引。在最后一个等式中，潜在的簇成员是潜在混合比例的多项式实现，由相应簇的数据可能性加权（例如。，c（c）我，=（0，1，0）产生于π= (0.4, 0.5, 0.1)). 我们省略了一些细节，如最后一行中的归一化常数，因为解释如何组装有限混合模型并不是我们的直接目标。相反，本练习强调了将这些结果推广到无限混合物中的挑战。重点是显示各种概率分布，它们看起来像什么，以及如何通过迭代采样导出它们。对吉布斯采样细节感兴趣的读者可以阅读我们关于有限组件模型的教程(Li等人，2018年).

如果…怎么办K（K）不是有限的和已知的吗？上面的算法不再有效，因为它需要一个固定的K（K）。更具挑战性的是，如果K（K）接近无穷大？不可能从具有无穷向量的Dirichlet分布中采样π的。这似乎没有希望了。然而，弗格森（1973）通过计算Dirichlet过程先验来解决这个问题，这是一种不依赖于无穷大的先验分布π的。下面将对此进行说明。

三。 无限混合模型

在这一节中，我们将讨论无穷多个簇的问题。我们的目的是使数学解释明确而具体。方程式的推导出来后，它们将与一个称为中餐厅流程（CRP，奥尔德斯，1985年;Blackwell&MacQueen出版社，1973年)更直观地表示DP。研究该主题的读者可能会遇到CRP，但不一定会遇到解释CRP工作原理的方程式。我们的目标是将它们联系起来，以便对DP有更全面、更全面的了解。CRP并不是DP的唯一代表。还存在其他表示（例如，通过Sethuraman（1994）和Pólya urnBlackwell&MacQueen（1973）)，相当于CRP。

3.1. 无穷多簇问题的处理

首先要解决的问题是混合比例的无限向量。回想一下π是的函数α(方程式（3）)还有那个C类我是向量的函数π(方程式（6）)。如果我们把这两个功能结合起来π，然后我们可以C类我直接产生于α。以下是中显示的内容拉斯穆森（2000）:

第页(c（c）1,…,c（c）k个∣α)=∫第页(c（c）1,…,c（c）k个∣π1,…,πk个)第页(π1,…,πk个∣α)d日π1⋯d日πk个=Γ(α)Γ(α∕K（K）)k个∫∏k个=1K（K）πk个n个k个+α∕K（K）−1d日πk个=Γ(α)Γ(n个+α)∏k个=1K（K）Γ(n个k个+α∕K（K）)Γ(α∕K（K）).
(7)

小写k个用作特定集群的索引。大写字母K（K）用于表示集群的总数。在下一节中，我们将看到DP允许K（K）表示不同数量的簇，无论是有限的还是无限的。A固定K（K）DP中不再需要。在整个教程中都使用了这种表示法。

请注意π从最后一个等式中消失。这里的一个重要见解是，可以通过将有问题的参数集成到模型之外来解决它。从第1行到第2行的转换相对简单。对称Dirichlet先验方程式（3）乘以多项式混合比例方程式（6）由于这两个方程都涉及πk个，乘法意味着将指数相加，以及Γ(α)Γ(α∕K（K）)k个是一个常数，所以它被移到积分符号的前面。然而，第2行和第3行之间的转换并不明显。技术专家可能会认为省略的步骤微不足道，但初学者需要这些步骤才能更全面地理解之前的DP。拉斯穆森（2000）解释说它只是一个“标准狄利克雷积分”，没有提供进一步的解释。尼尔（2000）完全跳过了这个，直接转到方程式（8）如下所示。Navarro等人（2006年）提供了更多细节，但不容易访问。因此，我们在附录A.4(方程式（A.6）–（A.7）具体来说）。鼓励读者循序渐进地阅读，以了解方程式（7）.

方程式（7）使我们离完全解决从无限数量的混合物比例中采样的问题又近了一步。然而，它仍然需要集群的总数K（K）以获取特定值。如果K（K）是允许接近无穷大，然后我们需要一种方法来解决它。接下来的几步可能是最具挑战性的。我们将非常明确地处理阐述和注释中的细节，比介绍性文本中通常提供的步骤更具体地解释步骤。

3.2、。当K接近无穷大时

方程式（7）可以进一步利用以生成一个完整的解决方案K（K）我们已经计算出单个指标的条件优先权C类我考虑到除c（c）我.

第页(c（c）我=k个∣c（c）−我,α)=n个−我,k个+α∕K（K）n个−1+α,
(8)

其中下标−我在里面c（c）−i表示除我th和n个−i、 k个表示嵌套在k个第个集群除外年我. Then个分母中是观察的总数。详细推导很长，因此在附录A.4(方程式（A.7）–（A.8）)。鼓励读者仔细遵循这些步骤，以了解原因方程式（8）是真的。

3.3、。可交换性

的直接结果方程式（8）是互换性属性-簇分配与单个数据点的顺序不变的概率。每个数据点都可以被视为此条件分布中的最后一个数据点。将哪个数据条目分配给哪个集群并不重要。等式右侧的任何内容都不受观测顺序的影响。概率主要受嵌套在k个第个集群。集群成员身份是按任意顺序一次确定一个数据条目（以现有集群安排为条件）。这还意味着，只要两个集群大小相同，那么它们的概率就相同n个−i、 k个交换性在DP中至关重要，因为它解决了从无限多类别的Dirichlet分布中采样的问题。

如果我们允许K（K）→ ∞ 在方程式中(8)，然后是α/K分子接近零的极限，因此消失。

集群，其中(n个−我,k个>0:第页(c（c）我∣c（c）−我,α)=n个−我,k个n个−1+α.
(9)

假设在一个假设的例子中，前10个观测值已经被分为两个集群，分别有3个和7个观测值。第11个观测被分配给任何一个集群，其概率与已经占据各个集群的观测数量成比例，即3:7的比例。如果我们总结一下n个−i、 k个在已经被占用的集群中，我们得到1010+α，第11个观测值被分配给任何已占用集群的概率。注意，如果我们添加α10+α到1010+α，然后我们得到10+α10+α，正好是1。是什么α10+α? 这是第11次观测的概率不属于任何已被占用的集群。这是第11个观测值被指定为新的，以前为空群集。一般来说，如果我们总结一下n个−i、 k个在已经被占据的星团的分子中，我们得到n个−1，不包括最后一个观察的总数。与上述假设示例类似，1−n个−1n个−1+α=αn个−1+α给出了最后一个观测值被分配给新簇的概率。

我们详细讨论这一重要点，以确保它绝对清楚。对于任何个人观察年我，它可以放在现有集群中，也可以是新的初始空集群的第一个成员。只有这两种情况，所以它们的概率总和必须为1。场景一涉及将第i个观测值放入现有集群之一（其中n个−i、 k个> 0). 在已占用集群中吸引新成员的概率与n个−i、 k个，每个集群中已有的观察数。情景一的总概率为∑k个=1K（K）n个−我,jn个−1+α注意分子，n个−我,1+n个−我,2+ · · · +n个−i、 K（K）总计n个−1，除尚未分类的单个观测外，所有现有集群中聚集的观测总数。Thu场景一的总概率为n个−1n个−1+α。如果我没有将观测值放置在任何已占用的簇中，则必须创建一个新簇。这是场景二。但创建新集群的概率是多少？它正好是1减去场景一的概率。如前所述，场景一和场景二的概率之和正好为1。因此，第二种情况的概率为

1−∑k个n个−我,k个n个−1+α=n个1+αn个−1+α−n个−1n个−1+α=(n个−1+α)−(n个−1)n个−1+α=αn个−1+α.

因此，综合所有因素，集群分配的概率为：

集群，其中n个−我,k个>0:第页(c（c）我∣c（c）−我,α)=n个−我,k个n个−1+α,所有其他集群合并:第页(c（c）我≠c（c）k个∀j≠我∣c（c）−我,α)=αn个−1+α.
(10)

DP的一个重要特性体现在这些方程中。创建新集群的概率与浓度参数α(Neal，2000年)，如第二行所示。更大的值α倾向于鼓励创建新的集群，尤其是在集群规模较小的情况下。例如，如果α=3和n个−i、 k个=3，则两个概率相等。如果混合权重给定Dirichlet的对称Dirichle先验(α/K（K）,…,α/K)然后，Dirichlet后验函数产生可交换性，这反过来又允许生成新簇的迭代过程，每次一个观察值。中给出了等效的证明Gershman&Blei（2012年，他们的方程式（8）).

3.4. 中餐厅流程

DP通常用一个叫做中餐厅流程（CRP）的烹饪隐喻来解释Gershman&Blei，2012年)，使用与中类似的图表图3.

保存图片、插图等的外部文件。对象名为nihms-1528335-f0003.jpg

中餐厅流程图。在CRP比喻中，想象一家中国餐馆有无数张桌子。顾客（个人数据条目，显示为菱形）一个接一个地进入餐厅，并按照进入餐厅的顺序坐在桌子旁（离散的簇，显示为圆圈）。有与集群关联的参数，表示为φk个= [μk个,τk个]用于聚类平均值和精度。第一个进入餐厅的顾客总是坐在第一张桌子上。第二位顾客以1/（1）的概率进入并坐在第一张桌子上+α)，和第二个概率表α/(1+α)，其中α是一个正实数（最上面一行，其中i=2）。当第三位顾客进入餐厅时，他或她坐在每一张有人的桌子旁，概率与已经坐在那里的前几位顾客的数量成正比，而坐在下一张没有人的桌子上，概率与α.

想象一下，一家中国餐厅有无限的座位。客户（数据条目，显示为菱形）坐在桌旁（集群，显示为圆圈）。表格由参数标识ϕk个= [μk个,τk个]聚类的平均值和精度（为了便于比较，我们使用了它们的符号）。第一位顾客走进餐厅，他或她坐在一张最初空着的桌子旁。当第二位顾客到达时，他或她可以和第一位顾客坐在已经有人的桌子旁，也可以坐在没有人的桌子上。假设她坐在一张空桌旁。现在我们有两张桌子被占用了。当第三位顾客到达时，她坐在表1和​和22概率与已经坐在那里的顾客数量成正比，而在一张空位上概率与浓度参数α。相同的迭代过程一次循环一个客户，直到为所有客户分配一个表。CRP从（原则上）无限多的表开始。然而，在有限数据集中，模型复杂性由DP决定。并不是所有无限多的表都是在给定数据的情况下实现的。

表1：

DP中的关键方程与如何在R中实现它们之间的比较。最后一列涵盖了执行计算的主要R函数。通过dmvnormO函数对新观测到的数据项进行评估，以获得其对数密度，并为预测分布提供适当的输入参数。喜欢不同计算机编程语言的读者可以使用这些语言来指导其实现。

DP方程行R代码
n个−我,k个n个−1+αN个(年~我;年‒k个n个k个τk个+μ0τ0n个k个τk个+τ0,1n个k个τk个+τ0+σ年2)
方程式（12）		108 – 110平均_p<-sig_p%*%（tau_y%*%总和数据+
tauO%*%t（muO））
logp[ci]<-log（nk[ci]）+dmvnorm（数据[n，]，
平均值=平均值_p，σ=sig_p+σ_y，log=真）
αn个−1+αN个(年~我;μ0,σ02+σ年2)
方程式（13）		117 – 118logp[Nclust+1]>-log（alpha）+dmvnorm（数据[n，]，
平均值mu0，sigma=sigma0+sigma_y，log=TRUE）

表2：

分类结果在各种模型假设下恢复一个模拟的4个集群，每个集群60个。在行下，我们有越来越多的α，这将鼓励创建新的集群。从左到右的纵列中，我们的测量误差越来越大σ年2，这将阻止创建新集群。

σ年2=0.50σ年2=1σ年2=1.50σ年2=3
α= 0.01c（c）我12三4512三412三1
n个k个6433636020796735591453065240
α= 1.00c（c）我12三4512三412三41
n个k个65603022636036687687306657240
α= 3.00c（c）我12三4512三412三41
n个k个66602427636860347836637437240
α=5.00c（c）我124512三41三561
n个k个7561337184616431337911315240

为了使CRP起作用，顾客进入餐厅的顺序必须无关紧要。这是互换性上述财产。CRP比喻很吸引人，因为它很直观，并且很好地映射到方程和每次一个观测值的迭代采样上。然而，对尼尔（2000）和拉斯穆森（2000）显示CRP未完成。必须进一步处理集群分配概率。这里我们离题讨论一个重要的贝叶斯概念，称为后验预测分布。

3.5. 先验和后验预测分布

在贝叶斯统计中，有两种类型的预测分布。我们首先转向后验预测分布或简称PPD。PPD是新观察到的基于您已经看到的数据的值(Christensen等人，2011年)。第i个人的新观察值通常表示为年~我。直觉上，如果一个新观测点距离某个星团很远，那么它属于该星团的概率应该很低，反之亦然。请注意公式（10）取决于所有群集指示符不包括的年~我，它没有考虑新的观测结果与任何已被占据的星团的吻合程度。因此，我们需要将年~我关于每个k个第个集群。

已占用集群的PPD。

PPD通常在以下一般设置中引入（例如。，Gelman等人，2003年第3.6节；Lee，2012年第2.3节；和Novick&Jackson，1974年，第6.4节）。假设一个变量年正态分布，平均值未知μ测量误差的形式为固定且已知方差σ年2也就是说，

年∼N个(μ,σ年2).

此外，如果未知平均值μ表达了先前的信念

μ∼N个(μ0,σ02),

然后是一个新观测值的后验预测分布年~根据你已经看到的数据

年~∼N个(μ0,σ年2+σ02).

新观测数据点的密度年~通过平均值的正态分布进行评估μ0和方差σ年2+σ02这里的直觉是，组合协方差是两个变量来源的总和，即测量误差引起的方差，σ年2以及未知平均值的方差，σ02。在继续之前，理解这一点很重要。如果仍然有点混淆，那么Lee（2012年（第2.3节）可能有助于解释后验分布和后验预测分布之间的差异。

我们现在准备将PPD的一般情况应用于k个第个集群。我们已经知道了方程式（4）第k个聚类平均值的条件后验分布是正态的，平均值为μ第页=年-我n个k个τk个+μ0τ0n个k个τk个+τ0精度为τ第页=n个k个τk个+τ0因此，未知平均值的方差为σ第页2=1∕(n个k个τk个+τ0)换句话说，聚类平均值遵循分布N个(μ第页,σ第页2)此外，如果我们假设观测值是以固定的已知误差测量的σ年2，然后我们可以从一般情况中推断出新观察到的PPDk个第个簇：

年~∼N个(μ第页,σ第页2+σ年2),=N个(年‒我n个k个τk个+μ0τ0n个k个τk个+τ0,1(n个k个τk个+τ0)+σ年2).
(11)

与一般情况类似，组合方差是测量误差的总和σ年2以及未知簇平均值周围的方差，σ第页2这里我们假设测量误差σ年2在每个潜在集群中都是已知且相同的，这是一种粗略但实用的近似，使示例更容易理解。在Bernardo&Smith（2000）中可以找到各种共轭情况下的后验预测分布目录。电子源可以方便地位于https://enwikipedia.org/wiki/Conjugate_prior网站.

先验预测分布。

与PPD相关的一个概念是先验预测分布（Gelman等人，2004年，第1章；Lee，2012年; 第2章；Novick&Jackson，1974年第142页）。假设一个变量年遵循具有未知平均值的正态分布μ，以及一个固定的已知测量误差σ年2此外，μ遵循平均值的先验分布μ0和方差σ02.前面的参数μ0,σ02以及数据的不确定性，σ年2，假设已知。Novic&Jackson（1974）表明年~遵循平均值的正态分布μ0和组合方差σ02+σ年2本质上，未来观察的不确定性年~是关于聚类平均值的综合不确定性(σ02)集群周围数据的剩余不确定性意味着当已知时(σ年2).

在实现新的聚类之前，没有定义其平均值和精度。因此，新观察到的年~我由先验引起的结果必须由先验预测分布来评估。我们现在准备将两个基本结果编织在一起：1）CRP用于指导集群分配；和2）先验和后验预测分布，以评估由聚类产生的单个观测的概率。接下来，我们将两者结合起来完成DP混合模型。

c的条件后验分布我.

当结合后验和先验预测分布时，CRP分配概率方程式（10）指标的条件后验分布如下：

  • 集群，其中n个−i、 k个> 0 :
第页(c（c）我∣c（c）−我,μk个,τk个,α)∝第页(c（c）我∣c（c）−我,α)第页(年~我∣μk个,τk个,c（c）−我)∝n个−我,k个n个−1+αN个(年~我;年‒k个n个k个τk个+μ0τ0n个k个τk个+τ0,1n个k个τk个+τ0+σ年2).
(12)
  • 所有其他集群的组合：
第页(c（c）我≠c（c）k个∀j≠我∣c（c）−我,μ0,τ0,γ,β,α)∝第页(c（c）我≠c（c）k个∀j≠我∣c（c）−我,α)∫第页(年~我∣μk个,τk个)第页(μk个,τk个∣μ0,τ0,γ,β)d日μk个d日τk个∝αn个−1+αN个(年~我;μ0,σ02+σ年2).
(13)

方程式（12）包含已占用集群的PPD。新观测到的年~我（即要分类的最后一个数据条目）通过相应PPD的手段和精度进行评估。方程式（13）根据CRP将先前的预测分布合并到创建新集群中。

方程式（13）包含一个看似可怕的积分。然而，它可以分为更简单的组件。第一，第页(μk个,τk个∣μ0,τ0,γ,β)表示可能由先验值产生的聚类平均值和精度矩阵的值。接下来，每次实现集群参数时μk个,τk个，计算相应的第页(年~我∣μk个,τk个)，新观测到的年~。这一点在所有可能出现的由先验值产生的簇上重复。最后，积分将所有可能实现的簇的结果相加。这也是如何正式定义先验预测分布的(Gelman等人，2003年第1.3节）。在其他教程中，有时会对集群的所有随机实现进行总结(Neal，2000年;拉斯穆森，2000年)作为∫N个(年~我;ϕ)G公司0(ϕ)d日ϕ，其中ϕ= (μk个,τk个)，是相同概念的缩写形式，其中G公司0(ϕ)表示先验分布的函数，所有簇平均值和协方差都可能由此产生。

一个更具体的例子可能有助于解释这些方程式。在图3，当第9位顾客进入餐厅时，他或她有可能坐在已经客满的位置表1–4与矢量成比例n个−i、 j个分别=（2，4，1，1）。此外，我们需要考虑第9位客户与已入住餐桌的匹配程度。例如，如果发现该客户距离表2（假设(年1,年2)=（1.3，−1.1），最接近μ2=（1.5，−1.5）英寸图1)然后，将他/她分配到表2是很直观的。这个新客户与每个已经占用的表格的接近程度也是一个向量，由PPD根据之前的8个观察结果进行评估方程式（12）。两个向量的乘积是第9位客户坐在已被占用的桌子上的概率。然而，事实证明，第9位顾客坐在图3。该概率取决于方程式（13），按产品α/(n个−1+α)以及先验预测分布。重要的是要记住，CRP是一个随机过程，因此特定的座位安排是DP的概率实现。它不是固定的。

3.6. DP、CRP和DPMM之间的关系

方程式（12）–(13)总结DPMM的必要推导。我们准备在一个工作的计算机程序中实现这些结果。然而，这一点很好地解释了DP、CRP和DPMM之间的关系。DP在引言中描述为分布上的分布因为它将具有固定类别数的普通Dirichlet分布扩展为具有可变类别数的无穷多分布k个类别（回忆以下步骤方程式（7）和(10))。CRP可以理解为分区上的分布（如何对样本进行分区），而不包括预测分布。DPMM是一种混合模型，它在考虑理论上无限数量的混合组分之前使用DP，其中包括预测分布(方程式（12）和(13))。接下来，它们将在一个工作的计算机程序中实现。

4 DPMM算法

算法1一次对一个数据条目的集群分配概率进行重复采样。它相当于中的算法3尼尔（2000）算法从初始参数和最大迭代次数开始。然后进入类似CRP的流程，其中年我，的我由于将在当前迭代中对新值进行采样，因此首先从其簇中删除第个观测值。接下来，方程式（12）–(13)用于计算年我属于其中一个K（K）+1个群集，带有K（K）+1表示新创建的集群。注意，采样算法不同于通常的吉布斯采样，因为采样是一次输入一个数据。此外，一个人可以在每次迭代中坐在不同的CRP表上。算法在达到最大迭代次数时结束。将返回迭代期间的集群成员身份。然后我们可以计算我第个观察结果属于集群1，k个在迭代过程中。我们可以选择使用概率最高的集群来表示我第th次观察。请注意，第8-9行基于Equaions公司（12）和第10-11行基于方程式（13）.

这引发了一个问题。例如，如果在迭代100时将一个观测值分配给集群2，然后在迭代800时再次分配给集群2中，那么这真的是同一个集群吗？通常是这样的。下面是使用CRP隐喻的直观解释。因为CRP一次只应用于一个客户，所以坐在同一张桌子上的其他客户可以提供相对于桌子成员分配的稳定性。随着更多客户参与CRP，拥有更多客户的桌子往往会吸引更多客户。此外，如果客户具有属于特定组的高后验预测概率，并且该组中已经有许多客户，那么根据所有其他客户的组成员分配情况，将一个客户置于完全不同的组中的可能性相对较小。靠近两个重叠簇边界的观测可能会在迭代期间在它们之间切换。但集群名称仍应保持稳定。当然，也有例外。CRP分配在开始时可能不稳定，因为只有前几个客户坐着。被分配到表1、表2和表3的概率可能大致相当。表赋值的顺序可以采用不同的排列，并且总体概率分布不受影响。这些最初的观察可能会遇到类似于众所周知的问题标签交换贝叶斯混合模型中的问题(斯蒂芬斯，2000;Jasra等人，2005年;Richardson&Green，1997年)。然而，当更多的客户通过CRP时，这在以后变得不太可能。PPD提供了额外的稳定性。显然，样本量不足加剧了问题。那么问题实际上是由于数据不足，而不是算法中的固有限制造成的。

保存图片、插图等的外部文件。对象名为nihms-1528335-f0004.jpg

5 DPMM计算的实现

5.1. R程序逐行

附录B显示了如何在R中实现算法1，根据下面给出的示例7进行了修改布罗德里克（2017）。更改原始程序中的变量名称以匹配我们的符号。这个练习的目的是使方程具体而明确。它还显示了计算机编程如何补充方法开发。因此，我们将重点关注最不重要的代码行。在整个示例中添加的注释应该使编程代码更容易理解。

前28行包含许可证信息。第29行定义了一个名为crp_gibbs（）的函数及其接受的输入参数。它接受原始数据、浓度参数α未知聚类均值的先验均值和方差(μ0,σ02)，测量误差方差(σ年2、假定固定和已知）、初始集群分配和最大迭代次数。下一个具有实质重要性的行是第63行，CRP从这里开始。接下来，在第64行，我们开始一次一个客户的CRP，从for（nin1:n）循环开始。接下来我们删除年我因为我们要画一个新的样本c（c）我对于年我。如果这导致表变空，那么我们将用最后一个集群（第75–79行）替换新的空集群，以便Nclust始终指向已占用的表，而下一个群集（Nclus+1）始终指向要创建的下一个新表（如果需要）。这些在DPMM算法1中进行了描述。下一行的实质重要性是第89行，我们在这里计算τ第页= (n个k个τk个+τ0)根据方程式（11）。然后我们反转τ第页得到σ第页2稍后使用。

第107行计算后验均值mean_p。请注意，如果Nclust包含3个已占用的簇，则mean-p将是长度为3的向量。将后验均值输入到第108–109行中，以计算CRP权重（log（n_k[c_i]）），并使用dmvnorm（）函数将其乘以对数后验预测密度，并根据方程式（12）。对数使计算更容易，因为它们将乘法转化为加法。第116–117行根据方程式（13）最后，第124行使用我们刚才计算的概率绘制了这个新客户应该属于哪个集群的样本。此客户可以属于任何1:（NClust+1）集群。其余行收集并返回结果。

然后，我们可以使用接下来的几行来拟合模型，并找到类成员估计的最大后验概率（后验模式）。

>多元正态函数库（“MASS”）#mvrnorm（）>set.seed（11）#再现性的随机种子集>n<-60#分别从4个独立分布中抽取60个样本>m1<-c（1.5，1.5）#右上x，y表示>S1<-矩阵（c（0.3，0.05，0.05，0.3），ncol=2）#c1的方差>#根据平均值mu和方差Sigma从每个簇中抽取n个样本>clus1<-mvrnorm（n=n，mu=m1，Sigma=S1）>m2<-c（1.5，−1.5）#右下角>S2<-矩阵（c（0.5，−0.08，−0.05，0.2），ncol=2）>clus2<-mvrnorm（n=n，mu=m2，Sigma=S2）>m3<-c（−1.5，1.5）#左上>S3<-矩阵（c（0.1，0.03，0.03，0.1），ncol=2）>clus3<-mvrnorm（n=n，mu=m3，Sigma=S3）>m4<-c（−1.5，−1.5）#左下>S4<-矩阵（c（0.8、0.50、0.50和0.8），ncol=2）>clus4<-mvrnorm（n=n，mu=m4，Sigma=S4）>datc<-rbind（clus1，clus2，clus3，clus4）总共240个观察结果>#运行附录B中的CRP Gibbs函数。>α<-0.01>mu0<-矩阵（rep（0，2），ncol=2，byrow=TRUE）>σ0<-诊断（2）*3^2>σ_y<-diag（2）*1>c_init<-rep（1，nrow（datc））>结果<-crp_gibbs（data=datc，alpha=alpha，mu0=mu0，σ0=σ0，σy=σy，cinit=rep（1，nrow（datc））> #####>#Gibbs采样迭代保存在“results”中。每行包含>#MCMC上一个人的潜在阶级成员分布>#迭代次数。如果某人被视为c=2的成员，例如>#最常分配给c2。所以我们将频率制成表格>#计数，频率最高的集群是>#人的潜在类。> #####>tab<-应用（结果，1，乐趣=函数（x）{tab<-表格（x）ans<-名称（选项卡[which.max（选项卡）]）返回（ans）}）>表（选项卡）选项卡1 2 3 476 39 65 60

表1提供了以下各项之间的并排比较方程式（12）和(13)和相应的R命令。最后一个数据条目通过适当正态密度的预测分布进行评估，并通过CRP分配概率进行加权。将适当预测分布的参数输入dmvnorm（）函数中，以获得预测对数密度，然后用CRP赋值概率对其进行加权，从而得出高达K（K）+1个集群。使用矩阵乘法（%*%），使tau_y（一个2乘2矩阵）符合sum_data的2乘1矩阵，从而为后聚类平均值生成正确的x-y坐标。然后将其左乘以sig_p，sig_p是τ第页=n个k个τk个+τ0，得出与除以相同的结果τ第页。计算不必在R中实现。Python程序员可以使用numpy.random.normal（）函数，并通过设置输入参数loc=mean_p，scale=sig_p+sigma_y。

5.2. 模拟练习

初学者通常通过玩模型的设置来学习更多关于模型行为的知识。两个输入参数，α和σ年2，在DPMM中发挥重要作用。我们将研究它们是如何影响聚类的。回想一下，假设数据是使用来自四个不同多元正态分布（共240个数据点）的60个模拟观测值构建的。虽然此程序在没有得到指示的情况下成功恢复了四个不同的集群，但每个集群的不同成员大小表明存在错误分类。误分类部分是由于算法的概率性质。其他因素也起作用。例如α影响创建新集群的概率。较大的值往往会鼓励创建新集群。每个集群中的测量误差，σ年2对于所有集群，当前设置为相同的1.0，也会影响集群（在第3.5节中解释）。如果数据点距离聚类中心的距离大于先验方差，并且靠近分布的重叠边界，则数据点更有可能被错误分类。虽然当前的程序很简单，但它将能够从不同且部分较远的分布中对样本进行分类。下面我们将说明读者如何尝试不同的设置来查看它们如何影响集群。这里要说明的最重要的思想是，实现并不困难。R语法和数学之间的相似性也有助于更全面地理解该理论。

表2显示了参数的输入值如何影响CRP结果。接下来的几行在以下所有组合上循环α=（0.01，1.00，3.00，5.00）和测量误差σ年2=（0.50，1.00，1.50，3.00），以检查它们如何影响聚类。

>设定种子（11）>α<-c（0.01、1.00、3.00、5.00）>σ_y<-c（0.50、1.00、1.50、3.00）>for（i in 1：长度（α））+for（j in 1：长度（sigma_y））+    {+结果<-crp_gibbs（数据=datc，α=α[i]，+mu0=矩阵（rep（0，2），ncol=2，byrow=TRUE），+sigma0=诊断（2）*3^2，sigmay=诊断（二）*sigmay[j]，+c_init=代表（1，nrow（datc））+tab<-应用（结果，1，FUN=函数（x）+{制表符<-表格（x）+ans<-名称（选项卡[which.max（选项卡）]）+返回（ans）}）+cat（“alpha=”，alpha[i]，“sigma_y=”，sigma_y[j]，“\n”）+打印（表格（选项卡））+    }

输入参数的选定值提供了算法在实际数据分析中的行为一瞥。的值α被选中，以便我们有一系列大小值。安α=1.0在创建新簇时是中性的，因此可能是一个好的默认值。更大的α价值观往往鼓励创建新的集群。

最大的σ年2=3之所以选择，是因为它近似于在x轴或y轴方向上240个观测值的整个样本的观测方差图1它反映了一种最坏的情况，在这种情况下，整个样本被假定在3.0的固有测量误差范围内，有效地将整个样本视为一个单独的簇。它还意味着，除非（粗略地说）它们之间的距离超过测量误差或3.0，否则无法区分簇。因此σ年2可以被认为是发现可识别簇的假定精度。一个σ年2=1近似为中扩散最广的簇的方差图1.最小测量误差σ年2=0.50接近中前两个簇的观测方差图1中等规模的集群。我们以中等规模集群的大小为指导，表达了预期的精确度。在这里，我们用一种粗略的、务实的视觉检查来指导设置。然而，分析员应努力使用可靠的来源来估计测量误差（例如Lee，2012年第2.2节）。

转向表2，我们首先查看最右边的列。什么时候？σ年2=3设置时，我们只得到一个集群，因为这是我们指示模型执行的操作，而不考虑的值α。当我们将测量误差减半到1.5时，我们开始在不同的α值。然而，错误分类是常见的。如果我们设置σ年2=1，然后我们得到了4个聚类的更一致提取，不同聚类之间的错误分类更少α值。对于σ年2=0.50我们倾向于得到五个簇。更仔细的检查（未显示）表明，集群4（分布最广）的一些成员被分组到第五个集群中。我们在中观察到跳过的类成员身份c（c）我当我们设定时=1,3,5,6α=5.00和σ年2=1.50可能是因为在迭代过程中会创建多达6个集群，但最终只有4个集群出现。总的来说，表2表明当前实现对σ年2、和σ年2=1（接近集群中观察到的差异，且分布明显最广）似乎提供了最理想的结果。分析师可以使用集群中的视觉检查或其他观察到的经验特征来指导设置。

下面几行代码说明了如何获得有关错误分类的更多信息。下面的变量c_true包含真正的类成员身份。假设我们设置σ年2=1和α=1.0时，我们运行crp_gibbs（）程序并获得变量c_modal，该变量包含DP迭代中最可能的集群成员身份。apply（）函数获取结果（一个240人乘以1000次迭代的矩阵），并将每个人在1000次迭代中最可能的集群成员分配列表，一次一行。

>设定种子（11）>结果<-crp_gibbs（数据=datc，α=1.0，+mu0=矩阵（rep（0，2），ncol=2，byrow=TRUE），+σ0=诊断（2）*3^2，σy=诊断（1）*1.0，+c_init=代表（1，nrow（datc））>c_modal<-apply（结果，1，FUN=函数（x）{tab<-表格（x）ans<-名称（选项卡[which.max（选项卡）]）返回（ans）} )>c_true<-rep（1:4，每个=60）>表格（c_true，c_modal）c模态c真1 2 3 41 60  0  0  02  0 60  0  03  0  0  0 604  0  7 36 17

这个交叉制表提醒我们：1）DP的随机性产生了c_modal分布，与之前的解决方案略有不同表2; 2） 聚类的顺序没有特别的重要性，它们只是离散的标签（真正的聚类1、2、3和4被不同地标记为1、2、4和3）；和3）错误分类发生在第四个簇中（回想一下，簇4是图1)，其中7个个体被错误归类为第2组（中的右下角图1)17个错误归类为第3组（左上角）。然而，集群4中的错误分类并不奇怪，因为它具有很高的可变性，并且靠近相邻的集群2和集群3。请注意中重叠的簇2和4椭球体图1。对于重叠区域中的数据点，它们的簇成员身份本质上是不确定的。本练习的目的是表明表2仅提供有关DPMM程序如何执行的有限信息。将错误分类分解为更多细节是有意义的。

5.3. 估计聚类均值和协方差

DPMM的用户可能对集群特征感兴趣（例如，均值和协方差）。集群意味着μk个使用计算起来很简单方程式（4）。簇协方差的参与程度稍高，部分原因是方程式（5）仅适用于单变量正态分布。幸运的是，我们可以找到计算二元正态协方差的公式（见Bernardo&Smith（2000）和https://en.wikipedia.org/wiki/Conjugate_prior网站).

下面的前5行对先验均值、先验协方差和先验精度进行了必要的准备。接下来，我们使用for（）循环一次计算一个簇的后验均值和协方差。请注意，我们在中重用了一些R代码附录B为了计算平均值pk个第个集群，μk个然后我们转向集群协方差。根据上面的维基百科页面，如果先前的协方差遵循逆Wishart分布v（v）=2个自由度和一个比例矩阵S公司0=σ0，则后验逆Wishart分布为n个k个+v（v）自由度和后标度矩阵S公司0+∑我=1n个k个(年我−μk个)(年我−μk个)T型. The∑我=1n个k个(年我−μk个)(年我−μk个)T型项是聚类平均值的平方偏差之和，通过首先取偏差（从中减去平均值年1和年2使用sweep（）函数），以便可以使用crossprod（）函数对它们进行平方和求和。此外，逆Wishart分布的平均值为S公司第页ν第页−第页−1，其中S公司第页,v（v）第页分别是比例矩阵和自由度，以及第页=2，因为它是二维（二元）法线。下面我们得到S公司第页=（σ0+交叉prod（y_ss））和v（v）第页=n_k[c_i]+2，剩下的应该很简单。

mu0<-矩阵（rep（0，2），ncol=2，byrow=TRUE）σ0<-矩阵（c（3，1，1，3），ncol=2）σy<-矩阵（c（1,0.2，0.2，1），ncol=2）tau0<-solve（sigma0）#逆先验协方差tauy<-solve（sigmay）#逆测量误差nk<-table（cmodal）#每个集群中的人数for（ci in 1：长度（nk））{tau_p<-tau0+n_k[c_i]*tau_y#R程序行89sigp<-求解（taup）y<-datc[其中（c_modal==名称（n_k）[c_i]），]sum_data<-colSums（y）#R程序行95平均值p<-sig_p%*%（tau_y%*%sum_data+tau0%*%t（mu0））#108y_ss<-扫描（y，MARGIN=2，STATS=mean_p，FUN=“-”）协方差<-（sigma0+交叉prod（yss））/（nk[ci]+2−2−1）cat（“\n cluster”，名称（n_k）[c_i]，“\n”）打印（圆形（平均p，3））打印（圆形（协方差，3））}

结果（可通过运行代码获得）很粗糙，但与中的真实参数值相当接近图1考虑到计算机程序的简短性和可访问性。精度可以通过更多的迭代次数来提高，或者在计算程序的更高版本中包括更新测量误差的步骤σ年当前已修复。可能包括其他步骤。我们可以将这些代码行添加到附录B并返回参数值。然而，这样做不仅增加了编程代码的复杂性，而且由于额外的计算和内存存储，可能会大大降低程序的速度。当然，这并不禁止读者修改代码和设置，以查看它们如何影响结果。事实上，我们鼓励对计算机程序进行修补，这无疑是学习任何计算算法的最佳方法。

6 讨论

本教程旨在涵盖四个主题：1）Dirichlet过程之前的严格推导；2） 将DP与中餐厅流程施工方法联系起来；3） 关于如何在二元混合模型中实现DP的假设示例；和4）简单的模拟研究，向读者展示DP的特性。下面我们将讨论这些目标的实现程度，如果没有实现，读者可以做些什么来弥补这些局限性。

一个复杂而技术性的主题本质上很难访问。在网上搜索“贝叶斯非参数”会产生大量关于这个主题的教程、视频、网页和博客。这本文献太多了，对于一个想提前学习DP的初学者来说，它们是压倒性的，甚至可能会让人瘫痪。他/她可以从哪里开始？迎合新手的教程往往主要侧重于直观的介绍，使用与我们的图3。统计的严谨性保持在与我们相似的复杂程度方程式（10），通常省略中间推导。解释了交换性等重要概念。然而，对于初学者如何将这些主题扩展到新模型，通常只有有限的讨论。我们需要对BNP有更深入的了解，例如。，方程式（12）–(13)并详细解释了如何将这些概念付诸实践。一些介绍超越了这些基本解释，试图填补文献中的第一个空白。更高级的展览（例如。，尼尔（2000）;拉斯穆森（2000）;布罗德里克（2017）;格什曼律师事务所（2012）)覆盖这些。它们帮助读者一步一步地了解统计的严谨性。然而，许多关键推导要么被省略，要么没有足够详细的解释，要么只是假设理解。这是另一个重要的差距，也许也是行为科学家特别难以克服的障碍。还有其他更具挑战性的介绍引用了度量理论(2017年;Orbanz&Teh，2011年;Jara，2017年;Müller等人，2015年;Hjort等人，2012年;Ghoshal&van der Vaart，2015年)非技术人员可能无法理解。

本教程有助于填补所有这些空白。统计的严谨性逐步提高，并辅以直观的解释。基本概念被赋予了统计严谨性，尤其是关于附录A一些推导一开始可能看起来很吓人，例如在Dirichlet密度上的积分方程式（A.6）然而，它们毕竟不是不可穿透的。我们表明，没有必要通过整合进行劳动。这些推导为CRP隐喻提供了一个透彻的解释。初学者可以通过了解方程是如何推导出来的，以及如何在计算机程序中实现这些方程，从而超越基础知识。我们希望我们已经实现了目标1和目标2，使读者满意。对于技术专家来说，我们的详细推导和解释可能显得多余，甚至可能迂腐。然而，我们希望他们同意我们的观点，即这样的教程有助于解决文献中的差距。

更深入的理解也可以帮助读者更好地欣赏其他教程。例如，一个初学者首先阅读Austerweil等人（2015）本章可能并没有完全理解他们对贝叶斯非参数的定义。BNP是从一般统计角度来看的，其中“非参数”指的是不限于某些固定参数类（例如，法线的有限混合）的模型。本教程显示BNP将有限混合模型扩展到理论上无限个参数。将BNP应用于特定数据集可能会产生一组有限的参数，从而使BNP模型的复杂性随着数据的增加而增加。

目标3通过解释R计划来实现(布罗德里克，2017年)一行接一行。在这样做的过程中，我们将方程与相应的计算算法连接起来，并应用程序恢复中假设数据的已知簇图1编程迫使学习者非常清楚，因为所有步骤都必须准确明确地表达出来。我们认识到CRP的轻微并发症。随着模拟的进行，集群可能会变为空。必须采取措施考虑到这一点。本练习通过DP previor提供了对BNP的关键见解。值得注意的是，新方法只能从两个基本成分发展而来：1）方程（12）-(13); 以及2）根据模型的先验和后验预测分布。BNP并不像最初看起来那样不可逾越。当然，一个完全安全的计算工具需要大量的工作。我们的目标3没有提供完全安全的工具。然而，我们希望读者看到示例R程序会感到鼓舞，该程序可以很容易地开始成形。所需的技能并不比编写循环和调用库存函数以从共轭分布中采样复杂得多。有了这一洞察力（只需要两个要素），定量心理学的高级研究生可以将简单的程序扩展到完整的解决方案，或者将其扩展到一种全新的聚类技术。

我们知道有几种强大的DP计算工具可用(Jara等人，2011年;Pedregosa等人，2011年;Eisentein，2012年;Christensen等人，2011年)。一个只想快速解决问题的研究人员不必费心编写吉布斯采样代码。然而，我们认为编程练习有助于加深对BNP的总体理解。本教程没有提供任何新颖的方法。另一个限制是，我们没有研究协方差矩阵的共轭逆Wishart分布中的已知困难(Gelman等人，2003年第19.2节）。感兴趣的读者可以在Lewandowski等人（2009）尽管存在这些局限性，但我们填补了文献中的重要空白，这些空白限制了BNP方法在行为科学中的广泛接受。我们希望学生和他们的教授会发现我们的教学法有助于他们的研究和教学，尤其是在训练方法学家时更容易理解DPMM等复杂概念。他们可能能够将基本技能和深入理解扩展到新方法的开发中。

DP混合模型只是一大类灵活的BNP方法的一个例子。本教程为行为科学家提供了进入BNP的途径，他们将知识和数学推导放在首位，使BNP能够解释个体差异。

​

集锦

  • 本文提供了狄里克莱过程发展所必需的数学推导。
  • DP建模中复杂和抽象数学的易用教学法。
  • R中的一个可公开访问的计算机程序逐行解释。
  • 晚期癌症患者的心理治疗效果因DP混合模型确定的心理社会特征而异。

致谢

这项工作得到了美国国立卫生研究院（NIH）向斯隆-凯特琳纪念癌症中心（Memorial Sloan Kettering Cancer Center）拨款P30 CA008748的支持。

我们感谢塔玛拉·布罗德里克（Tamara Broderick）对R程序的帮助，该程序可从以下网址下载：https://people.csail.mit.edu/tbroderick/tutorial_2017/mitll.html

我们感谢《教程》编辑Maarten Speekenbrink和Joseph Houpt的建设性意见。我们还感谢另外两位匿名评论员的宝贵意见。

附录A。技术细节推导

本附录总结了DPMM中基本方程的数学推导，其他教程中通常会省略这些推导。总体计划是用细节补充正文。我们表明，可以使用介绍性文本中的基本贝叶斯统计来严格解释DP。我们围绕着度量理论工作，这往往让初学者难以理解。希望心理学定量领域的研究生在使用测量理论处理主流文献之前，可以使用本教程首先对DP有一个清晰的了解，

本附录分为4小节。我们首先描述高斯混合中的数据似然函数，因为它将用于证明其他方程。然后我们从聚类平均数的后验分布开始，提供详细的推导来证明所有基本方程(附录A.2)。按照正文的流程，我们接下来逐个解释其他方程式。我们解决后精度问题(附录A.3)和其他详细信息。我们以中餐厅流程中潜在的集群成员指标结束(附录A.4).

附录A.1。高斯混合中的数据似然函数

正态分布的概率密度是其平均值的函数μ和标准偏差σ:

12πσ2经验(−年我−μ2σ2),

当乘以嵌套在每个潜在簇内的独立观测值时，我们得到数据似然函数：

我(年∣μk个,σk个2)∝1(σ2)n个k个经验(−∑我=1n个k个(年我−μk个)22σk个2),

其中，如果标准偏差与精度互换τk个= 1/σ2，然后我们得到

我(年∣μk个,τk个)∝τn个k个∕2经验(−τk个∑我=1n个k个(年我−μk个)22).
（A.1）

本节使用聚类精度的数据似然函数来更新其他参数的先验值。注意，上面的符号是基于一元高斯密度的，因为使用它要简单得多。然而，这同样适用于多元高斯密度。

附录A.2。簇的后验密度平均值μk个

数据的可能性方程式（A.1）当与潜在聚类均值的高斯先验相结合时，N个(μ0,τ0)，得出平均值的后验分布μk个.

第页(μk个∣年)∝第页(年∣μk个)第页(μk个)∝经验(−τk个∑我=1n个k个(年我−μk个)22)×经验(−τ0(μk个−μ0)22),
（A.2）

其中，第一项是在k个聚类平均值和精度。第二项是概率μk个事先进行评估。请注意

∑1n个(年我−μk个)2=∑1n个(年我−年‒+年‒−μk个)2=∑1n个(年我−年‒)2+n个(年‒−μk个)2因为2(年‒−μk个)∑1n个(年我−年‒)=0,

因此方程式（A.2）变为（下标k个删除以使方程式更容易理解）：

第页(μ∣年)∝经验(−12[τ(∑我n个(年−年‒)2+n个(年‒−μ)2)+τ0(μ2−2μ0μ+μ02)])∝经验(−12[τn个(年-2−2年‒μ+μ2)+τ0(μ2−2μ0μ)])∝经验(−12[−2τn个年‒μ+τn个μ2+τ0μ2−2τ0μ0μ])∝经验(−12[μ2(τn个+τ0)]−μ(2τn个年‒+2τ0μ0))∝经验(−12[μ2(τn个+τ0)+μ(τn个年‒+τ0μ0)).
（A.3）

如果我们称之为后验均值和精度μ第页和τ第页，分别写入

τ第页=τn个+τ0μ第页=1τ第页(τn个年‒+τ0μ0),
（A.4）

我们把这些插回最后一行方程式（A.3），我们得到

第页(μ∣年)∝经验(−12μ2τ第页+μμ第页τ第页).

添加到指数−12μ第页2τ第页，这是一个常量μ我们看到了

第页(μ∣年)∝经验(−12μ2τ第页−12μ第页2τ第页+μμ第页τ第页).

重新排列术语，指数现在是

经验(−12(μ2−2μμ第页+μ第页2)τ第页),

因此

第页(μ∣年)∝经验(−12(μ−μ第页)2τ第页).

因此，聚类平均值的后验分布是一个高斯分布，在方程式（A.4）。重要的是，这些参数与前面显示的参数完全相同方程式（4）.

我们有一种常见的形式，即后验精度等于先验精度加上数据精度，后验平均值等于先验平均值和数据平均值的加权平均值(Gelman等人，2003年第2.6节）。

附录A.3。聚类精度的后验密度τk个

聚类精度是给定Gamma先验的形状参数γ和率参数β,第页(τk个∣γ,β)=G公司(γ,β)∝τk个γ−1经验(−βτ)。中的数据可能性方程式（A.1）组合以获得

第页(τk个∣μ,年,γ,β)∝τk个γ−1经验(−βτk个)×τk个n个k个∕2经验(−τk个∑我=1n个k个(年我−μk个)22)=τk个γ+n个k个∕2−1经验(−τk个(β+∑我=1n个k个(年我−μk个)22)).
（A.5）

因此，合并数据似然后，聚类精度遵循具有形状参数的伽马后验分布γ+n个k个/2和速率参数β+∑我=1n个k个(年我−μk个)22.

附录A.4。潜在簇的后向密度

这里我们计算出了方程式（7），首先将混合比例积分πk个输出：

第页(c（c）1,…,c（c）k个∣α)=∫第页(c（c）1,…,c（c）k个∣π1,…,πk个)第页(π1,…,πk个)d日π1⋯d日πk个=∫∏k个=1K（K）πk个n个k个︸等式.6Γ(α)Γ(α∕K（K）)k个∏k个=1K（K）πk个α∕K（K）−1︸等式.3d日π1⋯d日πk个=Γ(α)Γ(α∕K（K）)k个∫∏k个=1K（K）πk个n个k个+α∕K（K）−1d日π1⋯πk个,
（A.6）

其中两个产品术语∏k个=1K（K）πk个n个k个和∏k个=1K（K）πk个α∕K（K）−1与指数组合成一个πk个简单地加在一起；Γ(α)Γ(α∕K（K）)k个是一个常数，因此放在积分之外（因为它不涉及πk个)。然而，∫∏k个=1K（K）πk个n个k个+α∕K（K）−1d日π1⋯d日πk个很难集成，因为它涉及K（K）条款。但请注意它与狄利克雷密度有多相似。如果我们在它前面加上一个常数

Γ(∑k个=1K（K）(n个k个+α∕K（K）))∏k个=1K（K）Γ(n个k个+α∕K（K）)∏k个=1K（K）πk个n个k个+α∕K（K）−1,

那么这个新公式就是带参数的Dirichlet分布的密度n个k个+α/K.

如果我们将这个Dirichlet密度积分到πk个,

∫Γ(∑k个=1K（K）(n个k个+α∕K（K）))∏k个=1K（K）Γ(n个k个+α∕K（K）)∏k个=1K（K）πk个n个k个+α∕K（K）−1d日π1⋯d日πk个,

我们知道它精确地积分为1，因为一个基本事实是，适当的概率密度必须积分为1。因此，

Γ(∑k个=1K（K）(n个k个+α∕K（K）))∏k个=1K（K）Γ(n个k个+α∕K（K）)∫∏k个=1K（K）πk个n个k个+α∕K（K）−1d日π1⋯d日πk个=1.

注意，上面的公式可以分为两部分，一个常数项Γ(∑k个=1K（K）(n个k个+α∕K（K）))∏k个=1K（K）Γ(n个k个+α∕K（K）)和一个积分∫∏k个=1K（K）πk个n个k个+α∕K（K）−1d日π1⋯d日πk个。因为它们相乘后正好为1，所以后一项必须是前者的倒数。因此，

∫∏k个=1K（K）πk个n个k个+α∕K（K）−1d日π1⋯d日πk个=∏k个=1K（K）Γ(n个k个+α∕K（K）)Γ(∑k个=1K（K）(n个k个+α∕K（K）)).

我们完全避免了复杂的集成。把这个插回方程式（A.6），我们得到

第页(c（c）1,…,c（c）k个∣α)=Γ(α)Γ(α∕K（K）)k个∫∏k个=1K（K）πk个n个k个+α∕K（K）−1d日π1⋯d日πk个=Γ(α)Γ(α∕K（K）)k个∏k个=1K（K）Γ(n个k个+α∕K（K）)Γ(∑k个=1K（K）(n个k个+α∕K（K）)),而且因为Γ(α∕K（K）)k个=∏k个=1K（K）Γ(α∕K（K）),还有Γ(∑k个=1K（K）(n个k个+α∕K（K）))=Γ(n个+α),我们得到第页(c（c）1,…,c（c）k个∣α)=Γ(α)Γ(n个+α)∏k个=1K（K）Γ(n个k个+α∕K（K）)Γ(α∕K（K）),
（A.7）

这正是方程式（7）。请注意Γ(∑k个=1K（K）(n个k个+α∕K（K）))=Γ(n个+α)因为单个集群的大小，n个k个，加上数据点总数n个此外，求个人α/K结束K（K）简单收益率α.

现在我们转向方程式（8）计算出省略的推导。我们需要潜在集群成员的概率我这个人认为所有其他人潜在的阶级成员身份是固定的和已知的。这可以追溯到CRP隐喻，其中所有之前的我−1名客户已就座，然后我这位顾客走进餐厅，需要一张桌子。该客户可以坐在已经有人的桌子旁，也可以坐在新桌子旁。下面的等式规定了如何将客户分配给已占用的表k个.

方程式图像
（A.8）

在第二行中，相同的术语相互抵消，重新排列术语后，我们得到最后一行。它基于伽玛函数Γ的递推性质(x个+ 1) =x个Γ(x个)，所以

∏k个=1K（K）Γ(n个k个+α∕K（K）)Γ(n个k个,−我+α∕K（K）)=(n个k个,−我+α∕K（K）)Γ(n个k个,−我+α∕K（K）)Γ(n个k个,−我+α∕K（K）)=n个k个,−我+α∕K（K）1.

更直观的解释是，因为Γ(x个)=============================================================================================================================================================(x个−1）！因此

Γ(n个+α−1)∕Γ(n个+α)=(n个+α−2)!∕(n个+α−1)!=1n个+α−1.

通过求解一个看似无法穿透的方程中的复杂项，一个简单的模式出现在方程式（A.8）-概率我第个客户坐在桌旁k个主要受以下因素的影响n个k个,−我，每个k个第个表。

附录B.R代码

1 #2#Dirichlet过程混合物模型示例。3#上次修改时间：2018年4月19日，由李月林（Yuelin Li）根据以下示例4#塔玛拉·布罗德里克（Tamara Broderick），其原始源代码位于：5 #https://raw.githubusercontent.com/tbroderick/bnp_tutorial6#/2016mlss/ex7_采样器。R（右）7#####以下是原始源代码文件中的标题######用于CRP-高斯混合模型的8#A吉布斯采样器Neal 2000中的9#算法310#版权所有（C）2015，Tamara Broderick11个#www.tamarabroderick.com12 #13#此程序是自由软件：您可以重新发布和/或修改14#根据GNU通用公共许可证的条款15#自由软件基金会，许可证版本3，或16#（由您选择）任何更高版本。17 #18#发布这个程序是为了希望它会有用，19#但无任何担保；甚至没有20#针对特定目的的适销性或适用性。请参阅21#GNU通用公共许可证了解更多详细信息。22个#23#您应该已经收到GNU通用公共许可证的副本24#与此程序一起使用。如果没有，请参见<-ttp://www.gnu.org/licenses（网址：www.gnu.org/licenses）/>.25#用于处理多元正态分布26 #27#####以下由YL修改。变量名称更改自28号原始符号与我们在教程中所能找到的符号相匹配。29 crp_gibbs<-函数（数据，alpha=0.01，mu0，sigma0，simma_y，c_init，maxIters=1000）30 {31#数据：数据点的N x D矩阵32#小alpha鼓励更少的集群33#α<-0.0134#σy：假设已知的数据y的测量误差，对于所有簇都相同。35#西格玛y<-diag（data_dim）*136#mu0，sigma0：未知平均值mu周围的先验均值和方差37#mu0<-矩阵（rep（0，data_dim），ncol=data_dim，byrow=TRUE）38#西格玛0<-diag（data_dim）*3^239#c_init：将数据点初始分配给集群40 #41要求（mvtnorm）42#数据点维度43 data_dim<-ncol（数据）44#数据点数45 N<-nrow（数据）46 #####47#前48 #####49#未知平均mu的先验精度。50 tau0<-solve（sigma0）#$mu$上的先验精度，先验协方差的逆51#集群特定精度，假设已知，所有集群假设为52#共享y~N（mu，sigma_y）的相同测量误差。53τy<-解（σy）54#初始化CRP Gibbs采样器55 z<-cinit#初始集群成员分配56 n_k<-as.vector（表（z））#每个簇的初始数据计数，等式（4）。57 Ncluster<-length（n_k）#初始簇数58 ##59号中餐厅工艺（CRP）吉布斯采样器开始60 ##61 res<-矩阵（NA，nrow=N，ncol=maxIters）#集群成员存储62 pb<-txtProgressBar（最小值=0，最大值=maxIters，样式=3）63 for（iter in 1:maxIters）{#maxIters还可以防止无休止循环64 for（n in 1:n）{#一次一个数据点（客户）65#当第n位顾客进入中餐厅时，我们需要首先66#取消分配他/她的初始集群成员，然后使用等式（12）67#[已占用表]和（13）[新表]计算68号更新了表赋值的概率。69 ci<-z[n]#第n个人表赋值是什么？70 n_k[ci]<-n_k[ci]-1#从表中删除第n个人71#如果删除第n个人时表格变空，72#然后删除该表/集群。73如果（n_k[c_i]==0）74 {75 n_k[c_i]<-n_k[Nclust]#替换此空集群的最后一个集群76 loc_z<-（z==Nclust）#谁在最后一个集群中？77 z[loc_z]<-ci#向上移动它们以填充刚刚清空的集群78 n_k<-n_k[-Nclust]#取出最后一个集群，现在为空79 Nclust<-Nclust-1#将集群总数减少180 }81 z[n]<-−1#确保z[n]不会被计算为集群#####82个##83#现在，我们准备通过公式（12）和（13）更新表赋值。84 ##85#集群的日志概率，添加之前未占用的表86 logp<-rep（北美，Nclust+1）87#回路在已占用的表1:J上，并按式（13）计算pr。88代表（1中的c_i：Nclust）{89 tau_p<-tau0+n_k[c_i]*tau_y#根据等式（4）的聚类精度90 sig_p<-求解（tau_p）#簇方差，tau_c的逆91#找到该簇中的所有点92 loc_z<-其中（z==c_i）93#将该簇中的所有点相加94如果（长度（loc_z）>1）{95 sum_data<-colSums（数据[z==c_i，]）}96其他{97汇总数据<-data[z==c_i，]98   }99#100#我们需要使用每个预测分布101#占据了桌子，预测坐在那里的下一位客户。102 #103#使用等式（4）找出104#集群表示（y*n_k*tau_j+mu0*s0）/tau_p，以及105#然后使用等式（11）中y_j的预测分布106从c-i预测新数据值ci。107 #108平均值p<-sig_p%*%（tau_y%*%sum_data+tau0%*%t（mu0））109 logp[ci]<-log（n_k[ci]）+110 dmv范数（数据[n，]，平均值=平均值_p，σ=sig_p+σ_y，log=真）}111 #112#我们完成了对已经占用的表的循环。接下来，我们使用113#等式（12）计算之前114#空闲，“新”表。从本质上讲，它是预先预测的115#DP分配。116 #117 logp[Nclust+1]<-log（alpha）+118 dmvnorm（数据[n，]，平均值=mu0，西格玛=sigma0+sigma_y，对数=TRUE）119#将非正规对数概率转换为概率120最大logp<-最大（logp）121 logp<-logp-最大logp122个loc_probs<-exp（logp）123个loc_probs<-loc_probs/sum（loc_props）124#抽取此新客户应属于哪个集群的样本125 newz<-sample（1:（Nclust+1），1，replace=TRUE，prob=loc_probs）126#必要时生成新集群127如果（newz==Nclust+1）{128 n_k<-c（n_k，0）129 Nclust<-Nclust+1130   }131 z[n]<-新z132 n_k[newz]<-n_k[newz]+1#更新集群n_k133 }134 setTxtProgressBar（pb，iter）#在每个iter之后更新文本进度条135 res[，iter]<-z#N个观测的聚类成员136 }137 close（pb）#关闭文本进度条138不可见（res）#返回结果，N by maxIters矩阵139 }

脚注

出版商免责声明：这是一份未经编辑的手稿的PDF文件，已被接受出版。作为对客户的服务，我们正在提供这份早期版本的手稿。手稿将经过编辑、排版和校对，然后才能以最终形式出版。请注意，在制作过程中可能会发现错误，这可能会影响内容，所有适用于该杂志的法律免责声明都适用。

工具书类

  • Aldous DJ（1985）。Hennequin PL（Ed.）中的互换性和相关主题，埃科尔·德·埃特塞·德·圣·弗洛尔第十三讲数学笔记（第1-198页）。柏林：斯普林格。[谷歌学者]
  • Anderson JR（1991）。人类分类的适应性.心理回顾,98, 409–429. doi:10.1037/0033-295X.98.3.409。[交叉参考][谷歌学者]
  • 安东尼亚克C（1974）。Dirichlet过程的混合及其在贝叶斯非参数问题中的应用.统计年鉴,2, 11521174.[谷歌学者]
  • Austerweil JL、Gershman SJ、Tenenbaum JB和Griffiths TL（2015）。贝叶斯认知模型的结构和灵活性牛津计算和数学心理学手册（第187–208页）。[谷歌学者]
  • Austerweil JL和Griffiths TL（2013）。构建柔性特征表示的非参数贝叶斯框架.心理回顾,120，817–51网址：https://www.ncbi.nlm.nih.gov/pubmed/24219850.doi:10.1037/a003419[公共医学] [交叉参考][谷歌学者]
  • Bernardo JM和Smith AFM（2000年）。贝叶斯理论英国奇切斯特：John Wiley&Sons，Ltd。[谷歌学者]
  • Blackwell D和MacQueen J（1973）。通过Polya urn方案的Ferguson分布.统计年刊,1, 353–355.[谷歌学者]
  • Broderick T（2017）。非参数贝叶斯教程。2017年7月25日检索自https://people.csail.mit.edu/tbroderick/tutorial_2017_mitll.html.麻省理工学院林肯实验室。[谷歌学者]
  • Christensen R、Johnson W、Branscum A和Hanson TE（2011年）。贝叶斯思想与数据分析佛罗里达州博卡拉顿：CRC出版社。[谷歌学者]
  • Collins AGE和Frank MJ（2013年）。对学习的认知控制：创建、聚类和概括任务集结构.心理回顾,120, 190–229. doi:10.1037/a0030852。[PMC免费文章][公共医学] [交叉参考][谷歌学者]
  • Eisentein J（2012）。Dpmm公司.https://github.com/jacobeisenstein/DPMM网站.
  • Escobar M和West M（1995）。混合贝叶斯密度估计与推理.美国统计协会杂志,90, 577–588.[谷歌学者]
  • Ferguson T（1973）。一些非参数问题的贝叶斯分析.统计年刊,1, 209–230. doi:10.1214/aos/1176342360。[交叉参考][谷歌学者]
  • Gelman A、Carlin JB、Stern HS和Rubin DB（2003年）。贝叶斯数据分析纽约：查普曼和霍尔[谷歌学者]
  • Gershman SJ和Blei DM（2012年）。贝叶斯非参数模型教程.数学心理学杂志,56, 1–12. doi:10.1016/j.jp.2011.08.004。[交叉参考][谷歌学者]
  • Ghoshal S和van der Vaart A（2015）。非参数贝叶斯推理基础英国剑桥：剑桥大学出版社。[谷歌学者]
  • Hjort NL、Holmes C、Mller P和Walker SG（2012年）。贝叶斯非参数纽约：剑桥大学出版社。[谷歌学者]
  • Hogg RV和Craig AT（1995年）。数理统计导论Englewood Cliffs，新泽西州：普伦蒂斯·霍尔。[谷歌学者]
  • Houpt JW（2018）。高斯过程.http://www.wright.edu（赖特）/~约瑟夫·霍普特/。访问日期：2018年6月10日。
  • Jara A（2017）。Dirichlet过程和相关模型的理论和计算：综述.国际近似推理杂志,81, 128–146. doi:10.1016/j.ijar.2016.11.008。[交叉参考][谷歌学者]
  • Jara A、Hanson TE、Quintana FA、Müller P和Rosner GL（2011年）。DPpackage:r中的贝叶斯半参数和非参数建模.统计软件杂志,40, 1–30. 网址：<GotoISI>：//WOS:000289228600001。[PMC免费文章][公共医学][谷歌学者]
  • Jasra A、Holmes CC和Stephens DA（2005年）。马尔可夫链蒙特卡罗方法与贝叶斯混合建模中的标签切换问题.统计科学,20, 50–67. 网址：<GotoISI>：//WOS:00022990600006。doi:doi 10.1214/08834230500000016。[交叉参考][谷歌学者]
  • Karabatsos G（2017）。用于回归分析和密度估计的贝叶斯非参数（和参数）混合模型的菜单驱动软件包.行为研究方法,49, 335–362. 网址：10.3758/s13428-016-0711-7。doi:10.3758/s13428-016-0711-7。[公共医学] [交叉参考] [交叉参考][谷歌学者]
  • Karabatsos G和Walker SG（2008年）。测试等值的贝叶斯非参数方法.心理测量学,74，211网址：10.1007/s11336-008-9096-6。doi:10.1007/s11336-008-9096-6。[交叉参考] [交叉参考][谷歌学者]
  • Karabatsos G和Walker SG（2009年）。基于贝叶斯非参数的一致性心理测量模型.英国数学与统计心理学杂志,62，1–20。地址：10.1348/000711007X246237。[公共医学] [交叉参考][谷歌学者]
  • 库鲁克利德斯一世（2018），贝叶斯非参数.kourouklides.wikia.com/wiki/Baysian_Nonparametrics。访问日期：2018年6月4日。
  • Lee PM（2012）。贝叶斯统计：简介（第四版）。英国西苏塞克斯：John Wiley&Sons。[谷歌学者]
  • Lewandowski D、Kurowicka D和Joe H（2009年）。基于藤蔓和扩展洋葱法的随机相关矩阵生成.多元分析杂志,100, 1989–2001. doi:10.1016/j.jmva.2009.04.008。[交叉参考][谷歌学者]
  • Li Y、Lord-Bessen J、Shiyko M和Loeb R（2018年）。贝叶斯潜在类分析教程.多元行为研究,35, 1–22. doi:10.1080/0273171.2018.1428892。[PMC免费文章][公共医学] [交叉参考][谷歌学者]
  • Liew SX、Howe PDL和Little DR（2016年）。平均法在语境效应研究中的适用性.心理通报与评论,23, 1639–1646. doi:10.3758/s13423-016-1032-7。[公共医学] [交叉参考][谷歌学者]
  • Müller P、Quintana FA、Jara A和Hanson T（2015）。贝叶斯非参数数据分析统计学中的斯普林格系列瑞士：施普林格国际出版公司。doi:10.007/978-3-319-18968-0。[交叉参考][谷歌学者]
  • Navarro DJ、Griffiths TL、Steyvers M和Lee MD（2006年）。使用Dirichlet过程建模个体差异.数学心理学杂志.50，pp.doi:10.1016/j.jmp.2005.11.006。[交叉参考][谷歌学者]
  • Neal RM（2000）。Dirichlet过程混合模型的马尔可夫链抽样方法.计算与图形统计杂志,9, 249–265. 网址：http://www.jstor.org/stable/1390653.doi:10.2307/1390653。[交叉参考][谷歌学者]
  • Novick MR和Jackson PH（1974年）。教育与心理研究的统计方法.纽约：麦格劳山。[谷歌学者]
  • Orbanz P（2018）。贝叶斯非参数教程.http://www.stat.columbia.edu/~porbanz/npb-tutorial.html。访问日期：2018年6月4日。
  • Orbanz P，&Teh TY（2011年）。贝叶斯非参数模型In机器学习百科全书施普林格美国；网址：http://www.stat.columbia.edu/~porbanz/papers/porbanz_OrbanzTeh2010_1.pdf.[谷歌学者]
  • 佩德雷戈萨F、瓦罗佐G、格拉姆福特A、米歇尔·V、蒂里昂B、格里塞尔O、布隆德尔M、普雷滕霍弗P、韦斯R、杜堡V、范德普拉斯J、帕索斯A、库纳波D、布鲁彻M、佩罗特M和公爵内伊E（2011）。Scikitslearn:Python中的机器学习.机器学习研究杂志,122825–2830.[谷歌学者]
  • 拉斯穆森CE（2000）。无限高斯混合模型，Solla S，L.T.K.，&Müller K（Eds.），神经信息处理系统研究进展 12（第554–560页）。剑桥：麻省理工学院出版社。[谷歌学者]
  • Richardson S和Green PJ（1997年）。未知组分数混合物的贝叶斯分析.英国皇家统计学会学报B辑方法学,59, 731–758. doi:10.1111/1467-9868.00095。[交叉参考][谷歌学者]
  • Sethuraman J（1994）。dirichlet先验的构造性定义.中国统计局,4, 639–650. 网址：http://www3.stat.sinica.edu.tw/statistica/j4n2/j4n216/j4n21.htm.[谷歌学者]
  • 斯蒂芬斯M（2000）。混合模型中标签切换的处理.英国皇家统计学会期刊B辑统计方法,62, 795–809. doi:10.1111/1467-9868.00265。[交叉参考][谷歌学者]
  • Teh YW（2017）。Sammut C和Webb GI（编辑）中的Dirichlet过程，机器学习和数据挖掘百科全书（第361-370页）。马萨诸塞州波士顿：Springer US。doi:10.1007/9781-4899-7687-1_219。[交叉参考][谷歌学者]