蒙特卡罗计算积分

点击上方“小白学视觉”,选择加"星标"或“置顶

重磅干货,第一时间送达

通常情况下,我们不能解析地求解积分,必须借助其他方法,其中就包括蒙特卡罗积分。你可能还记得,函数的积分可以解释为函数曲线下的面积。

蒙特卡罗积分的工作原理是在a和b之间的不同随机点计算一个函数,将矩形的面积相加,取和的平均值。随着点数的增加,所得结果接近于积分的实际解。

cb9ee3321974816a77116d50e09d13bc.png

蒙特卡罗积分用代数表示:

c49e8f08b7764d99c543c5f9e96a7a8e.png

与其他数值方法相比,蒙特卡罗积分特别适合于计算奇数形状的面积。

769d570be451a002bfbb892e9b96a49c.png

在上一节中,我们看到如何使用蒙特卡罗积分来确定后验概率,当我们知道先验和似然,但缺少规范化常数。

贝叶斯统计

后验概率是指贝叶斯公式中的一个特定项。

7e97a58c166dd637103498807502f82e.png

贝叶斯定理可以用来计算一个人在某一特定疾病的筛查测试中呈阳性的人实际上患有该病的概率。

60513015f068fab727780b3acd8f6a58.png

如果我们已经知道P(A),P(B)和P(B | A),但想知道P(A | B),我们就用这个公式。例如,假设我们正在检测一种感染1%人口的罕见疾病。医学专业人员已经开发出一种高度敏感和特异的检测方法,但还不够完善。

  • 99%的病人检测呈阳性

  • 99%的健康患者检测为阴性

贝叶斯定理告诉我们:

606ee0566a3f1c9485cc77c35d3b9467.png

假设我们有10000人,100人生病,9900人健康。此外,在给他们所有的测试后,我们会让99个生病的人测试生病,但是99个健康的人也测试生病。因此,我们将得到以下概率。

「p(sick) = 0.01」

「p(not sick) = 1–0.01 = 0.99」

「p(+|sick) = 0.99」

「p(+|not sick) = 0.01」

4df481f4fdddada7523ba7708f827638.png

Bayes定理在概率分布中的应用

在前面的例子中,我们假设一个人患病的先验概率是一个已知的量,精确到0.001。

然而,在现实世界中,认为0.001的概率事实上如此精确是不合理的。一个给定的人患病的概率会因其年龄、性别、体重等而有很大差异。一般来说,我们对给定先验概率的认识还远远不够完善,因为它是从以前的样本中得出的(这意味着不同的人群可能会对先验概率给出不同的估计)。

在贝叶斯统计中,我们可以用先验概率的分布来代替这个0.001的值,这个分布捕捉了我们关于其真实值的先验不确定性。包含先验概率分布最终产生的后验概率也不再是单一数量;相反,后验概率也变成了概率分布。这与传统的观点相反,后者假设参数是固定的量。

归一化常数

正如我们在Gibbs抽样和Metropolis-Hasting的文章中看到的,蒙特卡洛方法可以用来计算归一化常数未知时的后验概率分布。

让我们来探究一下为什么我们首先需要一个标准化常数。在概率论中,规范化常数是一个函数必须乘以的常数,因此它的图下面积为1。还是不清楚?让我们看一个例子。

回想一下正态分布的函数可以写成:

24cada4078bac0dcd4876ea72601b1c6.png

2*pi的平方根是归一化常数。

让我们来看看我们是如何确定它的。我们从以下函数开始(假设均值为0,方差为1):

e8dcaf686b75f361ddd45fca6d8c5489.png

如果我们能画出一个曲线的话。

cb15a15ca89d19c000d98e53789a3b5e.png

问题在于,如果我们取曲线下的面积,它不等于1,这要求它是一个概率密度函数。因此,我们将函数除以积分的结果(归一化常数)。

a82f6bf6799912dfdbd24c68e5ea50d5.png

回到手头的问题,即如何在没有归一化常数的情况下计算后验概率……事实证明,对于连续样本空间,规范化常数可以重写为:

4b416e949c085c60f6800fa2a77f4e6d.png

在这一点上,你应该考虑蒙特卡罗积分!

Python代码

让我们看看如何通过在Python中执行蒙特卡洛积分来确定后验概率。我们从导入所需的库开始,并设置随机种子以确保结果是可重复的。

import os
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as st

np.random.seed(42)

然后我们设置β分布和二项分布的参数值。

a, b = 10, 10
n = 100
h = 59
thetas = np.linspace(0, 1, 200)

概率密度函数的范围从0到1。因此,我们可以简化方程。

baf3e80a67dde278c5f2404830b71025.png

在代码中,前面的等式写如下:

prior = st.beta(a, b).pdf(thetas)
likelihood = st.binom(n, thetas).pmf(h)
post = prior * likelihood
post /= (post.sum() / len(thetas))

最后,我们将先验、后验和似然的概率密度函数可视化。

plt.figure(figsize=(12, 9))
plt.plot(thetas, prior, label='Prior', c='blue')
plt.plot(thetas, n*likelihood, label='Likelihood', c='green')
plt.plot(thetas, post, label='Posterior', c='red')
plt.xlim([0, 1])
plt.xlabel(r'$\theta$', fontsize=14)
plt.ylabel('PDF', fontsize=16)
plt.legend();
b67679306ddaeb8dba477166d9e655d6.png

结论

蒙特卡罗积分是求解积分的一种数值方法。它的工作原理是在随机点对函数求值,求和所述值,然后计算它们的平均值。

原文链接:https://towardsdatascience.com/monte-carlo-integration-db86b8d7beb3

下载1:OpenCV-Contrib扩展模块中文版教程

在「小白学视觉」公众号后台回复:扩展模块中文教程即可下载全网第一份OpenCV扩展模块教程中文版,涵盖扩展模块安装、SFM算法、立体视觉、目标跟踪、生物视觉、超分辨率处理等二十多章内容。

下载2:Python视觉实战项目52讲

在「小白学视觉」公众号后台回复:Python视觉实战项目即可下载包括图像分割、口罩检测、车道线检测、车辆计数、添加眼线、车牌识别、字符识别、情绪检测、文本内容提取、面部识别等31个视觉实战项目,助力快速学校计算机视觉。

下载3:OpenCV实战项目20讲

在「小白学视觉」公众号后台回复:OpenCV实战项目20讲即可下载含有20个基于OpenCV实现20个实战项目,实现OpenCV学习进阶。

交流群

欢迎加入公众号读者群一起和同行交流,目前有SLAM、三维视觉、传感器、自动驾驶、计算摄影、检测、分割、识别、医学影像、GAN、算法竞赛等微信群(以后会逐渐细分),请扫描下面微信号加群,备注:”昵称+学校/公司+研究方向“,例如:”张三 + 上海交大 + 视觉SLAM“。请按照格式备注,否则不予通过。添加成功后会根据研究方向邀请进入相关微信群。请勿在群内发送广告,否则会请出群,谢谢理解~

38d85a1feeb80564262c4374b9f5dc70.png

c39911548519f57c90ea52120e337ced.png

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值