让我们先了解一下基本的定义:马尔可夫链蒙特卡罗(MCMC)方法让我们从一个分布中计算样本，即使我们不能计算它.

这是什么意思?? 让我们回过头来谈谈蒙特卡罗采样.

蒙特卡罗抽样

什么是蒙特卡罗方法?

蒙特卡罗方法, 或蒙特卡罗实验, 是否有一类广泛的计算算法依赖于重复的随机抽样来获得数值结果.” (Wikipedia)

我们来分析一下.

例1:面积估计

Imagine that you have an irregular shape, like the shape presented below:

And you are tasked with determining the area enclosed by this shape. One of the methods you can use is to make small squares in the shape, 数方格数, and that will give you a pretty accurate approximation of the area. 然而，这既困难又耗时.

蒙特卡洛采样来拯救!

First, we draw a big square of a known area around the shape, for example of 50 cm2. Now we “hang” this square and start throwing darts randomly at the shape.

Next, 我们计算矩形中的省道总数和我们感兴趣的形状中的省道数量. 让我们假设使用的“飞镖”总数为100个，其中22个最终落在形状内. 现在可以用简单的公式来计算面积:

形状的面积=正方形的面积*(形状中的部分数)/(正方形中的部分数)

所以，在我们的例子中，这归结为:

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

如果我们把“飞镖”的数量乘以10倍, 这个近似值非常接近真实的答案:

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

这就是我们如何通过使用蒙特卡罗采样来分解复杂的任务，就像上面给出的那样.

大数定律

投掷的飞镖越多，面积近似越接近，这是因为 大数定律:

大数定律是一个定理，它描述了多次进行相同实验的结果. 根据法律, 从大量试验中得到的结果的平均值应该接近期望值, 并且会随着试验的进行而变得更加接近.”

This brings us to our next example, the famous 蒙蒂霍尔问题.

例2:Monty Hall问题

The 蒙蒂霍尔问题 是一个非常著名的脑筋急转弯:

“There are three doors, one has a car behind it, the others have a goat behind them. 你选择一扇门，主人打开另一扇门，告诉你门后面有一只山羊. 然后他问你是否想改变你的决定. Do you? Why? Why not?”

你首先想到的是，无论你换还是不换，获胜的机会都是一样的, 但事实并非如此. 让我们制作一个简单的流程图来演示相同的内容.

假设汽车在3号门后面:

Hence, if you switch, you win ⅔ times, and if you don’t switch, you win only ⅓ times.

我们通过抽样来解决这个问题.

wins = []

For I in range(int(10e6)):

	Car_door = assign_car_door ()

	选择=随机.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.追加(did_he_win)

print(总和(赢得)/ len(赢得))

The assign_car_door () function is just a random number generator that selects a door 0, 1, or 2, 后面有一辆汽车. Using assign_door_to_open 选择一扇门后面有一只山羊，而不是你选择的那扇门，主人打开它. win_or_lose returns true or false, 表示你是否赢得了这辆车, it takes a bool “switch” which says whether you switched the door or not.

让我们运行这个模拟1000万次:

• 如果你不换，赢的概率是0.334134
• 如果你换了，获胜的概率是:0.667255

这与流程图给出的答案非常接近.

In fact, 我们模拟得越多, 答案越接近真实值, 因此证实了大数定律:

从这张表中可以看出:

贝叶斯思维方式

“Frequentist, 被称为更经典的统计学版本, 假设概率是事件的长期频率(因此有了这个标题).”

贝叶斯统计是统计学领域的一种理论，基于贝叶斯对概率的解释，其中概率表达了对事件的信任程度. The degree of belief may be based on prior knowledge about the event, 比如之前的实验结果, 或者个人对事件的看法.” - from 黑客的概率规划和贝叶斯方法

这是什么意思??

In the frequentist way of thinking, we look at probabilities in the long run. 当一个频率学家说有一个0.001% chance of a car crash happening, it means, if we consider infinite car trips, 0.其中0.001%会以崩溃告终.

A Bayesian mindset is different, as we start with a prior, a belief. 如果我们讨论0的信念, it means that your belief is that the event will never happen; conversely, 信念为1意味着你确信它会发生.

然后，一旦我们开始观察数据，我们更新这个信念来考虑数据. 我们该怎么做呢? By using the Bayes theorem.

贝叶斯定理

让我们来分析一下.

• P(A | B) 给出了在给定事件B的情况下事件A的概率. 这是后部, B是我们观察到的数据, so we are essentially saying what the probability of an event happening is, 考虑到我们观察到的数据.
• P(A) 是先验，我们相信事件A会发生.
• P(B | A) 是可能性，我们观察到给定A为真数据的概率是多少.

让我们看一个例子，癌症筛查测试.

癌症的可能性

假设一个病人去做乳房x光检查，结果是阳性的. 病人真的得了癌症的概率是多少?

让我们定义一下概率:

• 1%的女性患有乳腺癌.
• Mammograms detect cancer 80% of the time when it is actually present.
• 9.6%的乳房x光检查错误地报告你患了癌症，而实际上你并没有.

So, 如果你说如果乳房x光检查呈阳性意味着你有80%的几率得了癌症, 这是错误的. You would not take into consideration that having cancer is a rare event, i.e.只有1%的女性患有乳腺癌. We need to take this as prior, and this is where the Bayes theorem comes into play:

P(c +| t +) =(P(t +| c +)*P(c +))/P(t +)

• P(C+ | T+)这是癌症发生的概率, 鉴于测试结果呈阳性, 这是我们感兴趣的.
• P(T+ | C+)这是测试结果为阳性的概率, 考虑到癌症的存在, this, 如上所述，等于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-)这是测试结果为阳性的概率 but there is no cancer, 这是0给出的.096.
• P(C-):这是不患癌症的概率，因为患癌症的概率是1%, 这等于99% = 0.99.
• P(T+ | C+)这是测试结果为阳性的概率, 因为你得了癌症, 这个等于80% = 0.8.
• P(C+): This is the probability of having cancer, 哪个等于1% = 0.01.

把这些代入原来的公式:

So, given that the mammogram came back positive, there is a 7.病人患癌症的几率是76%. 乍一看可能很奇怪，但这是有道理的. 测试结果是假阳性.6%的时间(相当高)，所以在给定的人群中会有很多假阳性. For a rare disease, most of the positive test results will be wrong.

Let us now revisit the 蒙蒂霍尔问题 and solve it using the Bayes theorem.

蒙蒂霍尔问题

先验可以定义为:

• 假设你一开始选择了门A.
• P(H) = ⅓, ⅓. 三扇门各占三分之一, which means that before the door was open and the goat revealed, there was an equal chance of the car being behind any of them.

可能性可以定义为:

• P(D|H), where event D is that Monty chooses door B and there is no car behind it.
• P(D|H) = 1 / 2，如果车在门A后面，因为他有50%的几率选择门B, 50%的几率选择门C.
• P(D|H) = 0，如果汽车在门B后面, since there is a 0% chance that he will choose door B if the car is behind it.
• P(D|H) = 1如果车在C后面，你选择A，他有100%的概率会选择门B.

我们对 P(H|D), 那辆车在一扇门后面的概率是多少如果它显示的是一只山羊在另一扇门后面.

从后面可以看到， P(H|D), that switching to door C from A will increase the probability of winning.

pmmh

Now, let’s combine everything we covered so far and try to understand how pmmh works.

pmmh使用贝叶斯定理得到复分布的后验分布, 从中直接取样是很困难的.

How? Essentially, 它从一个空间中随机选择不同的样本，并检查新样本是否比上一个样本更有可能来自后验, 因为我们在看概率之比, 式(1)中的P(B), 消掉了:

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

The “likelihood” of each new sample is decided by the function f. That’s why f 一定和我们想要采样的后验成正比吗.

决定θ '是被接受还是被拒绝, the following ratio must be computed for each new proposed θ’:

其中θ为旧样本， P(D| θ) 样本的似然是θ.

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

X ~ N(mean, std)

现在我们需要定义均值和标准差的先验. 为简单起见，我们假设两者都服从均值为1，标准差为2的正态分布:

均值~ N(1,2)
std ~ N(1,2)

Now, let’s generate some data and assume the true mean and std are 0.5 and 1.2,分别.

导入numpy为np

进口matplotlib.Pyplot为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 to model the data and find the posterior distribution for the mean and std.

将pymc3导入为PM

with pm.将()作为模型:

  mean = pm.正态(“平均值”，mu = 1, sigma = 2)

  std = pm.正态分布("std"， mu = 1, sigma = 2)

  x = pm.正态(“X”，mu =平均值，sigma =标准差，observed = X)

  step = pm.Metropolis()

  trace = pm.样品(5000，步长=步长)

让我们复习一下代码.

首先，我们定义均值和std的先验i.e.，平均值为1,STD为2.

x = pm.正态(“X”，mu =平均值，sigma =标准差，observed = X)

In this line, we define the variable that we are interested in; it takes the mean and std we defined earlier, 它也取我们之前计算过的观测值.

让我们看看结果:

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

_ = plt.Hist (跟踪(“的意思是”)， bins = 50, label = "mean")

_ = plt.legend()

std以1为中心.2，均值是1.55——非常接近!

Let’s now do it from scratch to understand pmmh.

From Scratch

生成提案

首先，让我们了解大都会黑斯廷斯是做什么的. 在这篇文章的前面，我们说过， “pmmh randomly selects different samples from a space,” 但是它怎么知道选择哪个样本呢?

它使用提案分发来实现这一点. 它是一个正态分布，以当前接受的样本为中心，STD为0.5. Where 0.5是一个超参数, we can tweak; the more the STD of the proposal, the further it will search from the currently accepted sample. 我们来编码一下.

因为我们需要找到分布的均值和标准差, this function will take the currently accepted mean and the currently accepted std, 并对两者都提出建议.

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

    return np.random.Normal (mean_current, proposal_width)， \

            np.random.正常(std_current proposal_width)

接受/拒绝提议

Now, let’s code the logic that accepts or rejects the proposal. 它有两个部分: prior and the likelihood.

First, let’s calculate the probability of the proposal coming from the prior:

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

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

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

它只需要 mean and std 然后计算这个的可能性 mean and std 来自 先验分布 我们定义的.

在计算可能性时, 我们计算我们看到的数据有多大可能来自提议的分布.

Def likelihood(mean, std, data):

    return np.prod(st.规范(意思是,std).pdf(X))

对于每个数据点, we multiply the probability of that data point coming from the proposed distribution.

现在，我们需要调用这些函数来获取当前的均值和std以及建议的 mean and std.

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.规范(prior_mean [0], prior_mean [1]).pdf(mean)* \

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

    

    Def likelihood(mean, std, data):

        return np.prod(st.规范(意思是,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)

    

    返回(prior_proposal * likehood_proposal) / (prior_current * likehood_current)

Now, 让我们来创建最后一个函数，它将把所有东西联系在一起，并生成我们需要的后验样本.

产生后路

Now, we call the functions we have written above and generate the posterior!

Def get_trace(data, samples = 5000):

    Mean_prior = 1

    Std_prior = 2

    

    Mean_current = mean_prior

    Std_current = std_prior

    

    trace = {

        “的意思是”:[mean_current],

        “性病”(std_current):

    }

    

    对于tqdm(range(samples))中的I:

        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() 

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.规范(prior_mean [0], prior_mean [1]).logpdf(平均)+ \

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

    

    Def likelihood(mean, std, data):

        return np.sum(st.规范(意思是,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)

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

_ = plt.Hist (跟踪(“的意思是”)， bins = 50, label = "mean")

_ = plt.legend()