作者都是各自领域经过审查的专家,并撰写他们有经验的主题. 我们所有的内容都经过同行评审,并由同一领域的Toptal专家验证.
Divyanshu卡尔拉的头像

By Divyanshu卡尔拉

Divyanshu是一名精通多种语言的数据科学家和全栈开发人员. 他在这个领域发表了三篇研究论文.

专业知识

以前在

公司
分享

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

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

蒙特卡罗抽样

什么是蒙特卡罗方法?

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

我们来分析一下.

例1:面积估计

想象你有一个不规则的形状,就像下面所示的形状:

面积估计

你的任务是确定这个形状所包围的面积. 你可以使用的一种方法是在形状上做小正方形, 数方格数, 这样就可以得到一个非常精确的面积近似值. 然而,这既困难又耗时.

蒙特卡洛采样来拯救!

首先,我们在形状周围画一个已知面积的大正方形,例如50cm2. 现在我们“悬挂”这个正方形,并开始向这个形状随机投掷飞镖.

面积估计

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

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

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

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

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

面积估计

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

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

大数定律

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

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

这就引出了下一个例子,著名的蒙蒂·霍尔问题.

例2:Monty Hall问题

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

蒙蒂·霍尔问题

“有三扇门,一扇后面有一辆汽车,其他的后面有一只山羊. 你选择一扇门,主人打开另一扇门,告诉你门后面有一只山羊. 然后他问你是否想改变你的决定. 你? 为什么? 为什么不?”

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

假设汽车在3号门后面:

蒙蒂·霍尔问题

因此,如果你换,你赢三分之二,如果你不换,你只赢三分之一.

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


胜利= []

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

	Car_door = assign_car_door ()

	选择=随机.r和int (0, 2)

	Opened_door = assign_door_to_open(car_door)

	did_he_win = win_or_lose(选择,car_door, opened_door,开关= False)

	赢了.追加(did_he_win)

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

assign_car_door () 函数只是一个随机数生成器,它选择一个门为0, 1, or 2, 后面有一辆汽车. 使用 assign_door_to_open 选择一扇门后面有一只山羊,而不是你选择的那扇门,主人打开它. win_or_lose 返回 真正的 or , 表示你是否赢得了这辆车, 它有一个bool“switch”,表示你是否开关了门.

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

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

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

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

大数定律

从这张表中可以看出:

模拟运行换手赢的概率如果你不换,获胜的概率
100.60.2
10^20.630.33
10^30.6610.333
10^40.66830.3236
10^50.667620.33171
10^60.6672550.334134
10^70.66687560.3332821

贝叶斯思维方式

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

贝叶斯统计是统计学领域的一种理论,基于贝叶斯对概率的解释,其中概率表达了对事件的信任程度. 相信的程度可能基于对事件的先验知识, 比如之前的实验结果, 或者个人对事件的看法.” ——从 黑客的概率规划和贝叶斯方法

这是什么意思??

在频率主义者的思维方式中,我们着眼于长期的概率. 当一个频率学家说有一个0.发生车祸的概率是0.001%,这意味着,如果我们考虑无限次的汽车旅行.其中0.001%会以崩溃告终.

贝叶斯思维是不同的,因为我们从一个先验,一个信念开始. 如果我们讨论0的信念, it 的意思是s that your belief is that the event will never happen; conversely, 信念为1意味着你确信它会发生.

然后,一旦我们开始观察数据,我们更新这个信念来考虑数据. 我们该怎么做呢? 通过使用 贝叶斯定理.

贝叶斯定理

贝叶斯定理

让我们来分析一下.

  • P(a | b) 给出了在给定事件B的情况下事件A的概率. 这是后部, B是我们观察到的数据, 我们本质上说的是一个事件发生的概率, 考虑到我们观察到的数据.
  • P(A) 是先验,我们相信事件A会发生.
  • P(b | a) 是可能性,我们观察到给定A为真数据的概率是多少.

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

癌症的可能性

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

让我们定义一下概率:

  • 1%的女性患有乳腺癌.
  • 乳房x光检查发现癌症的几率是80%.
  • 9.6%的乳房x光检查错误地报告你患了癌症,而实际上你并没有.

So, 如果你说如果乳房x光检查呈阳性意味着你有80%的几率得了癌症, 这是错误的. 你没有考虑到患癌症是一件罕见的事.e.只有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%. 乍一看可能很奇怪,但这是有道理的. 测试结果是假阳性.6%的时间(相当高),所以在给定的人群中会有很多假阳性. 对于一种罕见的疾病,大多数阳性的检测结果都会是错误的.

现在让我们回顾一下Monty Hall问题并使用贝叶斯定理来解决它.

蒙蒂霍尔问题

先验可以定义为:

  • 假设你一开始选择了门A.
  • P(h) = 1 / 3, 1 / 3. 三扇门各占三分之一, 也就是说在门打开,山羊出现之前, 汽车落在他们后面的几率是一样的.

可能性可以定义为:

  • P (D | H),事件D是Monty选择了B门,门后没有车.
  • P (D | H) = 1 / 2,如果车在门A后面,因为他有50%的几率选择门B, 50%的几率选择门C.
  • P (D | H) = 0,如果汽车在门B后面, 因为如果车在门B后面,他选择门B的几率是0%.
  • P (D | H) = 1如果车在C后面,你选择A,他有100%的概率会选择门B.

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

贝叶斯定理

从后面可以看到, P (H | D)即从A转到C会增加获胜的概率.

pmmh

现在,让我们结合到目前为止所讨论的内容,并尝试理解其中的原理 pmmh 作品.

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

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

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

每个新样本的“似然”由函数决定 f. 这就是为什么 f 一定和我们想要采样的后验成正比吗.

决定θ '是被接受还是被拒绝, 对于每个新提出的θ′,必须计算以下比率:

pmmh

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

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

X ~ N(的意思是, 性病)

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

均值~ N(1,2)
性病 ~ N(1,2)

现在,让我们生成一些数据并假设真实均值和性病为0.5和1.2,分别.


导入numpy为np

进口matplotlib.Pyplot为PLT

进口scipy.统计数据如下:

的意思是X = 1.5

性病X = 1.2

X = np.随机.normal(的意思是X, 性病X, size = 1000)

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

PyMC3

让我们首先使用一个库 PyMC3 对数据进行建模,找出均值和STD的后验分布.


将pymc3导入为PM

与点.将()作为模型:

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

  STD = PM.正态分布("性病", mu = 1, sigma = 2)

  X = PM.正态(“X”,mu =平均值,sigma =标准差,observed = X)

  步骤= PM.大都市()

  Trace = PM.样品(5000,步长=步长)

让我们复习一下代码.

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

X = PM.正态(“X”,mu =平均值,sigma =标准差,observed = X)

在这一行, we define the variable that we are interested in; it takes the 的意思是 和 性病 we defined earlier, 它也取我们之前计算过的观测值.

让我们看看结果:


_ = PLT.Hist (跟踪(“标准”), bins = 50, label = "性病")

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

_ = PLT.传奇()

性病以1为中心.2,均值是1.55——非常接近!

现在让我们从头开始理解大都会黑斯廷斯案.

从头开始

生成提案

首先,让我们了解大都会黑斯廷斯是做什么的. 在这篇文章的前面,我们说过, “大都会-黑斯廷斯从一个空间随机选择不同的样本,” 但是它怎么知道选择哪个样本呢?

它使用提案分发来实现这一点. 它是一个正态分布,以当前接受的样本为中心,STD为0.5. 在0.5是一个超参数, we can tweak; the more the STD of the proposal, 它将从当前接受的样本中搜索得越远. 我们来编码一下.

因为我们需要找到分布的均值和标准差, 这个函数将取当前接受的平均值和当前接受的STD, 并对两者都提出建议.


Def get_proposal(的意思是_current, 性病_current, proposal_width = 0.5):

    返回np.随机.Normal (的意思是_current, proposal_width), \

            np.随机.正常(性病_current proposal_width)

接受/拒绝提议

现在,让我们编写接受或拒绝提议的逻辑. 它有两个部分: 之前可能性.

首先,让我们计算来自先验的提议的概率:


Def 之前(的意思是, 性病, 之前_的意思是, 之前_性病):

        返回圣.规范(之前_的意思是 [0], 之前_的意思是 [1]).pdf(意味着)* \

                st.规范(之前_性病 [0], 之前_性病 [1]).pdf (性病)

它只需要 的意思是性病 然后计算这个的可能性 的意思是性病 来自 先验分布 我们定义的.

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


Def 可能性(的意思是, 性病, data):

    返回np.刺激(圣.规范(意思是,性病).pdf (X))

对于每个数据点, 我们乘以来自提议分布的数据点的概率.

现在,我们需要调用这些函数来获取当前的均值和性病以及建议的 的意思是性病.


Prior_current = 之前(的意思是_current, 性病_current, 之前_的意思是, 之前_性病)

可能性(的意思是_current, 性病_current, data)

    

Prior_proposal = 之前(的意思是_proposal, 性病_proposal, 之前_的意思是, 之前_性病)

可能性(的意思是_proposal, 性病_proposal, data)

整个代码:


Def accept_proposal(的意思是_proposal, 性病_proposal, 的意思是_current, \

                    Std_current, 之前_的意思是, 之前_性病, data):

    

    Def 之前(的意思是, 性病, 之前_的意思是, 之前_性病):

        返回圣.规范(之前_的意思是 [0], 之前_的意思是 [1]).pdf(意味着)* \

                st.规范(之前_性病 [0], 之前_性病 [1]).pdf (性病)

    

    Def 可能性(的意思是, 性病, data):

        返回np.刺激(圣.规范(意思是,性病).pdf (X))

    Prior_current = 之前(的意思是_current, 性病_current, 之前_的意思是, 之前_性病)

    可能性(的意思是_current, 性病_current, data)

    

    Prior_proposal = 之前(的意思是_proposal, 性病_proposal, 之前_的意思是, 之前_性病)

    可能性(的意思是_proposal, 性病_proposal, data)

    

    返回(之前_proposal * likehood_proposal) / (之前_current * likehood_current)

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

产生后路

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


Def get_trace(data, samples = 5000):

    Mean_之前 = 1

    Std_之前 = 2

    

    Mean_current = 的意思是_之前

    Std_current = 性病_之前

    

    Trace = {

        “的意思是”:[的意思是_current],

        “性病”(性病_current):

    }

    

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

        Mean_proposal, 性病_proposal = get_proposal(的意思是_current, 性病_current)

        Acceptance_prob = accept_proposal(的意思是_proposal, 性病_proposal, 的意思是_current, \

                                         Std_current, [的意思是_之前, 性病_之前], \

                                          [的意思是_之前, 性病_之前], data)

        如果np.随机.r和() 

改进空间

日志是你的朋友!

回忆 A * b = antilog(log(A) + log(b))A / b = antilog(log(A) - log(b)).

这对我们很有利因为我们要处理的概率很小, 乘以它会得到零, 所以我们要加上log问题, 问题解决了!

因此,我们之前的代码变成:


Def accept_proposal(的意思是_proposal, 性病_proposal, 的意思是_current, \

                    Std_current, 之前_的意思是, 之前_性病, data):

    

    Def 之前(的意思是, 性病, 之前_的意思是, 之前_性病):

        返回圣.规范(之前_的意思是 [0], 之前_的意思是 [1]).logpdf(平均)+ \

                st.规范(之前_性病 [0], 之前_性病 [1]).logpdf (性病)

    

    Def 可能性(的意思是, 性病, data):

        返回np.sum(圣.规范(意思是,性病).logpdf (X))

    Prior_current = 之前(的意思是_current, 性病_current, 之前_的意思是, 之前_性病)

    可能性(的意思是_current, 性病_current, data)

    

    Prior_proposal = 之前(的意思是_proposal, 性病_proposal, 之前_的意思是, 之前_性病)

    可能性(的意思是_proposal, 性病_proposal, data)

    

    返回(之前_proposal + likehood_proposal) - (之前_current + likehood_current)

由于我们返回的是接受概率的对数:

如果np.随机.r和()

就变成:

如果数学.日志(np.随机.r和())

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


_ = PLT.Hist (跟踪(“标准”), bins = 50, label = "性病")

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

_ = PLT.传奇()

如你所见, 性病 以1为中心.2、 的意思是 at 1.5. 我们做到了!

前进

到现在为止, 你可能知道大都会黑斯廷斯是如何工作的,你可能想知道你可以在哪里使用它.

好吧, 黑客的贝叶斯方法 是一本解释概率编程的好书吗, 如果你想了解更多关于贝叶斯定理及其应用的知识, 认为贝叶斯 是艾伦写的一本很棒的书. 唐尼.

谢谢你的阅读, 我希望这篇文章能鼓励你去发现贝叶斯统计的神奇世界.

了解基本知识

  • MCMC代表什么?

    MCMC代表马尔科夫链蒙特卡罗, 哪一类是从概率分布中抽样的方法.

  • 什么是贝叶斯决策方法?

    在贝叶斯决策方法中, 首先从先验开始, 这就是你的信念, 然后当数据进来时, 你结合这些数据来更新这些先验,从而得到后验.

  • 什么是贝叶斯模型?

    贝叶斯模型是一种统计模型,使用概率来表示模型中的所有不确定性.

就这一主题咨询作者或专家.
预约电话
Divyanshu卡尔拉的头像
Divyanshu卡尔拉

位于 新德里,印度

成员自 二零一七年十月三十日

作者简介

Divyanshu是一名精通多种语言的数据科学家和全栈开发人员. 他在这个领域发表了三篇研究论文.

Toptal作者都是各自领域经过审查的专家,并撰写他们有经验的主题. 我们所有的内容都经过同行评审,并由同一领域的Toptal专家验证.

专业知识

以前在

公司

世界级的文章,每周发一次.

订阅意味着同意我们的 隐私政策

世界级的文章,每周发一次.

订阅意味着同意我们的 隐私政策

Toptal开发者

加入总冠军® 社区.