Dirichlet分布的推导与理解

2023-10-26

1.概述

Dirichlet 分布与贝塔分布、伽马分布有着紧密的联系，在贝叶斯统计中经常被用作其它概率分布如多项分布的先验分布，且在LDA分析中得到了广泛应用，本文结合直观理解以及详细的数学推导得到狄利克雷分布具体形式，并结合可视化以加深理解。

2.直观理解

quaro问答网站上有一个举例非常直观，多项分布源于一个投掷骰子的过程，Dirichlet 分布可以视为是一个骰子工厂生产骰子的过程。也就是生产骰子是投掷骰子的前置步骤，对应上述的Dirichlet 分布(生产骰子）就是多项分布(投掷骰子）的先验分布。

比如投掷骰子N次，统计其中一个面朝上的次数具有相应的多项分布。如果Dirichlet的参数较大，则生产骰子对应的模具具有较高的精度，它生产的所有骰子都将具有接近平均多项式的系数，也就是各面无偏的骰子。如果模具精度较低，生产出的骰子将具有较高的方差(一个骰子可能有很高的概率落在1上，另一个可能有很高的概率落在5上，等等)。极端情况是当某些参数 $\small \theta<<1$ 时，方差非常大，pmf集中在角落，也就是骰子的投掷结果对某一面表现较高的倾向性和有偏性。

下面介绍如何使用 Dirichlet 分布来表征多项分布的随机可变性。还是以制造六面骰子为例，但允许抛掷的结果只能为1、2或3（这样可以简化后面的可视化，对于分析真实的六面骰子类同)。

如果骰子是完全无偏，那么三个结果的概率将相同且等于 1/3。我们可以将结果的概率表示为向量 $\small \boldsymbol \theta = (\theta_1,\theta_2,\theta_3)$ 。 $\small \boldsymbol \theta$ 有两个重要的属性：其一，每一项的概率之和必须等于 1,其二，任何概率都不能为负,也就是如下形式。

$\small \sum_{}^{}\theta_i=1$

$\small \small \theta_i\geq 0$

当这些条件成立时，我们可以使用多项分布来描述与掷骰子相关的结果。

如果我们投掷骰子 $\small N$ 次， $\small K$ 个面的投掷结果表示为向量形式 $\small \small D=\{x_1,...,x_K\}$ ，其中 $\small x_i$ 表示某个面出现的概率，则似然函数的形式为：

$\small p(\mathcal D|\theta)=\prod_{k=1}^{3}\theta_k^{N_k},N_k=\sum_{i=1}^{N}\mathbb{I}(y_i=k)$

其中， $\small N_k$ 是第 $\small k$ 面出现的次数， $\small k\in\{1, 2, 3\}$ 。我们期望生产的骰子的特性会有一些可变性，即便是我们尝试生产无偏的骰子，但由于生产过程中的各种因素的可变性，我们也不能完全期望特定骰子的每个面出现的概率正好就是 1/3，也就是完全无偏的骰子只存在理论上的可能性。

为了从数学上描述这种可变性，比如我们想知道在特定的制造过程中， $\small \boldsymbol \theta$ 的每个可能值的概率密度。为此，让我们将 $\small \boldsymbol \theta$ 的每个元素视为一个独立变量。也就是说，对于 $\small \boldsymbol \theta = (\theta_1,\theta_2,\theta_3)$ ，我们可以将 $\small \theta_1$ , $\small \theta_2$ , $\small \theta_3$ 视为一个个彼此独立的变量。

由于多项式分布要求这三个变量之和为 1，因此我们知道 $\small \small \boldsymbol \theta$ 的允许值被限制在一个平面内，也就是三维空间中三点(1,0,0)、(0,1,0)、(0,0,1)确定的一个平面。此外，由于每个值 $\small \theta_i$ 都必须大于或等于零，所有 $\small \boldsymbol \theta$ 允许值的向量集合被限制在一个三角形内(进一步限定了平面的范围局限于三角形内，超出三角形必然是至少一个值大于1，一个小于0)。我们想知道的是这个三角形上每个点的概率密度。这就是 Dirichlet 分布可以帮助我们的地方，我们可以将它用作多项分布的先验。

3.数学推导

在详细推导Dirichlet 分布之前，我们先介绍概率分布中常见的共轭先验，然后结合数学上几种常见概率分布及单纯形分布引出具体的Dirichlet 分布。

3.1 共轭先验

在贝叶斯概率论中，如果后验分布 $\small p(\boldsymbol \theta|\mathbf x)$ 和先验分布来自同一个概率分布，则先验和后验称为共轭分布，先验分布是似然函数的共轭先验。

如果我们考虑从给定数据集 $\mathbf x$ 推断分布参数 $\boldsymbol \theta$ 的问题，那么贝叶斯定理表明，后验分布等于似然函数 $\small p(\mathbf x|\boldsymbol \theta)$ 与先验函数的乘积，再由数据的概率密度进行归一化：

$\large p(\boldsymbol \theta|\mathbf x)=\frac{p(\mathbf x|\boldsymbol \theta)p(\boldsymbol \theta)}{\int p(\mathbf x|\boldsymbol \theta')p(\boldsymbol \theta')d\boldsymbol \theta'}$

由于似然函数通常是从数据生成过程中定义的，因此先验的不同选择会使积分的计算或多或少变得困难。如果先验与似然具有相同的代数形式，那么我们通常可以获得后验的闭式表达式，而无需进行数值积分。

3.2 相关概率分布

a) 二项分布

$\small Bin(k,n|p)=Pr(X=k))\triangleq \binom{n}{k}p^k (1-p)^{n-k}$

$\small \binom{n}{k} \triangleq \frac{n!}{(n-k)!k!}$

b) 多项分布

$\small Mu( \mathbf x|n,\mathbf p)=\begin{pmatrix} n\\ x_1,x2,...x_K \end{pmatrix}\prod_{j=1}^{K}p_j^{x_j}$

$\small \begin{pmatrix} n\\ x_1,x2,...x_K \end{pmatrix}=\frac{n!}{x_1!x_2!\begin{matrix}\cdot\cdot \cdot x_K! \end{matrix}}$

$\small x_1+x2+...+x_K =n$

c) 贝塔分布

$\small Beta(x|\alpha,\beta) =\frac{1}{B(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1}$

$\small B(\alpha,\beta) \triangleq \frac{\Gamma(\alpha) \Gamma(\beta) }{\Gamma(\alpha+\beta)}$

d)单纯形分布

单纯形分布一般指的是概率单纯形(probability simplex)，是一个数学空间概念，即每个点代表有限数量的互斥事件之间的概率分布。每个事件通常称为一个类别，我们通常使用变量 K 来表示类别的数量。

概率单纯形上的一个点可以用和为1 的 K 个非负数表示。这里给出一些具体的例子：

K=2 的1维单纯形中的一个点, (0.6, 0.4)

K=3 的2维单纯形中的一个点, (0.1, 0.1, 0.8)

K=6的5维单纯形中的一个点，(0.05, 0.2, 0.15, 0.1, 0.3, 0.2)

当K=2时，这个数学空间是 K-1维，是一条直线，当 K=3 时，它是二维单纯形，也就是一个三角形，当 K=4 时，它是三维单纯形，是一个四面体。在每种情况下，单纯形都是 (K-1) 维对象。（数字总和为 1 的求和项数量减少 1得到维度）。概率单纯形在贝叶斯推理中很常见。例如，假设我们在 A、B 和 C 之间做出决定，从而使我们的假设空间定义为{A, B, C}。我们对 A、B 或 C 为真的相对可能性的信念就落在 K=3 的概率单纯形上。

概率单纯形的每个“角”或“顶点”代表所有概率都放在一个类别上的情况。因此，例如在上述假设空间 {A, B, C} 的情况下，点 {A=0, B=1, C=0} 表示所有概率都放在 B 上的信念。如果单个假设被认定为不可能，其他 K-1 个假设的概率之和为 1，因此 K维单纯形的这个“边界”实际上是一个 K-1 单纯形。在博弈论中，代理人的 K 个不同潜在行动之间的“混合策略”存在于 K-单纯形上。

单纯形概率是定义在K-单纯形上的。也就是对于给定的参数K+1，如果满足如下条件：

$\small \small \theta_k\geq 0,k\in \{1,...,K,K+1\}$

$\small \small \sum_{k=1}^{K+1}\theta_k=1$

$\small \boldsymbol \theta \in \mathbb{R}^{K+1}$

则称 K+1维的向量 $\small \boldsymbol \theta=(\theta_1,..., \theta_K,\theta_{K+1})$ 定义一个K-单纯形或K维单纯形( unit K-simplex）。

单纯形：几何学上，单纯形或者n-单纯形是和三角形类似的n维几何体。精确的讲，单纯形是某个n维以上的欧几里得空间中的（n+1）个仿射无关（也就是没有m-1维平面包含m+1个点；这样的点集被称为处于一般位置）的点的集合的凸包。3维单纯形，也叫四面体。

例如，0-单纯形就是点，1-单纯形就是线段，2-单纯形就是三角形，3-单纯形就是四面体，而4-单纯形是一个五胞体（每种情况都包含内部）。正单纯形是同时也是正多胞形的单纯形。正n-单纯形可以从正(n − 1)-单纯形通过将一个新顶点用同样的边长连接到所有旧顶点构造。

注意：n-单纯形对应的顶点数量为n+1

标准n-单纯形的顶点为

e0 = (1, 0, 0, …, 0),

e1 = (0, 1, 0, …, 0),

...

en = (0, 0, 0, …, 1).

e）狄利克雷分布

具体来说，狄利克雷分布Dirichlet_distribution的定义形式不具有唯一性，下面给出一种定义形式，便于后面的推导。

令 $\small X$ 为一个 $K \times 1$ 的连续随机向量，它的定义域为

同时令 $\alpha_1,...,\alpha_{K+1} \in R_{++}$ . 我们就说 $\small X$ 具有以参数 $\alpha_1,...,\alpha_{K+1}$ 的 Dirichlet 分布，当且仅当它的联合PDF满足如下形式

其中归一化参数c定义如下，且 $\Gamma$ 是gamma函数.

$\small c=\frac{\Gamma (\sum_{k=1}^{K+1}\alpha_k)}{\prod_{k=1}^{K+1}\Gamma (\alpha_k)}$

或者另一种常见的定义形式如下：

$\small S_K =\{\mathbf x: 0\leq x \leq 1, \sum_{k=1}^{K}x_k=1\}$

$\small Dir(\mathbf x| \boldsymbol \alpha) =\frac{1}{B(\boldsymbol \alpha)}\coprod_{k=1}^{K}x_k^{\alpha_k-1}\mathbb{I}(\mathbf x \in S_K)\$

$\small {B(\boldsymbol \alpha)}=\frac{\coprod_{k=1}^{K}\Gamma (\alpha_k)}{\Gamma (\alpha_0)}=\frac{\Gamma (\alpha_1)\Gamma (\alpha_2)...\Gamma (\alpha_K)}{\Gamma (\alpha_1+\alpha_2+...+\alpha_K)}$

3.3 补充说明

在上面的定义中，向量 $\small X$ 的每一项 $\small x_1,...,x_K$ 都是一个概率值，各项求和小于或等于1：

$\small \sum_{k=1}^{K}x_k\leq 1$

如果我们想要一个概率向量恰好总和为 1，我们可以定义一个额外的概率:

$\small x_{k+1}=1-\sum_{k=1}^{K}x_k$ (1)

从而满足

$\small \sum_{k=1}^{K+1}x_{k}=1$ （2）

然而，没有办法严格定义向量的概率密度符合如下形式

$\small [x_1,...,x_K,x_{K+1}]$

因为(2)式中和为1的约束意味着在 $\small \mathbb{R}^{K+1}$ 空间上，除了 Lebesgue 测度等于 0 的子集外， $\small \mathbb{R}^{K+1}$ 空间上其它概率密度在任何地方都应该为零，而在这个特殊子集上，概率密度应该是无限的（这里涉及到 Dirac delta 函数相关内容）。

因此，处理概率总和为 1 的 K+1 个事件的正确方法是：

将上述定义的 Dirichlet 概率密度分配给 $\small K$ 个事件的概率为 $\small \(x_1,...x_K)$ 。
$\small x_{K+1}$ 按照等式 (1) 定义第 $\small K+1$ 个事件的概率。

3.4 详细推导

令 $\small Z_1,...Z_{K+1}$ 为gamma分布的 K+1 个独立随机变量，对应于gamma分布具有的均值分别为 $\small \alpha_1,..., \alpha_{K+1},$ , 自由度分别为 $\small \small 2\alpha_1,...,2\alpha_{K+1}$ ，定义

$\small S=\sum_{k=1}^{K+1} Z_k$

那么可以定义 $\small K\times 1$ 的随机变量

$\small X=[\frac{Z_k}{S},\frac{Z_2}{S},...,\frac{Z_k}{S}]^T$

那么我们就说随机向量 $\small X$ 具有以参数 $\small \small \alpha_1,..., \alpha_{K+1}$ 的Dirichlet 分布。

证明：Gamma 分布的随机变量是定义在区间 [0,1] 的正实数集合上的，也就是

$\small 0\leq \frac{Z_k}{S}\leq 1$

并且满足

$\small \sum_{k=1}^{K} \frac{Z_k}{S}\leq 1$

因此， $\small X$ 符合一个 Dirichlet 分布的随机变量所具有的定义域条件。而对于随机变量 $\small Z_k$ 的gamma分布，若具有均值为 $\small \alpha_k$ , 自由度为 $\small 2\alpha_k$ 的参数，则相应PDF表达式如下

既然变量 $\small z_1,...,z_K,z_{K+1}$ 是相互独立的，那么它们的联合PDF为

考虑一个一对一的变换（在前面一一对应的K项基础上再添加一个s项，构成K+1项）：

它的逆为：

$\small g^{-1}$ 的Jacobian 矩阵如下：

Jacobian 矩阵的行列式如下

1) 如果我们将前面的K行添加到第 K+1 行，行列式的值不会改变(线性变换不会改变行列式的值）

2) 三角矩阵(上三角或者下三角)的行列式等于其对角线元素的乘积。一对一变换的的联合PDF给出如下形式（作用在 $\small x_1,...,x_K,s$ ）：

积出s，我们得到

在上述最后步骤A中，我们使用了 Gamma 函数的定义。等式右边就是参数为 $\small \alpha_1,\alpha_2,...,\alpha_{K+1}$ 的狄利克雷分布的PDF。

特别注意(特别容易混淆）：对于概率值 $x_k$ 有和为1的条件

$\small \sum_{k=1}^{K+1}x_{k}=1$

但是对于参数向量 $\small \small \boldsymbol \alpha = (\small \alpha_1,\alpha_2,...,\alpha_{K+1})$ ，并没有和为1的条件

$\small \sum_{k=1}^{K+1}\alpha_{k} ?= 1$

3.5 一个特例

Beta 分布是 Dirichlet 分布的一个特例。

具体来说，如果我们定义上述维度K=1，那么

相应的PDF就变成

使用beta函数的定义形式

我们重写PDF为如下形式

这里的beta分布是以参数为 $\small \alpha_1,\alpha_2$ 关于随机变量 $\small X$ 的PDF。

4.数学解释

数学上的Dirichlet 分布是一种多变量连续概率分布，通常用于对未知概率向量的不确定性建模。注意：此处对向量建模反应了它的最根本特点，多项分布也是对向量建模。

Dirichlet 分布是 Beta 分布的多变量(多元)推广形式。对事件发生的概率 $X$ 建模。如果 $X$ 是未知的，我们可以把它当作一个随机变量，并分配给它一个 Beta 分布。如果是对互斥事件的未知概率向量建模，我们可以将 $X$ 视为随机向量，并为其分配一个 Dirichlet 分布。

Dirichlet 分布的参数可以视为是先验过程看到的次数，对应向量 $\small \boldsymbol \alpha=(\small \alpha_1,\alpha_2,...,\alpha_K)$ ,次数 $\small \alpha_i$ 对应的就是第 i 项(比如骰子的第i 面）先验过程已经出现的次数，那么次数越多，权重越大，意味这个面出现的可能性越高，在后验也就是投掷骰子试验中出现的可能性越大。

5. 可视化说明

根据上述分析，Beta 分布定义在[0,1]区间，由两个正形状参数参数化表示。如其所见，它是二项(包括伯努利)分布的共轭先验。下图显示了α和β具有不同取值的 Beta 分布PDF。

Dirichlet 分布视为是beta分布(两个形状参数)的推广，也就是两个参数推广为以向量 $\small \boldsymbol \alpha$ 为参数的Dirichlet 分布 $\small Dir(\boldsymbol \alpha)$ ，我们将这些形状概括为一个K-1单纯形。对于K=3，可视化其概率分布需要我们执行如下操作：

在我们的三角形上生成一组 x，y 坐标

将x，y坐标映射到二维单纯形坐标空间(二维对应就是三个点 K=3）

对每个点计算 $\small \small Dir(\boldsymbol \alpha)$

以骰子的投掷结果只有三个面为例(六个面不便于可视化展示)，选取不同的参数 $\small \small (\alpha_1,\alpha_2,\alpha_3)$ ，对应不同的狄利克雷分布，那么展示其概率分布(二维单纯形上的Dirichlet 分布)如下。

我们看到现在是控制分布的形状参数为向量 $\small \small \boldsymbol \alpha$ 。特别地，,控制分布的强度(达到峰值的尖锐程度）。如果对于所有 i 都有 $\small \alpha_i< 1$ ，我们会在单纯形的corner处得到尖峰。对于 $\small \alpha_i > 1$ ，概率分布趋向于单纯形的中心。随着不断增大，分布变得更紧密地集中在单纯形的中心周围。

在我们最初的骰子实验中，如果假定投掷的都是始终无偏的骰子，也就是 $\small \alpha_i \rightarrow \infty$ 。对于具有 $\small \alpha_i > 1$ 的对称 Dirichlet 分布，我们将生产出平均意义上的无偏骰子。如果目标是生产有负载(特意有偏)的骰子（例如，掷出 3 的概率更高），那么我们需要具有更高 $\small \alpha_3$ 的非对称 Dirichlet 分布。

6.实现代码

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import seaborn as sns
from math import gamma
from operator import mul
from functools import reduce
sns.set(style='white', font_scale=1.2, font='consolas')

def plot_mesh(corners):
    """Subdivide the triangle into a triangular mesh and plot the original and subdivided triangles."""
    triangle = tri.Triangulation(corners[:, 0], corners[:, 1])

    refiner = tri.UniformTriRefiner(triangle)
    trimesh = refiner.refine_triangulation(subdiv=4)
    
    plt.figure(figsize=(6, 4))
    for i, mesh in enumerate((triangle, trimesh)):
        plt.subplot(1, 2, i+1)
        plt.triplot(mesh)
        plt.axis('off')
        plt.axis('equal')    
        
        
class Dirichlet:
    """Define the Dirichlet distribution with vector parameter alpha."""
    def __init__(self, alpha):
        
        self._alpha = np.array(alpha)
        self._coef = gamma(np.sum(self._alpha)) / reduce(mul, [gamma(a) for a in self._alpha])
    
    def pdf(self, x):
        """Returns pdf value for `x`. """
        return self._coef * reduce(mul, [xx ** (aa-1) for (xx, aa) in zip(x, self._alpha)])
        
        
class PlotDirichlet:
    """
    Plot the Dirichlet distribution as a contour plot on a 2-Simplex.
    """
    def __init__(self, corners):
        self._corners = corners
        self._triangle = tri.Triangulation(corners[:, 0], corners[:, 1])
        # Midpoints of triangle sides opposite of each corner
        self._midpoints = [(corners[(i+1) % 3] + corners[(i+2) % 3]) / 2.0 for i in range(3)]
        
    def xy2bc(self, xy, tol=1.e-3):
        """Map the x-y coordinates of the mesh vertices to the simplex coordinate space (aka barycentric coordinates).
        Here we use a simple method that uses vector algebra. For some values of alpha, calculation of the Dirichlet pdf 
        can become numerically unstable at the boundaries of the simplex so our conversion function will take an optional 
        tolerance that will avoid barycentric coordinate values directly on the simplex boundary.        
        """
        s = [(self._corners[i] - self._midpoints[i]).dot(xy - self._midpoints[i]) / 0.75 for i in range(3)]
        return np.clip(s, tol, 1.0-tol)
        
    def draw_pdf_contours(self, ax, dist, label=None, nlevels=200, subdiv=8, **kwargs):
        """Draw pdf contours for a Dirichlet distribution"""
        # Subdivide the triangle into a triangular mesh
        refiner = tri.UniformTriRefiner(self._triangle)
        trimesh = refiner.refine_triangulation(subdiv=subdiv)
        
        # convert to barycentric coordinates and compute probabilities of the given distribution 
        pvals = [dist.pdf(self.xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]
    
        ax.tricontourf(trimesh, pvals, nlevels, **kwargs)
        #plt.axis('equal')
        ax.set_xlim(0, 1)
        ax.set_ylim(0, 0.75**0.5)
        ax.set_title(str(label))
        ax.axis('off') 
        return ax

        
if __name__ == '__main__':
    corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])
    plot_dirichlet = PlotDirichlet(corners)
   
    f, axes = plt.subplots(2, 3, figsize=(14, 8))
    ax = axes[0, 0]
    alpha = (0.85, 0.85, 0.85)
    dist = Dirichlet(alpha)
    ax = plot_dirichlet.draw_pdf_contours(ax, dist, alpha)

    ax = axes[0, 1]
    alpha = (1, 1, 1)
    dist = Dirichlet(alpha)
    ax = plot_dirichlet.draw_pdf_contours(ax, dist, alpha)

    ax = axes[0, 2]
    alpha = (5, 5, 5)
    dist = Dirichlet(alpha)
    ax = plot_dirichlet.draw_pdf_contours(ax, dist, alpha)

    ax = axes[1, 0]
    alpha = (1, 2, 3)
    dist = Dirichlet(alpha)
    ax = plot_dirichlet.draw_pdf_contours(ax, dist, alpha)

    ax = axes[1, 1]
    alpha = (2, 5, 10)
    dist = Dirichlet(alpha)
    ax = plot_dirichlet.draw_pdf_contours(ax, dist, alpha)

    ax = axes[1, 2]
    alpha = (50, 50, 50)
    dist = Dirichlet(alpha)
    ax = plot_dirichlet.draw_pdf_contours(ax, dist, alpha)


    f.savefig('../figures/dirichlet.png', bbox_inches='tight', transparent=True)

7.参考资料

Dirichlet distribution

The Dirichlet Distribution: What Is It and Why Is It Useful?

plot_dirichlet

simplex-distributions

Dirichlet Distribution, Dirichlet Process and Dirichlet Process Mixture

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)