Exact Bayesian Inference on Discrete Models via Probability Generating Functions: A Probabilistic Programming Approach
基于概率生成函数的离散模型精确贝叶斯推理:一种概率规划方法
https://proceedings.neurips.cc/paper_files/paper/2023/file/0747af6f877c0cb555fea595f01b0e83-Paper-Conference.pdf
摘要
我们提出了一种针对离散统计模型的精确贝叶斯推理方法，该方法能够为一大类离散推理问题找到精确解，即使这些问题具有无限支撑和连续先验。为了表达此类模型，我们引入了一种概率编程语言，该语言支持离散和连续采样、离散观测、仿射函数、（随机）分支以及对离散事件的条件设置。我们的关键工具是概率生成函数：它们为可通过程序定义的分布提供了紧凑的闭式表示，从而能够精确计算后验概率、期望、方差以及更高阶矩。我们的推理方法在名为Genfer的工具中得到了证明，该工具使用自动微分（特别是泰勒多项式）且无需计算机代数即可实现全自动化。实验表明，Genfer通常比现有的精确推理工具PSI、Dice和Prodigy更快。在一系列这些精确工具均无法解决的现实推理问题上，Genfer的性能可与近似蒙特卡洛方法相媲美，同时避免了近似误差。
1 引言
贝叶斯统计是一个在不确定性下进行推理的非常成功的框架，已广泛应用于人工智能/机器学习、医学与医疗保健、金融与风险管理、社会科学、气候科学、天体物理学等许多领域。其核心思想是将不确定性表示为概率，并通过贝叶斯定律根据观测数据更新先验信念，从而得出后验信念。贝叶斯统计中的一个关键挑战是计算这个后验分布：解析解通常不可能或难以求解，因此需要使用近似方法，如马尔可夫链蒙特卡洛（MCMC）或变分推理。在本文中，我们确定了一大类可以进行精确推理的离散模型，特别是计数数据的时间序列模型，如自回归模型、隐马尔可夫模型和切换点模型。我们通过利用概率生成函数（GFs）作为分布的表示来实现这一点。随机变量X的GF定义为函数G(x) := E[x^X]（这里X表示随机变量，x为其取值，E表示期望）。在概率论中，它是研究随机变量及其分布的一个众所周知的工具，与矩生成函数和特征函数相关[14，第4章]。在计算机科学中，GFs之前已被用于概率程序的分析[17, 4]以及某些类别图形模型的精确推理[26, 27]。
在这里，我们通过概率编程在更一般的上下文中统一应用它们，使得能够在更具表达力的贝叶斯模型上进行精确推理。
我们借助一种概率编程语言来刻画所支持的模型类别。
概率编程[25]最近已成为贝叶斯推理中的一个强大工具。概率编程系统允许用户以精确且直观的方式将复杂的统计模型指定为程序，并自动化贝叶斯推理任务。这使得从业者能够专注于建模，而将通用推理算法的开发留给专家。因此，诸如Stan[3]之类的概率编程系统在统计学家和数据科学家中越来越受欢迎。我们描述了一种名为SGCL（统计受保护命令语言）的编程语言，它扩展了pGCL（概率GCL）[17]。该语言经过精心设计，既简单又富有表现力，且限制性恰到好处，能够对用其编写的所有程序进行精确的贝叶斯推理。
贡献 我们为离散贝叶斯模型上的精确推理提供了一个新框架。
（a）我们的方法适用于一大类具有无限支撑的离散模型，特别是计数数据的时间序列模型，如人口动态的自回归模型、贝叶斯切换点模型、混合模型和隐马尔可夫模型。据我们所知，在之前，除了最简单的种群模型外，没有其他已知的确切推理方法。
（b）模型是在一种概率编程语言（PPL）中指定的，这提供了模型定义的灵活性，从而促进了模型构建。我们的PPL支持随机分支、连续和离散先验、离散观测以及对涉及离散变量的事件进行条件设置。
（c）用该语言编写的每个程序都可以自动且正确地转换为代表后验分布的生成函数。从这个生成函数中，可以提取后验均值、方差和更高阶矩，以及后验概率质量（对于离散分布）。
（d）我们构建了一个名为Genfer（“用于推理的生成函数”）的优化工具，它接受概率程序作为输入，并自动计算上述关于后验分布的描述性统计量集合。
（e）我们证明了（1）在具有有限支撑的基准测试中，Genfer的性能通常优于现有的精确推理工具，以及（2）在一系列现有精确工具不支持的现实世界示例中，Genfer与近似蒙特卡洛方法相媲美，同时实现了零近似误差。
相关工作 计算概率程序的确切后验分布通常是不可行的，因为它需要积分的解析解[9]。因此，现有系统要么限制编程语言只允许可处理的构造（这是我们的方法），要么无法保证成功推理。前者包括Dice[16]（仅支持有限离散分布）和SPPL[22]（支持一些无限支撑分布，但要求有限离散先验）。它们都基于概率电路[5]或其扩展，这些允许高效且精确的推理。后者包括PSI[9]和Hakaru[21]系统，它们不施加语法限制。它们依靠计算机代数技术来找到后验的闭式解。这样的形式在一般情况下可能不存在，即使存在，系统也可能找不到它，并且运行时间不可预测且不可扩展[22]。我们评估中的案例研究（第5.2节）都无法通过上述工具处理。
概率生成函数是研究具有无限支撑的随机变量的概率论中的有用工具，例如在分支过程的上下文中[7，第12章]。在贝叶斯推理的上下文中，概率生成电路（PGC）利用它们来提高概率电路的表达能力[28]，并实现高效推理[15]。然而，PGC无法处理具有无限支撑的随机变量，这是我们的方法所关注的重点。生成函数也被应用于概率图形模型：Winner和Sheldon[26]为泊松自回归模型找到了生成函数的符号表示，并从中提取了后验概率。随后，Winner等人[27]将该模型扩展到除泊松分布之外的其他潜在变量分布，其中符号操作不再可行。相反，他们使用自动微分来评估生成函数及其导数，这使得能够为图形模型进行精确推理。概率编程是图形模型的一种优雅泛化，允许更丰富的模型表示[10, 25, 12]。我们的贡献是通过概率编程为贝叶斯模型上的精确推理提供了一个新框架。
在离散概率程序未进行条件约束的上下文中，Klinkenberg等人[17]使用生成函数（手动）分析循环不变式并确定终止概率。在后续工作中，Chen等人[4]扩展了这些技术，以在特定限制下自动检查循环程序是否生成了指定的分布。他们的分析工具Prodigy最近还获得了执行贝叶斯推断的能力[18]。它支持具有无限支持的离散分布（但不支持连续先验），并且由于其依赖于计算机代数，因此在可扩展性上不如我们的自动微分方法（见第5节）。
局限性：对于仅涉及有限离散分布的概率程序，精确的后验推断已经是PSPACE-hard问题[16，第6节]。因此，我们的方法不可能始终具有高性能，并且必须限制所支持的概率程序类别。事实上，我们的编程语言禁止某些构造，如非线性变换和连续变量的观测，以保持生成函数的封闭形式表示（见第2节）。出于同样的原因，我们的方法无法计算连续参数的概率密度函数（但可以计算离散参数的概率质量和所有参数的精确矩）。
在性能方面，推断的运行时间与程序中观察到的数字数量成多项式关系，与程序变量的数量成指数关系（见第4节）。尽管存在这些局限性，我们的方法首先证明了精确推断是可能的，并且我们的评估（见第5节）表明，它支持许多现实世界中的模型，并且在实践中能够高效地进行精确推断。
2 贝叶斯概率编程
概率编程语言通过添加两个额外的构造来扩展普通编程语言：一个用于从概率分布中采样，另一个用于对观测值进行条件约束。我们首先通过一个基于种群生态学的简化示例（参见[26]）来讨论可以在我们的语言中编写的简单程序。
假设你是一位生物学家，试图估计迁移到新栖息地的动物种群数量。动物的迁移通常使用泊松分布进行建模。由于无法彻底计数（否则我们就不需要估算技术），因此我们假设每个个体被独立观测到的概率是一定的；也就是说，计数服从二项分布。为了简化，我们假设泊松分布的速率和二项分布的概率是已知的，分别为20和0.1。（对于更现实的种群模型，请参见第5节。）作为生成模型，这可以表示为XPoisson(20)；YBinomial(X, 0.1)。
3 生成函数
3.1 将程序翻译成生成函数
对于条件语句，如果Xk ∈ A则执行{P1}，否则执行{P2}，我们将生成函数G拆分为两部分：一部分是条件满足的情况（GXk∈A），另一部分是其余情况（G-GXk∈A）。前者通过then分支的JP1K gf进行变换，后者通过else分支的JP2K gf进行变换。理解GXk∈A的计算最好的方法是，将G视为幂级数，其中我们只保留xk的指数在A中的项。从具有恒定参数的分布中采样Xk ∼ D的操作，首先是通过边缘化Xk，然后乘以参数为xk的分布D的生成函数。
第一个新的构造是从复合分布中采样Xk ∼ D(Xj)（详见附录B）。观察事件Xk ∈ A与条件语句中一样，使用GXk∈A，如上文所述。就像程序定义的次概率分布必须在最后一步进行归一化一样，我们也必须对生成函数进行归一化。归一化常数是通过边缘化所有变量来计算的：G(1, ..., 1)。因此，我们通过使用逆值进行缩放，获得了表示归一化后验分布的生成函数：normalize(G) := G / G(1, ..., 1)。这些直觉可以在以下定理的形式中变得严谨，该定理在附录B.2中得到了证明。
新颖性：该语义建立在[17, 4]的基础上。据我们所知，复合分布Poisson(λ · Xj)和Bernoulli(Xj)的生成函数（GF）语义是新颖的，并且前者是支持第5节中大多数模型所必需的。虽然已有研究在特定模型[26]的上下文中考虑了观测值的生成函数，但在此之前的概率编程语言通用上下文中尚未有过此类研究。更广泛地讲，以往涉及生成函数的工作仅考虑了离散分布，而我们则允许从连续分布中进行采样。这是一项重大的推广，需要采用不同的证明技术，因为[26, 17, 4]中的证明所依赖的幂级数表示P{i∈N} P[X = i]x^i对于连续分布而言是无效的。
4 实现与优化
实现生成函数（GF）语义的主要难点在于计算偏导数。一种自然的方法（如[4, 18]所遵循的）是操作生成函数的符号表示，并使用计算机代数来计算导数。然而，Winner等人[27]的研究表明，这种方法通常扩展性较差，因为生成函数的大小通常随着所依赖的数据量的增加而迅速增长。之所以如此，是因为程序中的每个观察语句observe Xk = d都会被转换为d阶偏导数。由于概率程序往往包含许多数据点，因此总导数阶数达到数百是很常见的。函数符号表示的大小通常会（而且通常会）随着导数阶数的增加而呈指数级增长：两个函数f · g的乘积的导数是两个乘积f' · g + f · g'的和，因此表示的大小会加倍。因此，运行时间将是，其中d是所有观察值的总和，这显然是不可接受的。
相反，我们利用这样一个事实，即我们不需要生成生成函数的完整表示，而只需要评估它及其导数。我们为此实现了自己的自动微分框架，因为现有的框架并不适用于计算阶数大于4或5的导数。事实上，使用泰勒多项式而不是直接计算高阶导数更为高效。Winner等人[27]已经为总体模型（只有一个变量）做到了这一点，而我们则将其扩展到包含多个变量的更一般的环境中，这需要多元泰勒多项式。
在这种方法中，导数是容易的部分，因为可以从泰勒系数中读取它们，但泰勒多项式的组合是瓶颈。Winner等人[27]使用了一种朴素的方法，这对于他们的单变量用例来说已经足够快。
运行时间 对于 n 个变量，泰勒多项式的简单组合需要 时间，其中 d 是程序中所有观测值的总和，即微分的总阶数，也就是多项式的度数。请注意，这在 d 上是多项式级别的，与符号方法不同，但在变量数量 n 上是指数级的。这并不像看起来那么糟糕，因为在许多情况下，变量的数量可以保持在一两个，与数据点的值（如第 5 节中的模型）相比。实际上，我们利用生成函数的具体组合结构，在最坏情况下，对于 n = 1 时实现
，对于 n ≥ 2 时实现
，而在实践中通常更快。总的来说，我们的实现在最坏情况下需要
时间，其中 s 是程序中的语句数量，d 是所有观测值的总和，n 是程序变量的数量（见附录 C.3）。
减少变量数量 鉴于在 n 上的指数级运行时间，编写概率程序时减少变量数量至关重要。一方面，程序中不再需要的变量通常可以稍后重用。此外，赋值和采样可以与加法结合：Xk+=... 表示 Xn+1 := ...; Xk := Xk+Xn+1 和 Xk +∼ D 对于 Xn+1 ∼ D; Xk := Xk + Xn+1。这些语句的生成函数可以很容易地计算，而不需要引入临时变量 Xn+1。类似地，在观察语句和 if 条件中，我们使用简写 m ∼ D 与 m ∈ N 表示事件 Xn+1 = m，其中 Xn+1 ∼ D。概率程序通常包含许多这样的观察，特别是来自复合分布的观察。因此，优化生成函数以避免这个额外的变量 Xn+1 是值得的。Winner 和 Sheldon [26] 可以在特定模型的背景下为复合二项分布避免一个额外的变量。我们将其扩展到我们更通用的设置，并为连续变量，以及复合泊松、负二项和伯努利分布提供优化的语义。实际上，这种优化对于第 5 节中的许多示例实现良好性能至关重要。优化的翻译及其正确性证明可以在附录 C.2 中找到。
实现：我们的工具Genfer可以读取程序文件，并输出指定变量的后验均值、方差、偏度和峰度。对于在N上取值的离散变量，它还可以计算达到可配置阈值的后验概率质量。为了严格控制性能，特别是针对多元泰勒多项式，Genfer使用Rust[20]编写，Rust是一种安全的系统编程语言。我们的实现在GitHub上可获取：github.com/fzaiser/genfer
数值问题：Genfer可以使用多种数字格式进行计算：64位浮点数（默认且最快）、用户指定精度的浮点数和有理数（如果不出现无理数）。为了确保浮点结果数值稳定，我们还实现了区间算术来限制舍入误差。最初，我们发现包含连续分布的程序会导致灾难性的抵消误差，这是由于它们的生成函数（GFs）中的对数项导致的，这些对数项的泰勒展开表现不佳。我们通过略微修改生成函数的表示方式，避免了使用对数，从而解决了这个问题（详细信息见附录C）。对于第5.2节中的所有示例，我们的结果精确到至少5位有效数字。
5 实证评估
5.1 与精确推理方法的比较
我们将我们的工具Genfer与以下用于精确贝叶斯推理的工具进行比较：使用加权模型计数的Dice[16]、使用计算机代数操作密度函数的PSI[9]，以及基于生成函数（如Genfer）但使用计算机代数而非自动微分的Prodigy[18]。我们在PSI基准测试集[9]上对这些工具进行了评估，排除了PSI独有的测试（例如，由于连续分布的观察）。大多数基准测试仅使用有限离散分布（标记为“F”），但有三个测试具有连续先验（标记为“C”）。
我们测量了每个工具的实际推理时间（不包括启动时间和输入文件解析时间），并记录了在一小时超时限制下连续五次运行中的最小值。Dice和Genfer默认使用浮点数（FP），而Prodigy和PSI使用有理数，这虽然更慢但可以防止舍入误差。为了公平比较，我们在有理数模式下评估了所有工具，并单独比较了Dice和Genfer在浮点数模式下的表现。结果（表3）表明，即使在有限离散模型上，Genfer的速度也很快，尽管我们的主要关注点是具有无限支持的模型。
5.2 与近似推理方法的比较
在具有无限支持的现实模型上，我们的方法无法与其他精确推理方法进行比较（见第5.3节）：可扩展系统Dice[16]和SPPL[22]不支持此类先验，而符号求解器PSI[9]和Prodigy[18]在运行一小时后内存溢出或超时。
截断：作为一种替代方法，我们考虑了通过截断具有无限支持的离散分布来近似后验。这将问题简化为更适用于精确技术的有限离散推理。但我们没有采用这种方法，因为Winner和Sheldon[26]已经在他们的图形模型上证明了其相对于基于生成函数的精确推理的劣势。此外，截断通用概率程序更加困难，对于连续先验甚至是不可能的。
蒙特卡洛推理：因此，我们将我们的方法与蒙特卡洛推理方法进行了比较。具体来说，我们选择了Anglican[24]概率编程系统，因为它为具有许多最新推理算法的离散模型提供了最佳内置支持。其他流行的系统不太适合：Gen[6]专门用于可编程推理；Stan[3]、Turing[8]和Pyro[1]主要针对连续模型；而WebPPL[11]的离散推理算法不如Anglican广泛（例如，不支持交互粒子MCMC）。
方法论：蒙特卡洛推理算法的近似误差取决于其设置（例如，SMC的粒子数量），并且随着运行时间的延长和样本数量的增加而减小。为了确保公平比较，我们使用以下设置：对于每个推理问题，我们使用不同的配置（设置和采样预算）运行多个推理算法，并报告近似误差和经过的时间。为了衡量近似的质量，我们报告了精确解和近似后验分布之间的总变化距离（TVD），以及后验均值µ、标准差σ、偏度（第三标准化矩）S和峰度K的近似误差。为了确保我们的误差度量在分布的平移和缩放下具有不变性，我们计算均值的误差为|ˆµ−σµ|，其中µ是真实后验均值，ˆµ是其近似值，σ是真实后验标准差；标准差的相对误差（因为σ在平移下不变但在缩放下不变）；以及偏度和峰度的绝对误差（因为它们在平移和缩放下都不变）。为了减少噪声，我们对这些误差度量和计算时间进行了20次运行的平均，并报告了标准误差作为误差条。我们运行了Anglican中实现的几种知名推理算法：重要性采样（IS）、轻量级Metropolis-Hastings（LMH）、随机游走Metropolis-Hastings（RMH）、序贯蒙特卡洛（SMC）、粒子Gibbs（PGibbs）和交互粒子MCMC（IPMCMC）。为了可比性，在所有实验中，每个算法都使用两个采样预算（如果可能的话），以及两个不同的设置（一个是默认值），总共四个配置。采样预算为1000或10000，因为显著低于此的样本量会产生不可用的结果，而显著高于此的样本量则比我们的精确方法花费更多时间。我们丢弃了前20%的样本，这是一种称为“预热”的标准程序。
请注意，此设置对近似方法是有利的，因为我们只报告了每个配置一次运行的平均时间。然而，在实践中，人们并不知道最佳配置，因此需要以不同的设置多次运行算法。相比之下，我们的方法只需要运行一次，因为结果是精确的。
5.3 具有无限支撑分布的基准
种群生态学 我们的第一个基准来自[26, 27]，用于模拟动物种群。
我们在示例2.1中已经看过一个简化版本。在这里，我们模拟在时间步长k = 0, ..., m时的人口数量Nk。在每个时间步长，都有符合泊松分布的新到达个体数量，这些新个体被添加到前一个时间步长中符合二项分布的存活个体数量中。每个个体以固定的概率δ被观察到，因此观察到的个体数量符合二项分布：
其中，模型参数λk ∈ R, δ ∈ [0, 1]取自[26]；检测概率ρ被设定为0.2（[26]考虑了一系列的值，但出于空间考虑，我们只选择了一个）；在时间步长k时观察到的种群个体数量yk ∈ N是根据与[26]相同的真实情况模拟的。目标是推断种群中个体的最终数量Nm。我们将影响运行时间的种群规模模型参数和观察到的值设定为[26]中最大值的4倍（详见附录D.3）。结果（图1a）显示，由于我们的方法是精确的，因此在计算时间和准确性方面均优于MCMC方法。
修改 虽然这个模型已经在[26]中得到了精确求解，但我们的概率编程方法使得修改模型变得轻而易举，因为整个推理过程是自动化的，只需要改变程序的几行代码即可：（a）我们可以用一个条件来模拟自然灾害影响后代率的可能性：灾难~伯努利(0.1)；如果灾难=1，则{Newk~泊松(λ')}，否则{Newk~泊松(λ)}，或者（b）我们不再只模拟一个种群，而是可以模拟两种相互作用的个体种群，即一个多类型分支过程（见图1c）。这些修改中的任何一个都不能由[26]或[27]处理。第一个修改的结果（图1b）非常相似。更复杂的第二个修改虽然需要更长的时间来精确求解，但比使用10000个样本的近似推理花费的时间更少，并且实现了零误差。
转换点模型 我们的第二个基准是贝叶斯转换点分析，它涉及检测某些事件频率随时间的变化。我们使用[23]中的模型，该模型具有连续先验，以及关于煤矿事故频率的111个真实世界数据点。我们比较了矩误差（图3）和总变差距离（TVD）误差（图2a）。在这两种情况下，近似方法的准确性都较低，且耗时更长，不如我们的精确方法。
混合模型 我们考虑在同一个数据集上使用二元混合模型，混合权重相等，速率为几何先验：每个数据点都是从两个具有不同速率的泊松分布的混合中观察到的，任务是推断这些速率。由于混合模型的多模态性，近似推理方法对其来说非常困难，这在图2b中得到了证实。即使误差最低的运行也只能覆盖两种模式中的一种（见图2c中的样本直方图）。
隐马尔可夫模型 我们使用基于[22, 第2.2节]的隐马尔可夫模型，但涉及无限（几何）先验。它是一个具有已知转移概率的两状态系统，观察到的数据的速率取决于隐藏状态。我们在30个模拟数据点上运行了这个模型。对于这个模型，我们的方法也明显优于近似方法（图2d）。
据我们所知，除了[26]中未经修改的第一个问题外，我们的方法是首次为这些问题找到精确解。为了简洁起见，我们只介绍了这些基准的最重要方面，将其作为概率程序的编码留给了附录D.3。代码和复现说明见补充材料。
6 结论
通过利用生成函数，我们已经开发并验证了一个框架，该框架能够对离散模型进行精确的贝叶斯推理，即使这些模型具有无限支撑和连续先验。我们在一种富有表现力的概率编程语言中指定了一系列模型，并使用我们的工具Genfer对这些模型进行了自动处理，展示了其竞争性的性能。
未来工作
一个自然而然的问题是，如何将我们的方法与更通用的概率编程系统相结合。例如，在基于采样的推理算法中，如果可能的话，可以想象使用生成函数来精确求解子程序。更广泛地讲，我们希望能够探索如何改进GF（生成函数）翻译的组合性。
目前，我们的GF翻译描述了所有程序变量的联合分布——我们从未对变量的子集进行过“局部”推理。采用组合方法可能会促进GF方法在诸如Anglican等函数式概率语言中的应用。通过允许对程序的较小部分进行精确的推理和分析，然后再将这些结果组合起来以形成对整个程序的更广泛理解，组合性可以显著提高推理的效率和可扩展性。
B 关于生成函数的详细信息
可以在表 4a 和 4b 中找到支持的分布及其生成函数的详尽列表。
同样地，计算连续分布的后验密度也是值得期待的。事实上，存在通过拉普拉斯反变换从生成函数中恢复概率密度函数的数学方法。然而，这在实践中无法实现自动化，因为它需要求解积分，这是难以处理的。
实现细节 我们的工具Genfer是用Rust编程语言[20]实现的，Rust是一种安全的系统编程语言。选择Rust的主要原因是出于对其低级控制和性能的考虑：评估生成函数导数的泰勒多项式运算需要快速执行，并且已经针对由概率程序产生的GF（生成函数）的结构进行了优化（参见附录C.3）。虽然C或C++同样可以满足性能标准，但Rust的语言特性，如内存安全性、枚举（标记联合）和模式匹配，使得实现过程更加愉快、稳健且易于维护。第一作者对Rust的熟悉程度也是促成这一选择的因素之一。
泰勒多项式的系数被存储在由ndarray库（注：ndarray可能是指一个具体的库，但在此处作为示例，实际名称可能有所不同）提供的多维数组中。任意精度的浮点数和无界有理数由rug库（注：rug同样可能是指一个具体的库，用于提供高精度数学运算，但在此处仅为示例）提供，该库是GNU库GMP（用于有理数）和MPFR（用于浮点数）的接口。我们的实现在GitHub上可获取：github.com/fzaiser/genfer。
记忆中间结果 在计算函数G的泰勒展开时，如果G的子表达式的泰勒展开在计算过程中多次出现，那么对这些中间泰勒展开结果进行记忆化（缓存）是非常重要的。例如，条件语句被转换为使用先前生成函数两次的生成函数。如果反复评估这个生成函数，那么每个条件语句的运行时间都会翻倍。通过记忆化，可以在合理的时间内处理包含大量分支的程序（例如，在混合模型中路径数超过，参见表5）。
https://proceedings.neurips.cc/paper_files/paper/2023/file/0747af6f877c0cb555fea595f01b0e83-Paper-Conference.pdf
