python程式設計透過蒙特卡羅法計算定積分詳解

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

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

如上圖所示，計算區間[a b]上f(x)的積分即求曲線與X軸圍成紅色區域的面積。以下使用蒙特卡羅法計算區間[2 3]上的定積分：∫(x² 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中文網其他相關文章！