蒙特卡罗计算积分

做者|Cory Maklin
编译|VK
来源|Towards Datas Science

一般状况下，咱们不能解析地求解积分，必须借助其余方法，其中就包括蒙特卡罗积分。你可能还记得，函数的积分能够解释为函数曲线下的面积。

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

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

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

在上一节中，咱们看到如何使用蒙特卡罗积分来肯定后验几率，当咱们知道先验和似然，但缺乏规范化常数。

贝叶斯统计

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

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

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

• 99%的病人检测呈阳性
• 99%的健康患者检测为阴性

贝叶斯定理告诉咱们：

假设咱们有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

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

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

然而，在现实世界中，认为0.001的几率事实上如此精确是不合理的。一个给定的人患病的几率会因其年龄、性别、体重等而有很大差别。通常来讲，咱们对给定先验几率的认识还远远不够完善，由于它是从之前的样本中得出的（这意味着不一样的人群可能会对先验几率给出不一样的估计）。

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

归一化常数

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

让咱们来探究一下为何咱们首先须要一个标准化常数。在几率论中，规范化常数是一个函数必须乘以的常数，所以它的图下面积为1。仍是不清楚？让咱们看一个例子。

回想一下正态分布的函数能够写成：

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

让咱们来看看咱们是如何肯定它的。咱们从如下函数开始（假设均值为0，方差为1）：

若是咱们能画出一个曲线的话。

问题在于，若是咱们取曲线下的面积，它不等于1，这要求它是一个几率密度函数。所以，咱们将函数除以积分的结果（归一化常数）。

回到手头的问题，即如何在没有归一化常数的状况下计算后验几率……事实证实，对于连续样本空间，规范化常数能够重写为：

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

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。所以，咱们能够简化方程。

在代码中，前面的等式写以下：

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

结论

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

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

欢迎关注磐创AI博客站：
http://panchuang.net/

sklearn机器学习中文官方文档：
http://sklearn123.com/

欢迎关注磐创博客资源汇总站：
http://docs.panchuang.net/