为什么贝叶斯统计量需要蒙特卡罗方法

这篇帖子出自一系列围绕Twitter评论的问题，该评论提出了一些非常有趣的观点，涉及贝叶斯假设检验是如何工作的，以及分析解决方案甚至无法解决贝叶斯统计学中一些看似微不足道的问题。比较Beta分布随机变量是本博客(我的书中也是如此)经常出现的事情。设置相当简单：将A/B测试建模为来自两个Beta发行版的样本，从每个发行版抽取大量样本，然后比较结果。这种模拟方法最初常常看起来像是解决更复杂数学问题的一个聪明的小把戏，但实际上是蒙特卡洛积分的原始形式，并被证明是真正解决这个问题的唯一方法之一。通过在这篇文章中更深入地探讨这个话题，我们将看到许多人关于解析解的一些迷思，以及为什么蒙特卡罗方法对贝叶斯统计如此重要。

最近在推特上发生了一场有趣的对话。首先，我转发了一条关于内特·西尔弗(知名作家和选举预报员)的推文，发布了他对2020年总统选举的最新预测，显示拜登获胜的可能性为71%，而特朗普的获胜概率为29%。

我们的预测上升了！它认为乔·拜登获胜的几率为71%，唐纳德·特朗普获胜的几率为29%。Https://t.co/ajG88SznSA。

-Nate Silver(@NateSilver538)2020年8月12日。

在上次选举之后，内特·西尔弗的预测受到了相当多的批评，因为他的公司(538)预测希拉里·克林顿在2016年总统大选中获胜的可能性为71%。这种批评总是让我个人感到恼火，因为从统计数字来看，71%的人通常被认为不是对任何事情都有强烈的信念。因此，有人相信候选人有71%的成功机会，但他们仍然失败，这并不矛盾，也不令人惊讶。即使是在研究典型的p值时，我们也会等待95%的确定性，然后才会做出断言(许多人认为这是一个相当薄弱的信念)。但由于某种原因，每当选举民调出现时，似乎就连非常有统计意识的人都会突然认为51%的可能性是一个很高的概率。我转发了内特·西尔弗(Nate Silver)的预测，提到了我的烦恼，并提供了另一个胜算相似的案例的例子：

这一次，我们都能记住，在统计学上，我们很少会将P(H|D)=0.71作为对任何事情的强烈信念。作为比较，如果在A/B测试中，我们有这些结果：A在15次试验中有2次成功，B在14次试验中有3次成功，这大致是我们对B的信念有多强https://t.co/bB4PiB5Tao。

-威尔·库尔特(@will kurt)2020年8月12日。

如果您正在运行A/B测试，并且您的A变体在15次试验中有2次成功，而您的B变体在14次试验中有3次成功，那么您将有大约71%的信心认为B是更好的变体。即使是没有接受过太多统计培训的人也可能会非常怀疑这样的说法，但不知何故，在选举预测期间，即使是经验丰富的统计学家也能看到西尔弗的帖子，认为拜登获胜是板上钉钉的事。

推特用户@mbarras_ing问了一个非常重要的后续问题，要求解释这个结果：

我可能有点慢，但你能详细说明一下吗？如何根据信息计算0.71，A在15次试验中有2次成功，B在14次试验中有3次成功？

-马修·里斯·巴拉斯(Matthew Rhys Barras)(@mbarras_ing)2020年8月13日。

我非常坚信，所有的统计数字，即使是即席估计的数字，都应该是可重复和可解释的。我在这篇博客和我的书中都写了相当多关于处理类似问题的文章。总体情况是，我们将对A和B转换用户的速率进行参数估计，然后计算B大于A的概率。由于我们正在重新估计转换率，因此我们将使用Beta分布作为参数估计的分布。在本例中，我还假设我们的A和B变体在\(\text{Beta}(1，1)\)之前。A的可能性为\(\text{Beta}(2，13)\)，B的可能性为\(\text{Beta}(3，11)\)，因此我们可以将A和B表示为来自这些后验的两个随机变量样本：$$A\sim\text{Beta}(2+1，13+1)$B\sim\text{Beta}(3+1，11+1)$$。

最后，我们可以查看结果来计算B大于A的概率：

在这种情况下，内特·西尔弗的问题得到了0.7028，接近71%。

这解释了我们的概率是从哪里得到的，但当您看到这个结果时，有一个明显的问题会出现，这个问题是由@Little_Rocko提出的。

这是一个很棒的例子。有没有一种简单的方法(就像非暴力)来找到与概率匹配的贝塔参数？

-Rocko(@Little_Rocko)2020年8月16日。

这个问题很有趣，因为它提出了两个不一定相关的问题：-有没有简单的方法来解决这个问题-有没有一种非暴力的方法来解决这个问题，然后再继续下去我想说的是，这一小段R涉及的抽象比乍看起来要抽象得多。我们在这里真正做的等同于使用蒙特卡罗模拟来积分两个Beta分布随机变量之间的差值的分布。在接下来的两个部分之后，将会更清楚地看到，这里正在发生的是一个令人惊讶的复杂的操作，在我看来，这是解决这个问题的最简单的方法，也不是真正的暴力解决方案。

当我们看到数学问题的计算解决方案时，我们的第一反应通常是觉得我们在避免用解析的方式解决问题。解析解是使用数学分析来找到闭合形式解的解。举一个明显的例子，假设我想要找到使R中的\(f(X)=(x+3)^2\)最小化的值，我可以通过查看一系列答案来强行实现这一点，如下所示：

不出所料，我们得到的答案是-3，但是这个解决方案需要一些工作：我们需要知道如何编码，并且我们需要一台计算机。这也有点乱，因为如果我们迭代的增量没有精确地包含-3(比如0.031)，我们就不会得到确切的答案。如果我们知道一些基本的微积分，我们就知道我们的最小值必须是导数为0的地方。我们可以很容易地计算出$$f&(X)=2(x+3)$$，当$$x=-3$$时，$$2(x+3)=0$$。但即使是微积分的一部分也是困难的，经常一次解决它会让未来的解决方案变得容易得多。举个例子，如果你想找出均值为\(\µ\)、标准差为\(\sigma\)的正态分布的最大似然。要解决这个问题，我们先从正态分布的pdf开始：$$\varphi(X)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}$$现在计算它的导数不一定很容易。这当然是我们能做的事。我们真正关心的是$$\varphi'；(X)=0$$，我们可以发现发生在什么时候(计算派生当然留给读者练习)：$$\frac{\mux}{\sigma}$$。

这使我们认识到一个令人惊讶的事实，即对于我们遇到的任何正态分布，我们知道最大似然估计是样本在\(x=\µ\)！尽管我们的微积分可能需要我们做一些工作，但是一旦我们完成了这项工作，对任何正态分布进行最大似然估计的问题就真的变得很容易了！

让我们重新审视我们原来的问题，这一次试图找到一个解析的解决方案。这是一个非常有趣的案例，因为可以说这是你能想象到的最简单的贝叶斯假设检验。回想一下，我们在每次测试中都有两个随机变量代表我们的信念。这些是按照我们前面描述的后部分布的。$$A\sim\text{Beta}(2+1，13+1)$B\sim\text{Beta}(3+1，11+1)$$我在这里跳过了推理中的一些步骤。我们想知道的是：$$P(B&>A)$$。

它没有用一种特别有用的数学方式来表达。解决此问题的更好方法是将其视为两个随机变量的和(在本例中为差)。我们真正想知道的是：$$P(B-A>；0)$$为了解决这个问题，我们可以考虑一个新的随机变量\(X\)，它将是B和A之间的差：$$X=B-A$$最后我们假设我们有一个概率密度函数，我们将调用\(\text{pdf}_X\)。如果我们知道\(\text{pdf}_X\)我们的解决方案非常接近，我们只需要在0和此分布的最大域之间进行积分：$$P(B>；A)=P(B-A>；0)=\int_{x=0}^{\text{max}}\text{pdf}_{X}(X)$$这看起来有点复杂，但前面还有一个大问题。与正态分布随机变量不同，我们没有贝塔分布随机变量的正态和定理的等价物(我们将在稍后讨论这一点)。\(\text{pdf}_X\)是什么样子？首先，我们知道它本身不是Beta发行版。我们可以看到这一点，因为我们知道此分发的域(或支持)不是\([0，1]\)。因为它们是Beta分布的，所以A和B都可以采用从0到1的值，这意味着这种差异的最大结果是1，但最小结果是-1。所以不管这个分布是什么，它的域是\([-1，1]\)，这意味着它不可能是Beta分布。我们可以使用关于随机变量和的各种规则来确定该分布的均值和方差，但是如果不知道该分布的确切形式，我们就无法解析地求解该积分。在这里我们可以

这就引出了问题的下一个要点，那就是一种非蛮力方法。来自计算机科学背景的我通常会认为暴力是指在找到最好的解决方案之前不择手段地尝试解决方案。对于这个问题，该问题的代码将涉及通过概率进行迭代，然后提出一些方法来分析结果。我们在模拟中所做的实际上是蒙特卡罗积分的一种形式。我们已经确定\(P(B>；A)\)等同于\(\int_{x=0}^{\text{max}}\text{pdf}_{x}(X)\)。由于\(\text{pdf}_{X}\)未知，我们通过从B和A之间可能的差异中抽样来近似它。如果可能有助于可视化此分布：

现在我们可以在所有答案的比例中加上阴影，在我们的案例中，这些答案将由sum(b_sample>；a_sample)/N捕获，我们可以看到我们实际上是在做积分。

因此，即使我们只有一点R码，而且我们做了大量的计算，我们仍然相当有效地将整数\(\text{pdf}_{X}\)从0逼近到1。这与蛮力方法有很大的不同，后者尝试候选解，直到找到它满意的解。同样重要的是要注意，这不是一个巧妙的破解方法，而是解决这个问题的非常可靠的数值方法。它可能不像求解适当的积分那样准确和高效(感谢@FinnLindgren)，但它非常容易实现、理解并建立在坚实的数学基础上。随着我们转向日益复杂的贝叶斯模型，我们继续以类似的方式解决问题，但是由于集成本身变得更加棘手，我们需要更复杂的工具，如Stan和PyMC3。

到目前为止，我们已经为我们的解决方案的这种抽样方法进行了很好的辩护，但在实践中，我们仍然有很多次甚至无法负担运行简单模拟的费用。举个例子：在我最初的一份数据科学家工作中，我们想使用这些贝叶斯技术实现一个A/B测试工具，但计算是用JavaScript完成的(这意味着我们没有rbeta函数)，而且用户不喜欢在他们的结果中有任何随机性(即使这不会影响最终的结论)。当我们不想在这种情况下使用包含正态分布的两个非常好的性质的模拟时，我们可以做一个聪明的分析近似：正态分布近似贝塔分布的能力和正态和定理。

Beta分布定义为有两个参数\(\α\)和\(\β\)，其中\(\α+\β=n\)是观测总数。结果表明，随着(N)向初值(比这要早得多)和(α近似β)的方向发展(并且不需要那么接近)，Beta分布可以很好地用正态分布来近似。由于正态分布是用均值和方差(或标准差)参数化的，要进行这种近似，我们需要知道关于Beta分布的两条信息：它的期望和方差。这样就可以很容易地确定近似的正态分布。Beta分布的期望值为：$$E[\text{Beta}(\alpha，\beta)]=\frac{\alpha}{\alpha+\beta}$$，方差定义为：$$Var[\text{Beta}(\alpha，\beta)]=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$$这意味着我们现在可以通过以下方式近似任何Beta分布：$$。\beta)\约\mathcal{N}(\u=\frac{\alpha}{\alpha+\beta}，\sigma=\sqrt{\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}})$$我们可以通过比较Beta(50，50)和\(\mathcal{N}(0.5，0.05)\)来查看此处的示例。

回想一下，添加Beta分布随机变量的问题是，我们没有很好的方法来获得这个新随机变量的PDF，它表示它们的差异。值得庆幸的是，正常的分发没有这个问题。把两个正态分布的随机变量加起来，就会得到一个新的正态分布的随机变量，它的均值和标准差都很容易计算。正态和定理告诉我们：给定两个独立的正态分布的随机变量：$$A\sim\Mathcal{N}(\Mu_A，\sigma_A)$B\sim\Mathcal{N}(\MU_b，\Sigma_B)$$这些变量的和的分布定义为：$$A+B=\Mathcal{N}(\MU_A+\MU_B，\sqrt{\Sigma_A^2。

A的后验pdf为\(\text{Beta}(3，14)\)，这意味着我们可以将其近似为：$$\mathcal{N(\µ=0.18，\sigma=0.09)}$$我们可以看到此分布与后验分布的近似程度：

正如我们在这里看到的，这个近似值是..。好吧。这主要是因为\(\alpha+\beta\)仍然相当小。尽管如此，这仍然是我们提供的最好的工具

考虑到A的近似值不是很大，所以这个差值分布的近似值也不是很好也就不足为奇了。然而，最终的考验是概率是否相同。所有这些近似工作的好处是，我们可以给出一个清晰的积分。我们将使用R'；的集成函数来实现这一点，但是pnorm也可以工作(如果您愿意，您甚至可以对其进行z归一化，然后在旧的概率书中查找它！)。

我们得到的概率是0.64，考虑到我们的近似是多么粗略，这是相当接近的。因此，只要这项技术被用于估计，比如改善的可能性，它就是一个相当不错的把戏。随着观测次数的增加，蒙特卡罗积分和正态近似之间的差异将变得非常小。尽管如此，我们仍然可以看到蒙特卡罗方法是我们对解决方案最诚实的表示。

贝叶斯A/B检验的案例几乎是人们能为真实世界的贝叶斯假设检验想出的最简单的例子。让我着迷的是，即使在这个直截了当的例子中，数学细节也出人意料地复杂。尤其令人惊讶的是，这个测试在R中的完整实现只有4行代码，可以很好地直观地了解正在发生什么。在第一次构建此测试时，甚至不清楚您是否在做任何真正的数学计算(但您肯定是这样做的)。人们可能很容易指责贝叶斯方法细节的复杂性，但重要的是要认识到，当我们被教导微积分和分析方法的美丽时，我们往往被限制在相对较小的一组问题上，这些问题很好地映射到了计算101的解决方案。当试图用数学方法解决现实世界的问题时，复杂的问题随处可见，分析解决方案要么逃脱，要么失败。