這篇文章主要介紹了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()
# -*- 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='step', label=u'Importance Sampling'); plt.hist(Ivmc, 30, color='r',histtype='step', label=u'Vanilla MC'); plt.vlines(np.pi, 0, 100, color='g', linestyle='dashed') plt.legend() plt.show()
Python程式設計中NotImplementedError的使用方法_python
以上是python程式設計透過蒙特卡羅法計算定積分詳解的詳細內容。更多資訊請關注PHP中文網其他相關文章!

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

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

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

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

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

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

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


熱AI工具

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

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

Undress AI Tool
免費脫衣圖片

Clothoff.io
AI脫衣器

AI Hentai Generator
免費產生 AI 無盡。

熱門文章

熱工具

VSCode Windows 64位元 下載
微軟推出的免費、功能強大的一款IDE編輯器

SublimeText3 Linux新版
SublimeText3 Linux最新版

記事本++7.3.1
好用且免費的程式碼編輯器

EditPlus 中文破解版
體積小,語法高亮,不支援程式碼提示功能

禪工作室 13.0.1
強大的PHP整合開發環境