搜尋
首頁後端開發Python教學python程式設計透過蒙特卡羅法計算定積分詳解

這篇文章主要介紹了python程式設計透過蒙特卡羅法計算定積分詳解,具有一定借鑒價值,需要的朋友可以參考下。

想當初,考研的時候要是知道有這麼好東西,計算定積分。 。 。開玩笑,那時候計算定積分根本沒有這麼簡單的。但這確實給我開啟了一個思路,用程式語言解決更多更複雜的數學問題。下面進入正題。

如上圖所示,計算區間[a b]上f(x)的積分即求曲線與X軸圍成紅色區域的面積。以下使用蒙特卡羅法計算區間[2 3]上的定積分:∫(x2 4*x*sin(x))dx


##

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

def f(x):
  return x**2 + 4*x*np.sin(x) 
def intf(x): 
  return x**3/3.0+4.0*np.sin(x) - 4.0*x*np.cos(x)
a = 2;  
b = 3; 
# use N draws 
N= 10000
X = np.random.uniform(low=a, high=b, size=N) # N values uniformly drawn from a to b 
Y =f(X)  # CALCULATE THE f(x) 
# 蒙特卡洛法计算定积分:面积=宽度*平均高度
Imc= (b-a) * np.sum(Y)/ N;
exactval=intf(b)-intf(a)
print "Monte Carlo estimation=",Imc, "Exact number=", intf(b)-intf(a)
# --How does the accuracy depends on the number of points(samples)? Lets try the same 1-D integral 
# The Monte Carlo methods yield approximate answers whose accuracy depends on the number of draws.
Imc=np.zeros(1000)
Na = np.linspace(0,1000,1000)
exactval= intf(b)-intf(a)
for N in np.arange(0,1000):
  X = np.random.uniform(low=a, high=b, size=N) # N values uniformly drawn from a to b 
  Y =f(X)  # CALCULATE THE f(x) 
  Imc[N]= (b-a) * np.sum(Y)/ N;   
plt.plot(Na[10:],np.sqrt((Imc[10:]-exactval)**2), alpha=0.7)
plt.plot(Na[10:], 1/np.sqrt(Na[10:]), 'r')
plt.xlabel("N")
plt.ylabel("sqrt((Imc-ExactValue)$^2$)")
plt.show()


>>>

Monte Carlo estimation= 11.8181144118 Exact number= 11.8113589251

從上圖可以看出,隨著取樣點數的增加,計算誤差逐漸減少。想要提高模擬結果的精確度有兩個途徑:其一是增加試驗次數N;其二是降低方差σ2. 增加試驗次數勢必使解題所用計算機的總時間增加,要想以此來達到提高精度之目的顯然是不合適的。以下來介紹重要抽樣法來減少方差,提高積分計算的精確度。

重要性抽樣法的特點在於,它不是從給定的過程的機率分佈抽樣,而是從修改的機率分佈抽樣,使對模擬結果有重要作用的事件更多出現,從而提高抽樣效率,減少花在對模擬結果無關的事件上的計算時間。例如在區間[a b]上求g(x)的積分,若採用均勻抽樣,在函數值g(x)比較小的區間內產生的抽樣點跟函數值較大處區間內產生的抽樣點的數目接近,顯然抽樣效率不高,可以將抽樣機率密度函數改為f(x),使f(x)與g(x)的形狀相近,就可以保證對積分計算貢獻較大的抽樣值出現的機會大於貢獻小的抽樣值,即可以將積分運算改寫為:

#x是依照機率密度f(x)抽樣所得的隨機變量,顯然在區間[a b ]內應該有:

因此,可容易將積分值I看成是隨機變數Y = g(x)/f(x)的期望,式子​​中xi是服從機率密度f(x)的取樣點

下面的例子採用常態分佈函數f(x)來近似g(x)=sin(x )*x,並依據常態分佈選取取樣值計算區間[0 pi]上的積分個∫g(x)dx


# -*- coding: utf-8 -*-
# Example: Calculate ∫sin(x)xdx

# The function has a shape that is similar to Gaussian and therefore
# we choose here a Gaussian as importance sampling distribution.
from scipy import stats
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
mu = 2;
sig =.7;
f = lambda x: np.sin(x)*x
infun = lambda x: np.sin(x)-x*np.cos(x)
p = lambda x: (1/np.sqrt(2*np.pi*sig**2))*np.exp(-(x-mu)**2/(2.0*sig**2))
normfun = lambda x: norm.cdf(x-mu, scale=sig)

plt.figure(figsize=(18,8)) # set the figure size
# range of integration
xmax =np.pi 
xmin =0
# Number of draws 
N =1000
# Just want to plot the function
x=np.linspace(xmin, xmax, 1000)
plt.subplot(1,2,1)
plt.plot(x, f(x), 'b', label=u'Original $x\sin(x)$')
plt.plot(x, p(x), 'r', label=u'Importance Sampling Function: Normal')
plt.xlabel('x')
plt.legend()
# =============================================
# EXACT SOLUTION 
# =============================================
Iexact = infun(xmax)-infun(xmin)
print Iexact
# ============================================
# VANILLA MONTE CARLO 
# ============================================
Ivmc = np.zeros(1000)
for k in np.arange(0,1000):
  x = np.random.uniform(low=xmin, high=xmax, size=N)
  Ivmc[k] = (xmax-xmin)*np.mean(f(x))
# ============================================
# IMPORTANCE SAMPLING 
# ============================================
# CHOOSE Gaussian so it similar to the original functions

# Importance sampling: choose the random points so that
# more points are chosen around the peak, less where the integrand is small.
Iis = np.zeros(1000)
for k in np.arange(0,1000):
  # DRAW FROM THE GAUSSIAN: xis~N(mu,sig^2)
  xis = mu + sig*np.random.randn(N,1);
  xis = xis[ (xis<xmax) & (xis>xmin)] ;
  # normalization for gaussian from 0..pi
  normal = normfun(np.pi)-normfun(0)   # 注意:概率密度函数在采样区间[0 pi]上的积分需要等于1
  Iis[k] =np.mean(f(xis)/p(xis))*normal  # 因此,此处需要乘一个系数即p(x)在[0 pi]上的积分
plt.subplot(1,2,2)
plt.hist(Iis,30, histtype=&#39;step&#39;, label=u&#39;Importance Sampling&#39;);
plt.hist(Ivmc, 30, color=&#39;r&#39;,histtype=&#39;step&#39;, label=u&#39;Vanilla MC&#39;);
plt.vlines(np.pi, 0, 100, color=&#39;g&#39;, linestyle=&#39;dashed&#39;)
plt.legend()
plt.show()


#從圖中可以看出曲線sin(x)*x的形狀和常態分佈曲線的形狀相近,因此在曲線峰值處的取樣點數目會比曲線上位置低的地方要多。精確計算的結果為pi,從上面的右圖中可以看出:兩種方法均計算定積分1000次,靠近精確值pi=3.1415處的結果最多,離精確值越遠數目越少,顯然這符合常規。但是採用傳統方法(紅色直方圖)計算出的積分值方的差明顯比採用重要抽樣法(藍色直方圖)更大。因此,採用重要抽樣法計算可以降低方差,提高精確度。另外要注意的是:關於函數f(x)的選擇會對計算結果的精確度產生影響,當我們選擇的函數f(x)與g(x)相差較大時,計算結果的變異數也會加大。

相關推薦:


Python程式設計中NotImplementedError的使用方法_python


以上是python程式設計透過蒙特卡羅法計算定積分詳解的詳細內容。更多資訊請關注PHP中文網其他相關文章!

陳述
本文內容由網友自願投稿,版權歸原作者所有。本站不承擔相應的法律責任。如發現涉嫌抄襲或侵權的內容,請聯絡admin@php.cn
详细讲解Python之Seaborn(数据可视化)详细讲解Python之Seaborn(数据可视化)Apr 21, 2022 pm 06:08 PM

本篇文章给大家带来了关于Python的相关知识,其中主要介绍了关于Seaborn的相关问题,包括了数据可视化处理的散点图、折线图、条形图等等内容,下面一起来看一下,希望对大家有帮助。

详细了解Python进程池与进程锁详细了解Python进程池与进程锁May 10, 2022 pm 06:11 PM

本篇文章给大家带来了关于Python的相关知识,其中主要介绍了关于进程池与进程锁的相关问题,包括进程池的创建模块,进程池函数等等内容,下面一起来看一下,希望对大家有帮助。

Python自动化实践之筛选简历Python自动化实践之筛选简历Jun 07, 2022 pm 06:59 PM

本篇文章给大家带来了关于Python的相关知识,其中主要介绍了关于简历筛选的相关问题,包括了定义 ReadDoc 类用以读取 word 文件以及定义 search_word 函数用以筛选的相关内容,下面一起来看一下,希望对大家有帮助。

归纳总结Python标准库归纳总结Python标准库May 03, 2022 am 09:00 AM

本篇文章给大家带来了关于Python的相关知识,其中主要介绍了关于标准库总结的相关问题,下面一起来看一下,希望对大家有帮助。

分享10款高效的VSCode插件,总有一款能够惊艳到你!!分享10款高效的VSCode插件,总有一款能够惊艳到你!!Mar 09, 2021 am 10:15 AM

VS Code的确是一款非常热门、有强大用户基础的一款开发工具。本文给大家介绍一下10款高效、好用的插件,能够让原本单薄的VS Code如虎添翼,开发效率顿时提升到一个新的阶段。

python中文是什么意思python中文是什么意思Jun 24, 2019 pm 02:22 PM

pythn的中文意思是巨蟒、蟒蛇。1989年圣诞节期间,Guido van Rossum在家闲的没事干,为了跟朋友庆祝圣诞节,决定发明一种全新的脚本语言。他很喜欢一个肥皂剧叫Monty Python,所以便把这门语言叫做python。

Python数据类型详解之字符串、数字Python数据类型详解之字符串、数字Apr 27, 2022 pm 07:27 PM

本篇文章给大家带来了关于Python的相关知识,其中主要介绍了关于数据类型之字符串、数字的相关问题,下面一起来看一下,希望对大家有帮助。

详细介绍python的numpy模块详细介绍python的numpy模块May 19, 2022 am 11:43 AM

本篇文章给大家带来了关于Python的相关知识,其中主要介绍了关于numpy模块的相关问题,Numpy是Numerical Python extensions的缩写,字面意思是Python数值计算扩展,下面一起来看一下,希望对大家有帮助。

See all articles

熱AI工具

Undresser.AI Undress

Undresser.AI Undress

人工智慧驅動的應用程序,用於創建逼真的裸體照片

AI Clothes Remover

AI Clothes Remover

用於從照片中去除衣服的線上人工智慧工具。

Undress AI Tool

Undress AI Tool

免費脫衣圖片

Clothoff.io

Clothoff.io

AI脫衣器

AI Hentai Generator

AI Hentai Generator

免費產生 AI 無盡。

熱門文章

R.E.P.O.能量晶體解釋及其做什麼(黃色晶體)
2 週前By尊渡假赌尊渡假赌尊渡假赌
倉庫:如何復興隊友
4 週前By尊渡假赌尊渡假赌尊渡假赌
Hello Kitty Island冒險:如何獲得巨型種子
4 週前By尊渡假赌尊渡假赌尊渡假赌

熱工具

VSCode Windows 64位元 下載

VSCode Windows 64位元 下載

微軟推出的免費、功能強大的一款IDE編輯器

SublimeText3 Linux新版

SublimeText3 Linux新版

SublimeText3 Linux最新版

記事本++7.3.1

記事本++7.3.1

好用且免費的程式碼編輯器

EditPlus 中文破解版

EditPlus 中文破解版

體積小,語法高亮,不支援程式碼提示功能

禪工作室 13.0.1

禪工作室 13.0.1

強大的PHP整合開發環境