python求定積分的函式 蒙特卡羅計算積分

2021-10-13 20:43:49 字數 3463 閱讀 1578

作者|cory maklin 編譯|vk **|towards datas science

通常情況下,我們不能解析地求解積分,必須借助其他方法,其中就包括蒙特卡羅積分。你可能還記得,函式的積分可以解釋為函式曲線下的面積。

蒙特卡羅積分的工作原理是在a和b之間的不同隨機點計算乙個函式,將矩形的面積相加,取和的平均值。隨著點數的增加,所得結果接近於積分的實際解。

蒙特卡羅積分用代數表示:

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

在上一節中,我們看到如何使用蒙特卡羅積分來確定後驗概率,當我們知道先驗和似然,但缺少規範化常數。

貝葉斯統計

後驗概率是指貝葉斯公式中的乙個特定項。

貝葉斯定理可以用來計算乙個人在某一特定疾病的篩查測試中呈陽性的人實際上患有該病的概率。

如果我們已經知道p(a),p(b)和p(b | a),但想知道p(a | b),我們就用這個公式。例如,假設我們正在檢測一種感染1%人口的罕見疾病。醫學專業人員已經開發出一種高度敏感和特異的檢測方法,但還不夠完善。

貝葉斯定理告訴我們:

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

結論蒙特卡羅積分是求解積分的一種數值方法。它的工作原理是在隨機點對函式求值,求和所述值,然後計算它們的平均值。

python蒙特卡洛求定積分 蒙特卡洛定積分(一)

一 蒙特卡洛模擬法分類 蒙特卡洛法模擬法從其應用方面來劃分,可以分成以下三種形式 1 直接蒙特卡洛模擬。採用隨機數學咧來模擬複雜隨機過程的效應。2 蒙特卡洛定積分 間接蒙特卡洛模擬 利用隨機數序列計算積分的方法。積分維數越高,該方法的積分效率就越高。3 metropolis蒙特卡洛模擬。以 馬爾可夫...

python Matlab求定積分

計算 sympy庫中integrate函式 integrate f,x,lower bound,upper bound f 函式,x 變數,lower bound 下限,upper bound 上限 但是發現求不出來,如果是sin 2 x 就可以,為什麼?syms x f sin 2 x 1 x 2...

矩形法求定積分通用函式

題目 實現求sin,cos,exp的通用函式 思路 其實就是練習指向函式的指標 1 include2 include3 intmain 28 29 30double fsin float x1 33 double fcos float x1 36 double fexp float x1 39 vo...