蒙特卡羅計算積分

時間 2020-10-21

標籤 python dom 機器學習函數 post 學習測試 spa .net 3d 欄目 Python 简体版

原文原文鏈接

做者|Cory Maklin
編譯|VK
來源|Towards Datas Sciencepython

一般狀況下，咱們不能解析地求解積分，必須藉助其餘方法，其中就包括蒙特卡羅積分。你可能還記得，函數的積分能夠解釋爲函數曲線下的面積。dom

蒙特卡羅積分的工做原理是在a和b之間的不一樣隨機點計算一個函數，將矩形的面積相加，取和的平均值。隨着點數的增長，所得結果接近於積分的實際解。機器學習

蒙特卡羅積分用代數表示：函數

與其餘數值方法相比，蒙特卡羅積分特別適合於計算奇數形狀的面積。post

在上一節中，咱們看到如何使用蒙特卡羅積分來肯定後驗機率，當咱們知道先驗和似然，但缺乏規範化常數。學習

貝葉斯統計

後驗機率是指貝葉斯公式中的一個特定項。測試

貝葉斯定理能夠用來計算一我的在某一特定疾病的篩查測試中呈陽性的人實際上患有該病的機率。spa

若是咱們已經知道P（A），P（B）和P（B | A），但想知道P（A | B），咱們就用這個公式。例如，假設咱們正在檢測一種感染1%人口的罕見疾病。醫學專業人員已經開發出一種高度敏感和特異的檢測方法，但還不夠完善。.net

99%的病人檢測呈陽性
99%的健康患者檢測爲陰性

貝葉斯定理告訴咱們：3d

假設咱們有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();