MCMC方法：Metropolis-Hastings和贝叶斯推理

本文概述

• 蒙特卡洛采样
• 贝叶斯思维方式
• Metropolis-Hastings
• 从头开始
• 总结

让我们排除掉基本定义：马尔可夫链蒙特卡洛(MCMC)方法使我们能够从分布中计算样本, 即使我们无法计算它。

这是什么意思？让我们备份一下, 谈论蒙特卡洛采样。

蒙特卡洛采样

什么是蒙特卡洛方法？

“蒙特卡洛方法或蒙特卡洛实验是一类广泛的计算算法, 它们依赖于重复随机采样来获得数值结果。” (维基百科)

让我们分解一下。

示例1：面积估计

想象一下, 你有一个不规则的形状, 如下图所示：

你的任务是确定此形状所包围的区域。可以使用的方法之一是制作形状的小正方形, 计算正方形, 这将为你提供非常准确的面积近似值。然而, 这是困难且耗时的。

蒙特卡洛采样救援！

首先, 我们在形状周围绘制一个已知区域的大正方形, 例如50平方厘米。现在我们”悬挂”这个正方形, 开始向形状随机投掷飞镖。

接下来, 我们计算矩形中的飞镖总数和感兴趣的形状中的飞镖数量。假设使用的”飞镖”总数为100, 并且其中22枚最终落入该形状中。现在可以通过简单的公式来计算面积：

形状面积=正方形面积*(形状中的飞镖数)/(正方形中的飞镖数)

因此, 在我们的案例中, 这归结为：

形状面积= 50 * 22/100 = 11 cm2

如果我们将”飞镖”的数量乘以10倍, 则此近似值将非常接近真实答案：

形状面积= 50 * 280/1000 = 14 cm2

这就是我们通过使用蒙特卡洛采样分解复杂任务的方式, 如上面给出的那样。

大数定律

我们投掷的飞镖越多, 面积近似值就越接近, 这是由于大数定律：

“大数定律是一个定理, 描述了多次执行相同实验的结果。根据法律, 从大量试验中获得的结果平均值应接近预期值, 并且随着进行更多试验而趋于接近。”

这将我们带入下一个示例, 即著名的Monty Hall问题。

示例2：蒙蒂·霍尔问题

Monty Hall问题是一个非常著名的脑筋急转弯：

“有三扇门, 一扇门后面有辆汽车, 另一扇门后面有只山羊。你选择一扇门, 主持人会打开另一扇门, 并向你显示后面有只山羊。然后, 他问你是否要更改决定。你做？为什么？为什么不？”

首先想到的是, 无论你是否切换, 获胜的机会均等, 但这不是事实。让我们制作一个简单的流程图来说明这一点。

假设汽车在3号门后面：

因此, 如果你切换, 你将赢得⅔次；如果你不切换, 你将仅赢得win次。

让我们通过采样来解决这个问题。

wins = []

for i in range(int(10e6)):

	car_door = assign_car_door()

	choice = random.randint(0, 2)

	opened_door = assign_door_to_open(car_door)

	did_he_win = win_or_lose(choice, car_door, opened_door, switch = False)

	wins.append(did_he_win)

print(sum(wins)/len(wins))

Assign_car_door()函数只是一个随机数生成器, 它选择0、1或2号门, 后面是汽车。使用assign_door_to_open可以选择一扇门, 该门后面是山羊, 而不是你选择的门, 主机会打开它。 win_or_lose返回true或false, 表示你是否赢得了汽车, 它需要一个bool” switch”(开关)来表示是否切换门。

让我们运行此模拟一千万次：

• 如果你不切换, 则获胜的概率：0.334134
• 如果你进行切换, 则获胜的概率：0.667255

这与流程图给我们的答案非常接近。

实际上, 我们进行此模拟的次数越多, 答案就越接近真实值, 从而验证了大数定律：

x轴是运行的模拟次数, y轴是不切换时获胜的概率。

从此表可以看到相同的结果：

贝叶斯思维方式

“频率论者, 被称为统计的经典形式, 认为概率是事件的长期发生频率(因此被赋予标题)。”

贝叶斯统计是一种基于贝叶斯概率解释的统计学领域的理论, 其中概率表示对事件的信任程度。信任程度可以基于对事件的先前了解, 例如先前实验的结果, 或者对事件的个人信念。” -来自黑客的概率编程和贝叶斯方法

这是什么意思？

以常客主义的思维方式, 我们从长远来看概率。当常客说发生车祸的可能性为0.001％时, 这意味着, 如果我们考虑无限次出车, 那么其中0.001％的事故将以撞车事故告终。

贝叶斯的心态是不同的, 因为我们从先验的信念开始。如果我们将信念置为0, 则表示你的信念是该事件永远不会发生。相反, 相信1表示你确定会发生。

然后, 一旦我们开始观察数据, 便会更新此信念以考虑数据。我们如何做到这一点？通过使用贝叶斯定理。

贝叶斯定理

….(1)

让我们分解一下。

• P(A | B)给了我们事件B给定事件B的概率。这是后验, B是我们观察到的数据, 因此考虑到我们观察到的数据, 我们实质上是在说事件发生的概率是多少。
• P(A)是先验的, 我们相信事件A将会发生。
• P(B | A)是可能性, 假设A为真, 我们将观察数据的概率是多少。

让我们看一个例子, 癌症筛查测试。

癌症的可能性

假设患者去做乳房X线检查, 然后乳房X线检查恢复为阳性。患者实际患癌症的概率是多少？

让我们定义概率：

• 1％的女性患有乳腺癌。
• 乳房X光检查可以在实际存在的癌症中80％地检测出癌症。
• 9.6％的乳房X线照片错误地报告你实际上没有癌症。

因此, 如果你要说的是, 如果乳房X光检查结果呈阳性, 这意味着你有80％的机会患上癌症, 那是错误的。你不会考虑到癌症是罕见的事件, 即只有1％的女性患有乳腺癌。我们需要将此作为先验, 这就是贝叶斯定理起作用的地方：

P(C + | T +)=(P(T + | C +)* P(C +))/ P(T +)

• P(C + | T +)：这是癌症存在的可能性, 因为测试结果呈阳性, 这是我们感兴趣的。
• P(T + | C +)：鉴于存在癌症, 这是测试呈阳性的概率, 如上所述, 这等于80％= 0.8。
• P(C +)：这是先验概率, 即一个人患有癌症的概率, 等于1％= 0.01。
• P(T +)：这是无论什么情况测试均呈阳性的概率, 因此它具有两个成分：P(T +)= P(T + | C-)P(C-)+ P(T + | C +)P (C +)
• P(T + | C-)：这是测试返回阳性但没有癌症的概率, 由0.096给出。
• P(C-)：这是没有患癌症的可能性, 因为患癌症的可能性是1％, 等于99％= 0.99。
• P(T + | C +)：考虑到你患有癌症, 这是测试返回阳性的概率, 等于80％= 0.8。
• P(C +)：这是罹患癌症的概率, 等于1％= 0.01。

将所有这些插入原始公式中：

因此, 考虑到乳房X线检查结果呈阳性, 该患者患癌症的几率为7.76％。乍一看似乎很奇怪, 但是很有意义。该测试给出9.6％的时间出现假阳性(相当高), 因此在给定的总体中将有很多假阳性。对于罕见疾病, 大多数阳性测试结果将是错误的。

现在让我们重新审视Monty Hall问题, 并使用贝叶斯定理解决它。

蒙蒂·霍尔问题

先验可以定义为：

• 假设你一开始选择门A作为选择。
• P(H)=⅓, ⅓ ⅓对于所有三个门, 这意味着在门打开和山羊露出来之前, 汽车有相等的机会被挡在其中的任何一个后面。

可能性可以定义为：

• P(D | H), 其中事件D是Monty选择门B, 且门后没有汽车。
• 如果汽车在门A的后面, 则P(D | H)=½-因为他有50％的机会选择门B, 而他有50％的机会选择门C。
• 如果汽车在门B的后面, 则P(D | H)= 0, 因为如果汽车在门B的后面, 他有0％的机会选择门B。
• 如果汽车落后于C, 而你选择A, 则P(D | H)= 1, 则他有100％的可能性选择门B。

我们对P(H | D)感兴趣, P(H | D)是汽车在门后面的概率, 因为它向我们展示了另一扇门后面的山羊。

从后验P(H | D)可以看出, 从A切换到门C将增加获胜的可能性。

Metropolis-Hastings

现在, 让我们结合目前为止涵盖的所有内容, 并尝试了解Metropolis-Hastings的工作原理。

Metropolis-Hastings使用贝叶斯定理来获得复杂分布的后验分布, 从中很难直接采样。

怎么样？从本质上讲, 它是从一个空间中随机选择不同的样本, 并检查新样本是否比后一个样本更可能来自后验样本, 因为我们正在查看方程(1)中的概率比P(B), 取消：

P(接受)=(P(newSample)*新样本的可能性)/(P(oldSample)*旧样本的可能性)

每个新样本的”可能性”由函数f决定。这就是为什么f必须与我们要从中采样的后验成比例的原因。

为了确定是接受还是拒绝θ′, 必须为每个新提议的θ′计算以下比率：

在θ是旧样本的情况下, P(D |θ)是样本θ的可能性。

让我们用一个例子来更好地理解这一点。假设你有数据, 但是想找出适合它的正态分布, 所以：

X〜N(平均值, 标准)

现在我们需要定义均值和标准差的先验值。为简单起见, 我们将假定两者均遵循均值1和std 2的正态分布：

均值〜N(1, 2)

std〜N(1, 2)

现在, 我们生成一些数据, 并假设真实均值和标准差分别为0.5和1.2。

import numpy as np

import matplotlib.pyplot as plt

import scipy.stats as st

meanX = 1.5

stdX = 1.2

X = np.random.normal(meanX, stdX, size = 1000)

_ = plt.hist(X, bins = 50)

PyMC3

首先, 我们使用一个名为PyMC3的库对数据进行建模, 并找到均值和标准差的后验分布。

import pymc3 as pm

with pm.Model() as model:

  mean = pm.Normal("mean", mu = 1, sigma = 2)

  std = pm.Normal("std", mu = 1, sigma = 2)

  x = pm.Normal("X", mu = mean, sigma = std, observed = X)

  step = pm.Metropolis()

  trace = pm.sample(5000, step = step)

让我们来看一下代码。

首先, 我们定义均值和std的先验值, 即均值1和std 2的法线。

x = pm.Normal(” X”, mu =平均值, sigma = std, 观察到的= X)

在这一行中, 我们定义了我们感兴趣的变量；它采用了我们先前定义的均值和std, 还采用了我们先前计算的观测值。

我们来看一下结果：

_ = plt.hist(trace['std'], bins = 50, label = "std")

_ = plt.hist(trace['mean'], bins = 50, label = "mean")

_ = plt.legend()

性病中心位于1.2附近, 均值位于1.55-非常接近！

现在, 让我们从头开始, 了解Metropolis-Hastings。

从头开始

产生提案

首先, 让我们了解一下Metropolis-Hastings的作用。在本文的前面, 我们说过” Metropolis-Hastings从一个空间中随机选择不同的样本”, 但是如何知道要选择哪个样本呢？

这样做是使用提案分配。它是一个以当前接受的样本为中心的正态分布, 其STD为0.5。这里0.5是一个超参数, 我们可以进行调整；提案的STD越多, 它将从当前接受的样本中搜索得越远。让我们对此进行编码。

由于我们必须找到分布的均值和标准差, 因此此函数将采用当前接受的均值和当前接受的标准差, 并返回两者的建议。

def get_proposal(mean_current, std_current, proposal_width = 0.5):

    return np.random.normal(mean_current, proposal_width), \

            np.random.normal(std_current, proposal_width)

接受/拒绝提案

现在, 让我们对接受或拒绝提案的逻辑进行编码。它包括两个部分：先验和可能性。

首先, 让我们计算提案来自先前提案的可能性：

def prior(mean, std, prior_mean, prior_std):

        return st.norm(prior_mean[0], prior_mean[1]).pdf(mean)* \

                st.norm(prior_std[0], prior_std[1]).pdf(std)

它仅取均值和标准差, 并计算该均值和标准差来自我们定义的先前分布的可能性。

在计算可能性时, 我们计算出我们看到的数据来自建议分布的可能性。

def likelihood(mean, std, data):

    return np.prod(st.norm(mean, std).pdf(X))

因此, 对于每个数据点, 我们乘以该数据点来自建议分布的概率。

现在, 我们需要为当前均值和标准差以及建议的均值和标准差调用这些函数。

prior_current = prior(mean_current, std_current, prior_mean, prior_std)

likelihood_current = likelihood(mean_current, std_current, data)

    

prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std)

likelihood_proposal = likelihood(mean_proposal, std_proposal, data)

整个代码：

def accept_proposal(mean_proposal, std_proposal, mean_current, \

                    std_current, prior_mean, prior_std, data):

    

    def prior(mean, std, prior_mean, prior_std):

        return st.norm(prior_mean[0], prior_mean[1]).pdf(mean)* \

                st.norm(prior_std[0], prior_std[1]).pdf(std)

    

    def likelihood(mean, std, data):

        return np.prod(st.norm(mean, std).pdf(X))

    prior_current = prior(mean_current, std_current, prior_mean, prior_std)

    likelihood_current = likelihood(mean_current, std_current, data)

    

    prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std)

    likelihood_proposal = likelihood(mean_proposal, std_proposal, data)

    

    return (prior_proposal * likelihood_proposal) / (prior_current * likelihood_current)

现在, 让我们创建最终函数, 它将所有内容绑定在一起并生成我们所需的后验样本。

产生后验

现在, 我们调用上面编写的函数并生成后验！

def get_trace(data, samples = 5000):

    mean_prior = 1

    std_prior = 2

    

    mean_current = mean_prior

    std_current = std_prior

    

    trace = {

        "mean": [mean_current], "std": [std_current]

    }

    

    for i in tqdm(range(samples)):

        mean_proposal, std_proposal = get_proposal(mean_current, std_current)

        acceptance_prob = accept_proposal(mean_proposal, std_proposal, mean_current, \

                                         std_current, [mean_prior, std_prior], \

                                          [mean_prior, std_prior], data)

        if np.random.rand() < acceptance_prob:

            mean_current = mean_proposal

            std_current = std_proposal

        trace['mean'].append(mean_current)

        trace['std'].append(std_current)

    return trace

改进空间

日志是你的朋友！

回想a * b = antilog(log(a)+ log(b))和a / b = antilog(log(a)-log(b))。

这对我们是有益的, 因为我们将处理很小的概率, 相乘将得出零, 因此我们宁愿添加对数概率, 解决问题！

因此, 我们之前的代码变为：

def accept_proposal(mean_proposal, std_proposal, mean_current, \

                    std_current, prior_mean, prior_std, data):

    

    def prior(mean, std, prior_mean, prior_std):

        return st.norm(prior_mean[0], prior_mean[1]).logpdf(mean)+ \

                st.norm(prior_std[0], prior_std[1]).logpdf(std)

    

    def likelihood(mean, std, data):

        return np.sum(st.norm(mean, std).logpdf(X))

    prior_current = prior(mean_current, std_current, prior_mean, prior_std)

    likelihood_current = likelihood(mean_current, std_current, data)

    

    prior_proposal = prior(mean_proposal, std_proposal, prior_mean, prior_std)

    likelihood_proposal = likelihood(mean_proposal, std_proposal, data)

    

    return (prior_proposal + likelihood_proposal) - (prior_current + likelihood_current)

由于我们正在返回接受概率的对数：

如果np.random.rand()<accept_prob：

成为：

如果math.log(np.random.rand())<accept_prob：

让我们运行新代码并绘制结果。

_ = plt.hist(trace['std'], bins = 50, label = "std")

_ = plt.hist(trace['mean'], bins = 50, label = "mean")

_ = plt.legend()

如你所见, std的中心为1.2, 平均值为1.5。我们做到了！

总结

到目前为止, 你希望了解Metropolis-Hastings的工作原理, 并且可能想知道可以在哪里使用它。

好吧, 《贝叶斯黑客方法》是一本出色的书, 解释了概率编程, 如果你想了解有关贝叶斯定理及其应用的更多信息, 《 Think Bayes》是Allen B. Downey的一本好书。

感谢你的阅读, 我希望本文鼓励你发现贝叶斯统计数据的惊人世界。