Skip to content

Latest commit

 

History

History
445 lines (305 loc) · 46 KB

File metadata and controls

445 lines (305 loc) · 46 KB

前言:狄利克雷过程是狄利克雷分布的无限维的一般化,可以被用作设置未知分布的先验。此外,这些未知分布的密度函数可以将有限分量密度模型拓展为无限分量密度模型。

狄利克雷过程模型

狄利克雷过程定义

狄利克雷分布

​ 在概率统计中,狄利克雷分布(Dirichlet distribution)是一系列由正实数向量$\alpha$参数化的连续型多元概率分布,通常表示为$Dir(\alpha)$。它是beta分布的多元推广,因此,它也被称为多元beta分布(multivariate beta distribution/MBD)。狄利克雷分布通常用作贝叶斯统计中的先验分布,实际上狄利克雷分布是分类分布(Categorical distribution)和多项式分布(Multinomial distribution)的共轭先验。

狄利克雷分布的无限维推广是狄利克雷过程

狄利克雷过程

简介

​ 在概率论中,狄利克雷过程是一系列随机过程,其实现是概率分布。它通常用于贝叶斯推理中描述随机变量分布的先验,即随机变量的分布是依据一种或另一种特定分布的可能性。

​ **狄利克雷过程(Dirichlet process,DP)**由基分布H和称为浓度参数的正实数$\alpha$指定。其中,基分布是过程的期望值。然而,即使基分布是连续的,从狄利克雷过程中得出的分布几乎肯定是离散的,浓度参数则制定了这种离散化的强度。

​ 既然狄利克雷过程是狄利克雷分布的推广,与狄利克雷分布是分类分布的共轭先验一样,狄利克雷过程是无限非参数离散分布的共轭先验。狄利克雷过程的一个特别重要的应用是作为无限混合模型中的先验概率分布。

正式定义

​ 给定样本空间$\Omega$,划分成可测度的子集$B_1,B_2,...,B_k$。给定一个基本概率分布P和一个正实数$\alpha$,**狄利克雷过程$DP(P,\alpha)$**是一个随机过程,其样本路径(即从过程中抽取的随机变量的无限序列)是$\Omega$上的概率分布。因此,对于$\Omega$的任何可测有限分区${B_i}_{i=1}^n$,有:

如果$P\sim DP(\alpha P_0)$,那么$(P(B_1),...,P(B_n))\sim Dir(\alpha P_0(B_1),...,\alpha P_0(B_n))$。 其中,Dir表示狄利克雷分布,$P\sim D$表示随机变量P服从分布D。$P_0$是提供对P的初始猜测的基本概率度量。$\alpha$表示控制P向$P_0$收缩的程度的先验浓度参数。

​ 狄利克雷过程是关于分布的分布,即狄利克雷过程的每个采样本身即是一个随机分布,该随机过程之所以被称为狄利克雷过程,是因为其任意有限维边缘分布均是狄利克雷分布。这和高斯过程类似,其任意有限维边缘分布均是高斯分布。

其他定义

​ 除了正始定义外,狄利克雷过程还有几种等效的定义。比如中餐厅过程(Chinese restaurant process,CRP)、破棒过程(Stick-breaking process)和Polya urn过程。这些将在后面的部分详细探讨。

狄利克雷过程混合模型

狄利克雷过程三种构造

  1. Stick-breaking构造(破棒过程)

    ​ 基于相互独立的变量序列$(\beta_k){k=1}^{\infty}$和$(\phi_k){k=1}^{\infty}$的Stick-breaking构造:

$$ \beta_k|\alpha_0, G_0\sim Beta(1,\alpha_0), \phi_k|\alpha_0,G_0\sim G_0 $$ 定义随机概率分布G如下: $$ \pi_k =\beta_k\prod_{l=1}^{k-1}(1-\beta_l),G=\sum_{k=1}^{\infty}\pi_k\delta_{\phi_k} $$ 这样的分布G是服从狄利克雷过程$DP(\alpha_0,G_0)$的一个随机分布。其中,$\delta_{\phi}$是$\phi$点的概率测度。$\pi =(\pi_k){k=1}^{\infty}$以概率1满足$\sum{k=1}^{\infty}\pi_k =1$.因此,可以将π视为关于正整数的随机概率分布,用$\pi \sim GEM(\alpha_0)$表示权重系数的构造关系。

​ 破棒过程引出$P\sim DP(\alpha P_0)$是通过令: $$ P(·)=\sum_{h=1}^{\infty}\pi_h\delta_{\theta_h}(·),,,,\pi_h=V_h\prod_{l<h}(1-V_l),,,,V_h\sim Beta(1,\alpha),,,,\theta_h\sim P_0 $$ 其中,$δ_θ$表示所有概率质量都在$\theta$点的退化分布,原子$(θ_h)^∞_h=1$是独立于基础分布$P_0$生成的,$π_h$是原子$θ_h$的概率质量,这些概率质量生成一个所谓的破棒过程,保证权重和为1。

​ 为了描述破棒过程,我们从一根单位长度的棒开始,它的长度表示分配给所有原子的总概率。我们首先分离一个长度为$V_1$的随机片段,其长度由$Beta(1,α)$分布生成,并将这个$π1=V1$概率权值分配给随机生成的第一个原子$θ_1∼P_0$。然后剩下$1−V_1$的棒子被分配给其他原子。我们将$1−V_1$棒的比例$V_2∼Beta(1,α)$分离出来,并将概率$π_2=V_2(1−V_1)$分配给第二个原子$θ_2∼P_0$。随着我们的继续,该棒变得越来越短,因此分配给较高索引原子的长度随机减小,而减小的速率取决于DP精度参数α。因为$E(V_h)=\frac{1}{1+\alpha}$,α的值接近0会导致第一对原子权重很高,其余的原子被分配的概率很小。

image-20211020144911004

图26.1 来自不同精度参数α设置的狄利克雷过程的截断表示的样本。

图26.1显示了对应标准正态分布和不同α值的$P_0$的破棒过程的实现。从图中可以明显看出,DP不适合作为数据分布的直接先验,特别是在数据是连续的情况下。对于连续数据,每个新主题需要一个新的原子,因此需要一个大的α值,这意味着每个原子的权重接近零,因此p实现的概率很小。在$α→∞$的极限中,我们得到$y_i∼P_0$,因此对于大的α,以及观测中没有联系,DP先验可以有效地从参数基分布中提取数据。

  1. Polya urn scheme构造

DP先验分布作为数据分布的直接模型的失败并不意味着它在应用程序中并不有用。相反,DP更合适地用作未知混合物分布的先验。为了简单起见,再次关注密度估计的情况,密度的一般核混合模型可以指定为 $$ f(y|P)=\int \kappa (y|θ)dP(θ) $$ ,其中$\kappa$(·|θ)是一个核,θ是包括位置的可能的尺度参数,P是一个混合度量。在P是有限的k个原子的质量的离散的特殊情况下,得到高斯过程章节中讨论的有限混合模型。在无限核混合模型中,为未知混合P选择先验P∼πP,其中P表示(Ω,B)上所有概率测度的空间,πP表示该空间上的先验。混合测量的先验通过该式中的积分映射得到密度f(y)的先验。如果选择πP对应于DP先验,那么就得到一个DP混合模型。并且有,P上的DP先验: $$ f(y)=\sum_{h=1}^{\infty}\pi_h\kappa(y|\theta_h^*) $$ 其中$\pi \sim stick(\alpha)$是一个简记符号,表示概率权值是从具有参数$\alpha$和$\theta_h \sim P_0$独立同分布的DP破棒过程中采样得出的。$h=1,...,\infty$.

​ 上式类似于高斯过程章节中考虑的有限混合模型,但重要的区别是混合成分的数量被设置成无穷大。然而,这并不意味着样本中的被试者占据了无限多个成分;相反,模型允许通过添加受试者时引入新的混合成分来实现灵活性。考虑如下的层次结构规范: $$ y_i\sim \kappa(\theta_i),,,\theta_i\sim P,,,P\sim DP(\alpha P_0) $$ 该公式等价于从无限混合模型中采样$y_i$。一个关键问题是如何在这种DP混合物(DPM)下进行后验计算?因为混合测度P具有无穷多个参数,在P给定的$y^n=(y_1,...,y_n)$的后验有一个简单的形式,不再有联合共轭。

​ 解决这个问题的一个聪明的方法是边缘化P,以获得对特定主题的参数θn=(θ1,……,θn)的诱导先验分布。特别地,边缘化P,我们得到了Polya urn预测规则: $$ p(\theta_i|\theta_1,...,\theta_{i-1})\sim (\frac{\alpha}{\alpha+i-1})P_0(\theta_i)+\sum_{j=1}^{i-1}(\frac{1}{\alpha+i-1})\delta_{\theta_j} $$ 这种条件先验分布由基础测度$P_0$和前一个受试者参数值上的概率质量的混合物组成。

  1. 中餐厅过程

​ 中餐厅过程的隐喻通常被用来描述Polya urn方案。想象一家拥有无限多张桌子的餐厅。第一个顾客坐在一张带有$\theta_ 1$盘子的桌旁。第二个顾客有$\frac{\alpha}{1+\alpha}$的可能坐在第一张桌子,有$\frac{1}{1+\alpha}$的可能坐在新桌子旁。这个过程继续下去,第i个顾客坐在一个被占用的桌子旁的概率与先前坐这桌的顾客数成正比,坐一张新桌子的概率与$\alpha$成正比。在这家有无限张桌子的餐厅,每张被占的桌子表示不同的主题集,渐近极限下以与$αlogn$成比例的速率添加新的集群。聚类的数量取决于(概率上)取决于受试者的数量n,当额外的被试者被添加到样本中时,根据需要引入新的聚类。在典型的应用中,这比k不依赖于n的有限混合模型更有意义,并且可以被认为是模拟拟合有限混合模型时的正式程序。当拟合一个有限混合模型时,在必要时手动添加新的混合成分来拟合数据。

​ 从Polya urn过程中条件分布的简单形式为后验计算和预测提供了一个有用的想法。从受试者i=1,…,n的可交换性,可以得到给定$\theta_{(-i)}=(\theta_j,j\neq i)$后$\theta_i$的先验分布: $$ \theta_i|\theta_{-i}\sim (\frac{\alpha}{\alpha +n-1})P_0(\theta_i)+\sum_{h=1}^{k^{(-i)}}(\frac{n_h^{(-i)}}{\alpha +n-1})\delta_{\theta_h^{(-i)}} $$ 其中,$\theta_h^{},h=1,...,k^{(-1)}$是$\theta^{(-i)}$的唯一值,并且$n_h^{(-i)}=\sum_{j\neq i}1_{\theta_j=\theta_h^*}$.

​ 用数据更新完整的条件先验(上式),得到具有相同形式的但在分量上具有更新的权重的$\theta_i$的条件后验分布,以及,只要$P_0$与核K共轭,更新$P_0$.例如,当K(·|θ)是正态核时,θ=(µ,$\phi$)为平均值和精度,$P_0$为正态伽马先验分布。在实现吉布斯采样时,人们可能可以从这些完整的条件后验分布中一次更新$\theta_i$中的一个。然而,这种方法的混合效果往往较差。

​ 另一种边缘吉布斯采样器,它分别更新集群的分配和集群特定的参数,步骤如下。设$\theta^$=$(θ^_1,…,θ_k^)$表示θ的唯一值,当$\theta_i=\theta_c^$,令$S_i=c$,因此$S_i$表示受试者i对簇的分配。吉布斯采样器交替在以下两种情况:

  • 1通过从多项式条件后验中采样来更新分配S:

$$ Pr(S_i=c|-)\propto \begin{cases}n_h^{(-i)}\kappa(y_i|\theta_c^*),,,c=1,...,k^{(-1)}\\ \alpha \int\kappa(y_i|\theta)dP_0(\theta),,,c=k^{(-1)}+1 \end{cases} $$

如果$S_i=k^{(-1)}+1$,那么受试者i被分配给一个单例集群。

  • 2通过以下采样来更新唯一值$\theta^*$:

$$ p(\theta_c^_|-)\propto P_0(\theta_c^_)\prod_{i:S_i=c}\kappa(y_i|\theta_c^*) $$

这只是参数模型下的后验分布,该模型将先验分布$P_0$分配给参数$\theta_h^$,然后用集群h中那些受试者的可能性更新这个先验。 当$P_0$与核$\kappa$共轭时,可以计算步骤1中的积分,步骤2中的条件后验具有除更新参数外与$P_0$相同的参数形式。例如,当核是高斯核时,θ为均值和方差,$P_0$是一个共轭正态反伽马先验,$\theta_c^$的条件分布具有高斯过程章节描述的相同的正态反伽马形式。也有一些修改可以适应非共轭的情况。

​ 在上述吉布斯采样器的步骤1中,第一个受试者被分配到样本中其他受试者之一所占据的现有集群,或者受试者被分配到一个新的集群。分配到一个新聚类的条件后验概率与DP精度参数α乘以第一个受试者数据的边际似然成正比,这是通过将似然K(yi|θ)积分到先验的θ∼$P_0$得到的。如果α接近于零或者这个边际可能性相对于被分配到一个被占用集群的数据的可能性很小,那么被试者i将倾向于被分配到一个现有的集群。因此,α和$P_0$在控制聚类的后验分布方面都发挥着重要作用。随着α的减小,聚类受试者的趋势不断增加,在极限下获得的公共θ的参数模型$y_i∼K(θ)$为α→0。在实践中,通常是设置α=1以支持分配到少数集群,或者为α选择一个伽马超优先级,以允许更大的数据自适应,并包括一个额外的MCMC步骤来更新α。

​ 在控制聚类行为中$P_0$的作用更微妙,而且经常被忽视。人们可能会天真地尝试一个高方差$P_0$来表达对不同核的可能位置的先验分布的无知。然而,正如第7.4节中讨论的那样,$P_0$的平坦先验可以做出强有力的假设,在这种情况下,有效地对引入新的集群造成了沉重的惩罚。这是因为当$P_0$的方差变得高时,边际似然会减小。在这种情况下,先验的$P_0$使小概率位于θ值的置信区间。在作为$P_0\rightarrow \infty$的方差的极限下,后验将表现近似于似然为$y_i∼K(θ)$,其中$\theta$是针对所有个体公共的。

使用Gibbs抽样求解狄利克雷过程混合模型

​ 假定有服从狄利克雷过程混合模型的观测数据集$X={x_1,...,x_n}$,由于观测数据条件独立,在对观测数据进行聚类分析时,不考虑数据的出现顺序。令$z_i$为每个数据的指示因子,$z_{'i}$表示从$Z={z_1,...,z_n}$中移出$z_i$后剩余数据组成的数据集。则关于$z_i$的条件分布: $$ p(z_i|x_1,...,x_N,Z_{'i},\lambda,\alpha_0)\propto p(z_i|Z_{'i},\alpha_0)p(x_i|z_1,...,z_N,X_{'i},\lambda) $$ 式中,等号右边第1项可以用Dirichlet过程中的CRP表示,由于各个观测量之间是可交换的,可以把第i个观测数据视为最后一个观测量。设$Z_{'i}$有K个类别,每一类中观测数据的个数为$N_k^{'i}$,第一项为: $$ z_i|Z_{'i},\alpha_0 \sim \sum_k^K \frac{N_k^{'i}}{N-1+\alpha_0}\delta(z_i,\bar k)+\frac{\alpha_0}{N-1+\alpha_0}\delta(z_i,\bar k) $$ 若第i个观测数据的指示因子为$z_i=k$,则有 $$ p(x_i|z_i=k,\chi_{'i},\lambda)=p(x_i|{x_j|z_j=k,j\neq i},\lambda)\=\frac{\int_{\theta}d(x_i|\theta)\prod_{z_j=k,j\neq i}f(x_j|\theta)g(\theta|\gamma)d\theta}{\int_{\theta}\prod_{z_j=k,j\neq i}f(x_j|\theta)g(\theta|\lambda)d\theta} $$

若$z_i=\bar k$为一新类别,则有 $$ p(x_i|z_i=k,\chi_{'i},\lambda)=p(x_i|\lambda)\ =\int_{\Theta}p(x_i|\theta)g(\theta|\lambda)d\theta $$ 因此, $$ p(z_i|x_1,...,x_N,Z_{'i},\lambda,\alpha_0)\ \propto \sum_k^K \frac{N_k^{'i}}{N-1+\alpha_0}×p(x_i|{x_j|z_j=k,j\neq i},\lambda)\delta(z_i,k)+\frac{\alpha_0}{N-1+\alpha_0}\int_{\Theta}p(x_i|\theta)g(\theta|\lambda)d\theta\delta(z_i,\bar k) $$

分层分布

​ 在前面的部分中,我们讨论了单个未知数的狄利克雷过程混合模型。为了构建丰富的半参数分层模型,可能需要合并几个DP过程,设置成独立的先验分布。这个方法在先前的 药理学数据分析中得到了应用。然而,这个方法有明显的限制,并且在很多情况下这个方法倾向于使用依赖于未知分布的先验。为了这种泛化的需要,我们要描述一个需要这种灵活性的应用。

超密度估计

非参数残差分布

在这里,密度估计主要是为了简化一个困难主题的阐述。狄利克雷过程混合模型(DPM)的真正吸引力在于,它们可以更广泛地用于放宽层次模型中的参数假设。这部分旨在介绍一些可能性,并不全面。首先,考虑具有非参数误差分布的线性回归: $$ y_i=X_i\beta+\epsilon_i,,,\epsilon_i\sim f $$ 其中$X_i=(X_{i1},……,X_{ip})$是一个预测因子的向量,$\epsilon_i$是一个分布为f的误差项。如前所述,均值中的线性假设很容易放宽。在这里,我们考虑放宽误差分布f具有参数形式的假设的问题。

在第17章中,我们认为t模型是降低异常值影响的一种方法。通过令$\epsilon_i∼N(0,φ^{−1}i1σ^2)$和$\phi_i∼Gamma∼(\frac{ν}{2},\frac{ν}{2})$,将分布$t_v$表示为法线的尺度混合。虽然由于t分布的尾巴沉重,t分布可能超过正态,但它仍然有限制性的形状,我们可以使用法态的DP尺度混合物进行非参数建模: $$ \epsilon_i∼N(0,φ^{−1}i1σ^2),,,,\phi_i∼P,,,,P\sim DP(\alpha P_0) $$ 其中,选定$P_0$对应$Gamma∼(\frac{ν}{2},\frac{ν}{2})$将f的先验集中在t分布上。得到的先验是灵活的,但仍然限制为单模和对称约为零。消除单模态和对称约束的另一种方法是使用高斯位置混合的f。从$X_iβ$项中去除截距,允许f具有未知的平均值,令: $$ \epsilon_i∼N(\mu_i,\tau^{-1}),,,\mu_i\sim P,,,P\sim DP(\alpha P_0),,,\tau \sim Gamma(a{\tau},b{\tau}) $$ 其中,选定$P_0$为$N(0,τ^{−1})$。密度估计的计算可以很容易地适应于包括更新回归系数β的步骤,然后在前面的步骤中用$y_i−X_iβ$替换$y_i$.

按组变化的参数的非参数模型

​ 在第15章中,我们考虑了具有不同系数的层次线性模型。系数分布的不确定性可以通过考虑DP或DPM的分布来考量。作为一个简单的例子,考虑单因子方差模型: $$ y_{ij}=\mu_i+\epsilon_{ij},,,\mu_i\sim f,,,\epsilon_{ij}\sim g $$ 其中,$y_i$是第一项的重复测量向量,$\mu_i$是受试者特异性的平均值,$\epsilon_{ij}$是观察特异性的残差。典型的参数模型在$g≡N(0,σ^2)$时,让f对应一个$N(\mu,ψ^{−1})$密度。为了更灵活地描述受试者之间的可变性,我们可以让: $$ \mu_i∼P,P∼DP(αP_0) $$ ,其中P是变化参数的未知分布,为了简单起见,我们将残差密度g建模为$N(0,σ^2)$。在分布上放置DP先验诱导一个潜在类模型,在这个模型中,被对象被分组为未知数量的簇: $$ \mu_i=\mu^∗_{S_i},Pr(S_i=h)=π_h,h=1,2,… $$ 其中$Si∈{1,…,∞}$是一个潜在类指数,$π_h$是分配到h的概率,这些概率遵循破棒形式。对于有限的潜在类模型,这个公式假设不同参数的分布是离散的,因此不同的受试者可以有相同的参数值。这可能是一个简化的假设,并且每个受试者的后验均值将是不同的,因为聚类是软的、概率的,$\mu_i$的后验均值在聚类分配的后验分布上得到平均值。

​ 在考虑非参数层次模型时出现了一些实际问题。首先,数据是否包含信息,允许非参数估计P,因为模型参数$\mu_i$没有直接观察到任何受试者。这个问题的答案和对结果估计的解释取决于每个受试者的观察次数。首先假设所有主题的$n_i=1$。在这种情况下,我们在数据中没有任何信息来区分受试者之间的可变性和受试者内部测量之间的可变性。然而,在残差密度g的正态性假设下,我们仍然有大量关于P的信息,因为P适应缺乏正态残差分布的拟合。在假设n≥1和正态g的一般情况下,P在允许受试者之间残差和系统变异性的正态分布缺乏拟合方面具有双重作用。当每个主题有许多观察时,后来的角色占主导地位,但当$n_i$对P的解释很少时,需要考虑到双重角色。

​ 消除这种混淆的一种自然可能性是也使用高斯分布的狄利克雷过程混合模型来建模g。在这种情况下,数据包含的分布信息较少,准确估计可能需要包含每个项目和许多项目的数据集。在参数分布和残差分布都未知的情况下,确实存在可识别性问题,即在非参数贝叶斯模型中,很难将分布的均值限制为零。然而,我们可以对过度参数化的模型进行MCMC分析,而不限制平均值,然后通过后处理来估计总体平均值和平均值为中心的参数和残余密度。

分层依赖

依赖的狄利克雷过程

​ 狄利克雷过程为单一随机概率测度P∼DP(αP0)提供了一个先验。在专注于彗星分析的应用和第j组内DNA损伤密度的建模的案例中,一种自然的方法是使用DP位置尺度的高斯混合物, $$ f_j(y)=\int N(y|\mu,\phi^{-1})dP_j(\mu,\phi),P_j\sim DP(\alpha_jP_0),\ =\sum_{h=1}^{\infty}\pi_{jh}N(y|\mu_{jh}^,\phi_{jh}^{-1}),\pi_j \sim stick(\alpha_j), (\mu_{jh}^,\phi_{jh}^)\sim P_{0j} $$ ​ 另一种更一般的策略是定义一个随机概率度量集合${P_1,...,P_T}$的依赖狄利克雷过程(DDP)先验。对于集合中每个概率测度的边缘先验都是狄利克雷过程的随机概率度量集合来说,一个依赖狄利克雷过程是很广泛的一类先验。

​ 定义DDP的一种非常方便的做法是依据一个破棒过程的构造: $$ P_j=\sum \pi_{jh}\delta_{jh}^, \pi={\pi_{jh}}\sim Q,\theta^{jh}\sim P_0 $$ ​ 其中,Q和$P_0$被选定后,于是就有$\pi_j=(\pi{j1},\pi_{j2},...)\sim stick(\alpha_j)$并且$\theta_{jh}^*$与$P_{0j}$在j的边缘分布上独立同分布,$j=1,...,T$。

​ 定义合适的Q使得依赖狄利克雷过程得到正确的边缘分布很困难,因此大多数文献关注所谓的“固定$\pi$的DDP(fixed-$\pi$ DDPs)”,表示为: $$ P_j=\sum_{h=1}^{\infty}\pi_h \delta_{\theta_{jh}^},\pi \sim stick(\alpha),\theta_h^\sim P_0 $$ 于是在不同分量上的概率权重是相同的,并且只有原子位置改变。原子位置由P0的值控制。

​ 回到先前的应用,我们使用固定$\pi$的DDP混合高斯核: $$ f_j(y)=\sum_{h=1}^{\infty}\pi_h N(y|\mu_{jh}^*,\phi_h^{-1}),\pi\sim stick(\alpha),\phi_h\sim Gamma(a,b) $$ 于是剂量组的权重和带宽是相同的,但核的位置不同。

​ 尽管DDP混合模型可能看上去复杂,但它实际上是DPM模型的简单转变,因为一个单一密度和先验计算就能得到相似曲线。如果用吉布斯抽样,步骤与先前章节介绍的类似。令$y_i \sim N(\mu_i,\phi_i^{-1})$,有以下步骤:

  1. 更新混合分量索引$S_i$。通过从以下多元条件先验中抽样:

$$ Pr(S_i=h|-)=\frac{\pi_hN(y_i|\mu_{x_ih}^_,\phi_h^{_-1})}{\sum_{l=1}^N\pi_lN(y_i|\mu_{x_il}^_,\phi_l^{_-1})}, h=1,...,N $$

  1. 像先前章节那样通过beta全条件更新破棒过程的权重。
  2. 从伽马全条件后验分布更新核特定准确度$\phi_h^*$:

$$ (\phi_h^_|-)\propto Gamma(\phi_h^_|a,b) \prod_{i:S_i=h}N(y_i,\mu_i,\phi_h^{*-1}) $$

  1. 从$\mu_{1h}^$的高斯全条件更新$\mu_{1h}^$,从$\beta_{jh}$的全条件更新$\beta_{jh}$。

分层狄利克雷过程

​ 固定$\pi$的DDP在表征未知分布的层次依赖结构方面的灵活性有限。一个比较常用的DDP的替代选择是有层次的DP(HDP),通过令: $$ P_j\sim DP(\alpha P_0),P_0 \sim DP(\beta P_{00}) $$

对应于在未知基度P0条件下为每个$P_j$选择独立的DP先验,同时也分配了一个狄利克雷过程。作为缩写形式,表示为$P\sim HDP(\alpha,\beta,P_{00})$。在破棒过程中,有: $$ P_j=\sum_{h=1}^{\infty}\pi_{jh}\delta_{\theta_h^},\theta_h^ \sim P_{00} $$ 因此,由于DP实现的离散性,作为常见的向基测度$P_0$分配一个DP先验的自然结果,我们在所有组特定的随机概率度量中使用相同的因子,同时允许整个组的权重出现偏差。

​ 这就使得从固定$\pi$的DDP中得出一个不同的结构,当使用相同权重时,这个结构让因子也能变化。

​ 为了说明这个结构,假定: $$ f_j(y)=\int N(y|\mu,\phi^{-1})dP_j(\mu,\phi),P\sim HDP(\alpha,\beta,P_{00}) $$ 在这种情况下,该模型引入了一个具有不同位置和尺度的普通核的公共词典:$N(\mu_h^*,\phi_h^{-1}),h=1,2,...,\infty$.这里有一个中心密度$f_0(y)$,其特征是在这个词典与权重λ混合,特定组的密度使用相同的字典表示,但权重$\pi_j$是从以λ为中心的破棒过程中提取的。超参数$\alpha$控制权重中组间的可变性,α→0表示群特定密度是高斯的,在高斯核中群簇具有相同的平均值和精度。在另一个极端α→$\infty$时,用fj≡f0和f0建模作为一个DPM位置尺度高斯混合,来获得组之间的池。

​ 分层狄利克雷过程的另一个含义是分层聚类。为了阐述这一点,我们关注一个应用程序,其中有美国各州的索引i=1,...,n,在州i中各个医院的索引$j=1,...,n_i$。$y_{ij}$表示在州i的医院j的医疗服务质量。假设$y_{ij}\sim f_i$并且分配一个高斯先验的HDP位置混合给每个州的医疗服务质量概率密度的集合${f_i}$。我们可以得到以下分层模型: $$ y_{ij}\sim N(\mu_{S_{ij}}^,\phi_{S_{ij}}^{-1})\ Pr(S_{ij}=h)=\pi_{jh},(\mu_{S_{ij}}^,\phi_{S_{ij}}^{-1})\sim P_{00} $$ 其中$S_{ij}=h$表示在第i州的医院j的医疗服务质量集h。由于权重上的分层结构,在同一州的医院更有可能被分配进同一个集合,但是人们还可以获得跨州的医院集。这种软概率聚类可能在某些应用中很有用,同时描述了HDP如何跨组借用信息。由于超参数$\alpha$和$\beta$分别在控制聚类中的组内依赖性和聚类总数方面起着关键作用,因此,允许数据通过超先验来获取它们的值十分重要。一个常见的选择是$\alpha\sim Gamma(1,1)$,独立于$\beta\sim Gamma(1,1)$。

嵌套狄利克雷过程

​ HDP的工作原理是在不同的群体的特定分布中合并相同的原子,同时允许权重的变化,这导致跨组的依赖聚类。在许多应用中,使用具有相同分布的集群组更好。例如,在彗星实验的应用案例中,我们可以聚集在DNA损伤分布上没有差异的剂量组,而在医院案例应用中,我们可以聚焦于本州内医院的护理质量分布上没有差异的州。在前一个案例中,我们获得了多重治疗组比较的非参数贝叶斯方法,每一对剂量组和所有剂量组的后验概率都相同。这些概率可以形成关于治疗组平等性假设的贝叶斯检验的基础。

​ 为了完成这种分布聚类,我们可以依赖一个嵌套的狄利克雷过程(NDP)混合模型。这个NDP可以表示如下: $$ P_j\sim P,P\sim DP(\alpha P_0), P_0\equiv DP(\beta P_{00}) $$

这和HDP看上去很类似,它从具有DP先验的DP中提取特定组的随机概率度量。然而在这里,我们从一个共同的随机概率测度中提取$P_j$,该量来自一个DP,并且基是DP,而不是DP的实现。当我们检查该形式的破棒表示时,区别就清晰了起来: $$ P_j\sim P=\sum_{h=1}^{\infty}\pi_h\delta_{P_h^},\pi \sim stick(\alpha),P_h^\sim DP(\beta P_{00}) $$ 因此,P采取了寻常的DP破棒形式,但原子对应于从DP中提取的独立同分布的随机度量。这导致了随机概率测度的聚类,其中$P_i=P_{j'}$的先验概率是$\frac{1}{1+\alpha}$,是DP Polya urn方案的自动结果。然而,如果$P_j$和$P_{j'}$在不同的团簇中,它们将有独立于$P_{00}$的不同原子。与HDP不同,HDP使用一组共同的原子来形成$P_j$,因此$P_r(P_j=P_{j'})=0$.

在实践中,$P_j$通常被用作NDP混合模型中的混合分布,用于群体特定密度的集合,例如: $$ f_j(y)=\int N(y|\mu,\phi^{-1})dP_j(\mu,\phi), P\sim NDP(\alpha,\beta,P_{00}) $$ 其中,$P\sim NDP(\alpha,\beta,P_{00})$用作先验的简写形式。在这种情况下,通过$P_r(f_i=f_{j'})=\frac{1}{1+\alpha}$,组特定密度将被分配给簇。我们可以选择一个超先验分布$\alpha\sim Gamma(a,b)$,其中a,b是为了得到$P_r(H_{0jj'})$和$P_r(H_0)$的期望值引入的。这使得数据可以被告知给$\alpha$,而且随着T的增加,数据会被感知的更强烈,从而导致所谓的“维度祝福”。

​ 在机器学习文献中成功实现的一个自然的NDP-HDP修改是在$P_{00}$处放置一个DP,以便簇特定的密度共享一组共同的全局原子,但有着不同的权重。这结合了NDP和HDP的先验,潜在地利用了这两种方法的优势。

Convex混合模型

​ HDP和NDP都是DDP框架的特殊情况,因为它们在随机概率测度集中包含了依赖性,同时保持单个RPM的DP先验。尽管DP具有一些吸引人的特性,并且在某种意义上是一个典型情况,但它依然限制了模型的灵活性,即总是限制于对DDP的关注,而不考虑更广泛的先验类,这些先验在不保持DP边缘的情况下以不同的方式合并依赖性。诱导依赖的另一种方法是使用分量随机概率测度的随即凸组合。例如,假设兴趣的焦点集中在结合不同研究中心按密切相关的方案进行纵向研究的数据。特别地,令$y_{cij}$表示中心c对第i个个体的响应,令$x_{cij}$表示相应的协变量。然后,我们可以考虑如下的分层模型: $$ y_{cij}\sim N(x_{cij}\beta_{ci}+\epsilon_{cij},\sigma^2) $$

其中$\beta_{ci}$是中心c和个体i的特定系数的p×1维向量,$\epsilon$是一个残差,$\sigma^2$是残差的方差。问题于是变为如何从不同中心读取信息。如果假设一个参数模型,那么使用多层结构将$\beta_{ci}$分解为$\beta_{ci}=\alpha_c+\Psi_{ci}$.中心特异性效应$\alpha_c$建模为以α为中心的多变量高斯,个体特异性偏差$\Psi_{ci}$建模为以零为中心的多元高斯。

​ 用一种更复杂的半参数方法,我们可以令: $$ \beta_{ci}\sim P_c,{P_c}\sim \Pi $$ 其中,$P_c$是中心受试者间变异性的分布,而$\Pi$是不同分布的联合先验。HDP或NDP可能可以用于$\Pi$,但一个简单的选择是令: $$ P_c=\pi G_0+(1-\pi)G_c, G_c\sim DP(\alpha G_0),\pi\sim Beta(a,b) ,,,,, 26.3 $$ 于是c组中的分布表示为全局分布G0和组特定分布$G_c$的混合物。$G_0$上的概率权值$\pi$越高,各组间的分布就越相似。虽然有相似的动机,但这个先验不同于HDP,包括一组单独的全局和群特异性原子。

​ 分配到$G_0$内全局原子的受试者可以与其他组的受试者聚在一起,而分配到中心特定原子的受试者只能与同一中心的其他受试者聚在一起。后验计算可以简单地使用数据增强方法,该方法引入一个二进制指示器$z_{ci}=1$,表示对全局组件的分配。从这些指标的完整条件计算后,就有条件独立的DP先验,可以直接使用狄利克雷过程中的计算算法。边缘上,$P_c$没有DP先验,因此$\Pi$不是DDP。

​ 表达式(26.3)让人想起一个具有总体均值和中心偏差的高斯层次模型。然而,由于我们建模的是随机概率度量。而不是实值随机变量,所以用凸组合而不是加性结构。通过使用分量随机概率度量的凸组合,我们确保了结果是一个概率度量。相关的凸组合可以广泛地用于在度量集合中包含更多结构化的依赖性。例如,回顾先前讨论的彗星实验应用。在这种情况下,剂量组有一个自然的顺序,选择一个比$P_t$和$P_{t+1}$更依赖的先验是有意义的,可以通过定义一个一阶自回归模型来实现: $$ P_t=(1-\pi)P_{t-1}+\pi G_t,,,P_0=G_0,,,G_t\sim DP(\alpha P_0) $$ 因此,剂量组t的RPM等于前一个剂量组的随机概率测量和剂量组特异性偏差的混合物。这在动机上类似于高斯随机游走,但是是在度量空间中的随机游走,参数$\pi$控制了依赖水平。人们可能会让$\pi$依赖于t以增加额外的灵活性。这种DP模型的动态混合在时间序列应用中也很有用,但是尽管早期引入的原子不断被分配越来越少的权重,一个潜在缺陷是原子只能随着时间推移被引入,并且永远不会消失。一种规避这个问题的方法是从HDP中提取$P_0$,于是相同的原子才能随时间推移重复使用。这种结构已经成功地应用于音乐数据的分析。

密度回归

​ 前一节关注的是我们有一个有限的随机概率度量的集合,它们要么遵循简单的顺序,要么是可交换的。在许多应用中,设置并不那么简单,考虑不可数的集合更自然,对于: $$ P_{\chi}={P_x,x\in\chi},\chi⊂ℜ^p $$ 其中$x=(x_1,…,x_p)$是预测器的向量,$P_x$是特定于预测值x的随机概率度量,$\chi$表示预测器的域,$P_{\chi}$是在每个预测值上定义的随机概率度量的集合。其中一个动机来自于密度回归的设置。在先前章中,我们讨论了使用高斯过程的密度回归,接下来我们讨论一种基于狄利克雷过程的替代方法。

​ 理想情况下,我们允许给定预测因子的响应的整个条件密度p(y|x)随着x的变化而灵活地变化。一种已被广泛使用的方法是分层专家混合模型,即: $$ p(y|x)=\sum_{h=1}^Hπ_h(x)N(y|xβh,τ_h^{−1}) $$ 这对应于将条件密度表示为正态线性回归的混合,不同混合分量的权重随预测器灵活变化。在机器学习文献中,有限H的情况下,一种常见的方法是依赖于$π_h(x)$的概率树模型,尽管人们也可以使用一种更简单的方法,如逻辑回归。

依赖的破棒过程

作为一个依赖于混合物的非参数贝叶斯密度回归模型,我们可以令: $$ p(y|x)=\int N(y|x\beta,\tau^{-1})dP_x(\beta,\tau),,,P_{\chi}\sim \prod_{\chi} $$ 这是一种线性回归的混合模型,是依赖于预测器的线性回归,但是以更一般的形式。这里,$Π_{\chi}$表示不可计的混合测量集合${P_x,x∈\chi}$的先验。当参数模型近似正确时,将先验集中在一个合理的参数模型上,以有利于后验的扁平化是有用的。

​ 主要的问题是如何选择$Π_\chi$。在这方面,很自然地想到依赖预测的破棒过程: $$ P_x=\sum_{h=1}^{\infty}\pi_h(x)\delta_{\theta_{h^*}}(x),,,,\pi_h(x)=V_h(x)\prod{l<h}(1-V_l(x)) $$ 其中$π_h(x)$是预测值x特定的分量h的权重,$V_h(x)$是步骤h中断的概率棒的比例,$\theta_{h^*}(x)$是一个预测相关的原子。在这里,我们将为了简单起见关注$\theta{h^*}(x)=\theta{h^*}=(\beta{h^*},\tau{h^*})$的情况,从而有一个回归系数向量和精度的单一全局集合。然后我们有 $$ P(y|x)=\sum{h=1}^{\infty}\pi_h(x)N(y|x\beta_h^*,\tau_h^{-1}),,,\pi_h(x)=V_h(x)\prod_{l<h}(1-V_l(x)) $$ 我们希望能够对我们的大部分数据使用少量的专家,这可以通过选择一个$π_h(x)$在指数h中快速下降的先验来得到支持。如果$V_h∼Q$是由一个随机过程生成的,而对于所有x,$V_h(x)∼Beta(1,α)$都处于边缘,那么我们将得到一个线性回归的DDP混合物。Q提出了不同的选择,包括基于命令的DDP和本地DP。

​ 然而,在计算上有一些明显的优点,即不限制对DDP的注意。其中一个先验在各种应用中具有良好的实际性能,并已被证明在估计条件密度时具有很大的支持和后验一致性,那是核破棒过程, $$ V_h(x)=K_{\Psi_h}(x,\Gamma_h)V_h,,,,\psi_h\sim H,,,\Gamma_h\sim G,,,V_h\sim Beta(1,\alpha) $$ 其中$K_ψ(·,Γ)$是一个内核,它位于Γ,带宽为ψ。选择内核$K_{ψ_h}(x,Γ_h)$来获得$x=Γ_h$的最大值为1,在这种情况下,$V_h(x)=V_h$。当x远离$Γ_h$时,内核减少,导致相应的$V_h(x)$减少。我们可以把这个过程看作是生成一个随机位置$Γ_h$,具有相应的粘断随机变量Vh和原子$(β_h^*,τ_h^∗)$。由于内核的合并,当x位于靠近$Γ_h$时,$π_h(x)$将会趋向于更大,特别是当指数h较小时。在预测空间中权值的空间可变性是由核带宽控制的,通过核特定的带宽,我们允许在某些区域发生更快速的变化。当内核变得平坦,对于所有x和h的$K_{ψ_h}(x,Γ_h)≈1$,我们得到了一个DP破棒过程作为极限情况。

​ 核破棒过程混合物的后验计算往往是直接、有效的,例如,通过预先为带宽和位置指定一个潜在值的网格,以促进吉布斯采样。然而,另一种概率破棒规范有一些计算优势,令:

$$ \pi_h(x)=V_h(x)\prod_{l<h}(1-V_l(x)),,,,V_h(x)=\Phi(\alpha_h+\mu_h(x)),,,,\alpha_h\sim N(\mu,1) $$ 其中Φ(z)为标准法态cdf,$\mu_h:\chi→ℜ$为任意回归模型。为了激发单位概率破棒过程(PSBP),首先考虑没有预测因子的基线情况,以便在$\alpha_h\sim N(\mu,1)$时,$V_h(x)=V_h=Φ(α_h)$。这个模型类似于狄利克雷过程,但我们不是从独立同分布的$Beta(1,α)$生成$V_hs$,而是我们通过转换独立同分布的$N(µ,1)$通过标准的正态cdf绘制到单元区间(也可以使用逻辑转换并获得逻辑破棒过程)。在µ=0的特殊情况下,我们得到了$V_h∼Beta(1,1)$,因此µ=0和DP精度为1的PSBP是等价的。一般来说,µ超参数在控制DP中破棒随机变量的下降率和样本中簇数量的相关先验方面起着与精度α类似的作用。对于大的µ,我们得到了$V_h≈1$,因此在第一个分量上的大权重类似于DP中的$α≈0$。

狄利克雷过程模型的主要应用

机器学习中的应用

​ 狄利克雷过程经常用于贝叶斯非参数统计,这里的非参数并不代表无参数模型,而是表示随着观察得到的数据更多而增长的模型。由于这样的灵活性,在机器学习尤其是无监督学习中,贝叶斯非参数模型取得了相当程度的普及。

​ 在贝叶斯非参数模型中,先验分布和后验分布不是参数分布,而是随机过程。Dirichlet分布是概率总和为1的一组非负数的单纯形上的概率分布,这一事实使其成为对分布上的分布或者函数上的分布进行建模的良好选择。

​ 此外,在聚类问题中,如果不同数量的聚类事先未知,狄利克雷模型的非参数性质使其这类问题的良好选择。

​ 另外,狄利克雷过程还用于在监督学习背景下开发混合专家模型。例如,高斯过程专家的混合,其所需要的专家的数量必须从数据中推断出来。

作为先验

​ 由于狄利克雷过程的提取是离散的,因此还有个重要用途是作为无限混合模型中的先验概率。在这种情况下,${B_k}$是分量分布的参数集。因此,生成的过程是从狄利克雷过程中抽取一个样本,对于每一个数据点,依次从该样本分布中抽取一个值并用作该数据点的分量分布。在混合分量的数量没有预先明确定义时,由于狄利克雷过程可以生成的不同分量的数量没有限制,因此可以适用。例如,高斯模型的无限混合以及相关的混合回归模型。

自然语言处理

​ 狄利克雷过程模型的无限性质也使它们适用于自然语言处理应用程序。在这些应用程序中,通常需要将词汇表视为无限的离散集合。

非参数假设检验

​ 狄利克雷过程也可以用于非参数假设检验,即开发经典非参数假设检验的贝叶斯非参数版本,如:符号检验、Wilcoxon秩和检验、Wilcoxon符号秩检验等。

​ Wilcoxon秩和检验与Wilcoxon符号秩检验的贝叶斯非参数版本是通过使用不精确的狄利克雷过程开发的,这是一个不知道先验的狄利克雷过程。

几个应用案例

基因毒性应用

​ 假设我们有一个实验研究,在一个连续反应变量上对不同组的“受试者”进行观察。令$y_i$表示对主体i的响应,$x_i\in {1,...,T}$表示组。例如,在基因研究中,使用单细胞凝胶电泳(也被称为“彗星实验”)来测量DNA损失,$y_i$表示在细胞i中DNA的损伤量,$x_i$表示潜在的基因毒性暴露的剂量组。这类研究的重点是评估细胞间DNA损失的密度如何随剂量而变化。图26.2显示了暴露于不同水平的过氧化氢、已知的基因毒性化学物质的剂量组的直方图和核平滑密度估计,并且使用彗星实验在单个细胞水平上测试DNA损伤。从图中可以看出,响应密度较低的分位数并不随暴露水平变化而变化很大,而较高的分位数则随着剂量的增加而明显增加。这与预期的,由于细胞之间的可变性,实验中不可能对所有细胞获得相同的剂量。

image-20211019101927920

 图26.2 在基因毒性案例中,每个过氧化氢剂量组细胞间DNA损伤的直方图和核平滑密度估计

​ 继续以上实验,我们将大量细胞暴露于0、5、20、50或100µmol过氧化氢中,允许0、60或90分钟后的修复时间后测量单个细胞的DNA损伤。测量细胞时,使用索引$i=1,...,n$,测量的细胞i的响应$y_i$是橄榄尾矩,它是使用彗星实验获得的DNA链断裂频率的替代品。

​ 该研究的目的是评估彗星实验对检测已知基因毒性药物-过氧化氢引起的损伤的敏感性,同时也研究损伤修复的速度。设xi∈{1,……,K}是一个组指数,表示细胞i的过氧化氢水平和修复时间。每个剂量×修复时间值的xi值,以及组间已知的随机排序限制如图26.2所示。随着剂量组内修复时间的增加,DNA损伤密度会随机降低,而修复时间为零、剂量组增加,DNA损伤密度会增加。样本量为1400个,除第9组和第13组有50个细胞外,每组有100个细胞。

 图26.3 有向图说明了遗传毒性模型中的顺序限制。箭头指向随机上更大的组。$H_{1k}$的后验概率也展示了出来。

​ 我们希望评估DNA损伤是否随着过氧化氢剂量的增加而增加,以及损伤是否随着修复时间的增加显著减少。我们使用一个受限制的DDP来建模这些数据。每一组内的DNA损伤密度被描述为高斯分布的狄利克雷过程位置混合物。我们根据相邻的组差异进行参数化,并使用类似的先验,但稍作修改以限制集群特定的平均差异,以服从图26.3中的顺序。图26.3中的第k条有向边连接了两个未知密度,其特征是高斯分布的DPMs,具有相同的权值和核带宽,但可能存在不同的位置。设$d_k$表示由第k条边连接的组之间不同的簇(混合成分)的总概率权重。如果$d_k$很小,这意味着这两个密度是相似的,提供了一个简单的标量摘要。

​ 图26.3显示了每条边对应的组间比较的局部排序的后验概率。对操作特性的模拟研究表明,该测试程序在低I型错误率和高功率方面,即使在组间有细微差异的小样本中,也具有良好的频率性能。我们的结果提供了强有力的证据,表明在修复时间为0min,过氧化氢浓度在0、5和20µmol的剂量组之间,DNA损伤随过氧化氢浓度增加而增加,而在较高剂量水平下进一步增加的证据不太清楚。正如预期的那样,没有证据表明第1、6和11组之间的分布发生了变化,因为没有诱导损伤需要修复。然而,在修复时间为60min后,每个暴露组的DNA损伤都明显减少。增加的额外的30min的修复并没有显著改变分布。

​ 这些结果与对原始数据的主观检查一致,图26.3所示的贝叶斯非参数密度估计与各组分别获得的简单核平滑密度估计一致。如果图中数据表明相邻组有相似或相同的分布,这些组就会被有效地合并,获得较高的后验概率$H_{0k}$和接近相同的密度。这种方式显著减少了估计个体密度的均方误差,同时通过对组的比较进行推断。

糖耐受性预测

​ 我们将概率破棒过程应用于糖尿病的流行病学研究。重点是评估n=868例患者中$y_i$=糖耐量(GT)、$x_{i1}$=对数转换胰岛素敏感性(IS)和其他糖尿病危险因素$x_{i2}$=年龄、$x_{i3}$=腰臀比(WTH)、$x_{i4}$=体重指数(BMI)、$x_{i5}$=舒张压(DBP)和$x_{i6}$=收缩压(SBP)之间的关系。GT通过口服糖耐量试验中的2小时血糖水平(mg/dl)来测量,并表明血糖从血液中清除的速度。GT也被用于诊断2型糖尿病,GT<140(正常)、GT$\in$[140,200](糖尿病前期),GT>200(糖尿病)。IS提供了身体对胰岛素反应如何的指标,胰岛素是一种调节葡萄糖从血液细胞到体细胞运动的激素。

​ 图26.4绘制了2小时血糖水平与IS、年龄、腰臀比(WTH)、体重指数(BMI)、舒张压(DBP)和收缩压(SBP)的关系。葡萄糖分布有较大的右偏,分布形状随IS而变化。由于这些数据不支持线性或非线性平均或中位数回归模型,我们应用贝叶斯密度回归方法来允许2小时葡萄糖的分布灵活地改变不同的风险因素,同时也允许风险因素退出模型和影响局部的特定区域的预测空间。

​ 各预测因子的边际后验概率分别为1.0、1.0、0.03、0.02、0.03、0.03,说明IS和年龄是葡萄糖分布变化的重要因素,其他预测因子可以丢弃。图26.5显示了估计的条件密度p(y|x),IS和年龄在其第5、50和95个经验百分位数中发生变化。葡萄糖密度有一个沉重的右尾。但是,随着IS的增加,右尾就消失了。右尾是2小时血糖水平超过200mg/dl的人群(参考线为0.2,标准化)。随着年龄的增长,右尾变得更重,特别是对于那些IS较低的受试者,这意味着衰老也与GT较差有关。

图26.4 糖耐量研究数据:y=2小时血糖水平(mg/dl)x1=胰岛素敏感性x2=年龄x3=腰臀比x4=体重指数x5=舒张压;x6=收缩压。 

image-20211019095308938

图26.5 预测(虚线)条件反应密度p(y|x)和95%可信区间(虚线),归一化x1(胰岛素敏感性)和x2(年龄)在第5、50、95个经验百分位数之间发生变化。 

参考文献

  1. Anrew Geiman(2013),Bayesian Data Analysis,CRC Press