這篇文章主要介紹了python編程通過(guò)蒙特卡洛法計(jì)算定積分詳解，具有一定借鑒價(jià)值，需要的朋友可以參考下。

想當(dāng)初，考研的時(shí)候要是知道有這么個(gè)好東西，計(jì)算定積分。。。開玩笑，那時(shí)候計(jì)算定積分根本沒(méi)有這么簡(jiǎn)單的。但這確實(shí)給我打開了一種思路，用編程語(yǔ)言去解決更多更復(fù)雜的數(shù)學(xué)問(wèn)題。下面進(jìn)入正題。

如上圖所示，計(jì)算區(qū)間[a b]上f(x)的積分即求曲線與X軸圍成紅色區(qū)域的面積。下面使用蒙特卡洛法計(jì)算區(qū)間[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)

# 蒙特卡洛法計(jì)算定積分：面積=寬度*平均高度

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

從上圖可以看出，隨著采樣點(diǎn)數(shù)的增加，計(jì)算誤差逐漸減小。想要提高模擬結(jié)果的精確度有兩個(gè)途徑：其一是增加試驗(yàn)次數(shù)N；其二是降低方差σ2. 增加試驗(yàn)次數(shù)勢(shì)必使解題所用計(jì)算機(jī)的總時(shí)間增加，要想以此來(lái)達(dá)到提高精度之目的顯然是不合適的。下面來(lái)介紹重要抽樣法來(lái)減小方差，提高積分計(jì)算的精度。

重要性抽樣法的特點(diǎn)在于，它不是從給定的過(guò)程的概率分布抽樣，而是從修改的概率分布抽樣，使對(duì)模擬結(jié)果有重要作用的事件更多出現(xiàn)，從而提高抽樣效率，減少花費(fèi)在對(duì)模擬結(jié)果無(wú)關(guān)緊要的事件上的計(jì)算時(shí)間。比如在區(qū)間[a b]上求g(x)的積分，若采用均勻抽樣，在函數(shù)值g(x)比較小的區(qū)間內(nèi)產(chǎn)生的抽樣點(diǎn)跟函數(shù)值較大處區(qū)間內(nèi)產(chǎn)生的抽樣點(diǎn)的數(shù)目接近，顯然抽樣效率不高，可以將抽樣概率密度函數(shù)改為f(x)，使f(x)與g(x)的形狀相近，就可以保證對(duì)積分計(jì)算貢獻(xiàn)較大的抽樣值出現(xiàn)的機(jī)會(huì)大于貢獻(xiàn)小的抽樣值，即可以將積分運(yùn)算改寫為：

x是按照概率密度f(wàn)(x)抽樣獲得的隨機(jī)變量，顯然在區(qū)間[a b]內(nèi)應(yīng)該有：

因此，可容易將積分值I看成是隨機(jī)變量 Y = g(x)/f(x)的期望，式子中xi是服從概率密度f(wàn)(x)的采樣點(diǎn)

下面的例子采用一個(gè)正態(tài)分布函數(shù)f(x)來(lái)近似g(x)=sin(x)*x，并依據(jù)正態(tài)分布選取采樣值計(jì)算區(qū)間[0 pi]上的積分個(gè)∫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[ (xisxmin)] ;

# normalization for gaussian from 0..pi

normal = normfun(np.pi)-normfun(0) # 注意:概率密度函數(shù)在采樣區(qū)間[0 pi]上的積分需要等于1

Iis[k] =np.mean(f(xis)/p(xis))*normal # 因此,此處需要乘一個(gè)系數(shù)即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()

從圖中可以看出曲線sin(x)*x的形狀和正態(tài)分布曲線的形狀相近，因此在曲線峰值處的采樣點(diǎn)數(shù)目會(huì)比曲線上位置低的地方要多。精確計(jì)算的結(jié)果為pi，從上面的右圖中可以看出：兩種方法均計(jì)算定積分1000次，靠近精確值pi=3.1415處的結(jié)果最多，離精確值越遠(yuǎn)數(shù)目越少，顯然這符合常規(guī)。但是采用傳統(tǒng)方法(紅色直方圖)計(jì)算出的積分值方的差明顯比采用重要抽樣法(藍(lán)色直方圖)要大。因此，采用重要抽樣法計(jì)算可以降低方差，提高精度。另外需要注意的是：關(guān)于函數(shù)f(x)的選擇會(huì)對(duì)計(jì)算結(jié)果的精度產(chǎn)生影響，當(dāng)我們選擇的函數(shù)f(x)與g(x)相差較大時(shí)，計(jì)算結(jié)果的方差也會(huì)加大。

相關(guān)推薦：

以上就是python編程通過(guò)蒙特卡洛法計(jì)算定積分詳解的詳細(xì)內(nèi)容，更多請(qǐng)關(guān)注php中文網(wǎng)其它相關(guān)文章！

本文原創(chuàng)發(fā)布php中文網(wǎng)，轉(zhuǎn)載請(qǐng)注明出處，感謝您的尊重！

總結(jié)

以上是生活随笔為你收集整理的python计算定积分_python编程通过蒙特卡洛法计算定积分详解的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。