蒙特卡洛模擬方法

       花費兩天的時間去幫老師研究課題，完成那一刻瞬間感覺自己學到的數學知識很有用。於是想著以一種方式去記錄這次算法。思來想去寫了這個公眾號，作為第一篇自我感覺有價值的文章，我會一點一滴的去寫。但是由於老師的課題需要發表論文，所以本文不會給出具體關於課題的任何代碼。但是會用python去構思基本的東西出來。我希望此文可以記錄思想，整理筆記、知識，並將其中承載的價值傳播給他人。 
在本文中您可以學習到以下幾點知識

蒙特卡洛定積分算法原理講解（公式篇）

python實現算法（代碼篇）

蒙特卡洛算法起源

       蒙特卡羅方法於20世紀40年代美國在第二次世界大戰中研製原子彈的「曼哈頓計劃」計劃的成員S.M.烏拉姆和J.馮·諾伊曼首先提出。數學家馮·諾伊曼用馳名世界的賭城—摩納哥的Monte Carlo—來命名這種方法，為它蒙上了一層神秘色彩。在這之前，蒙特卡羅方法就已經存在。1777年，法國Buffon提出用投針實驗的方法求圓周率∏。這被認為是蒙特卡羅方法的起源。 
很簡單的一個小例子，在長寬都為1cm的正方形中有一個最大圓，如何用數理統計的思想去求解這個面積呢？實際上我們可以在正方形中扔豆子，扔很多很多次（n），然後計算圓內豆子和總豆子的比例，再乘以正方形面積，當n趨近於無窮大時，這時候正方形面積乘以豆子比例在理想情況下就是圓的面積了。基於這個思想，我們可以思考一下，如何去計算一個一維積分呢？

蒙特卡洛求解定積分原理（公式篇）

一維積分: 

     在x的定義域[0,1]中均勻隨機取點，該均勻分布的隨機變量記為η1,定義一個隨機變量η為：

則顯然有： 

     η1的期望值等於積分I。只要抽取足夠多的隨機點，即隨機點的點數足夠多大時，n的平均值f(ξ)：

就是積分Ι的無偏估計值。 
η1的方差為： 

      顯然V{η1 }依賴被積函數f(x)在積分域上的方差。當f(x)在x的定義域內變化曲率較小（平坦），就是和I的差值較小時，方差也小；反之，方差越大。 
如果想減小積分估計值得差值，就需要減小被積函數在積分區域上的方差。一個很重要的思想——重要抽樣法： 
將f(x)的方差吸收到一個特定函數g(x)中。假設f*(x)=f(x)/(g(x)在定義域內極其平緩，那麼積分公式變為： 

若選取η'為服從概率密度函數g(x)的函數f*(x)的抽樣值。則可以得到：

因此它的平均值為： 

給出一個I的無偏估計值。這時的方差為：

在實際計算中，方差可以通過以下方程得到： 

      方程中角型括號表示對括號內所有可能的[0,1]區間，按g(x)分布的隨機坐標數序列{x_i }對應的數值求平均值。方程右邊第一項對{f*²(x_i)}求平均記為a，第二項表示求{f*(x_i)}的平均值的平方,記為b。上面方程可以推導得到： 

    由此我們看出其誤差平方與f*在[0,1]區間的方差成正比，並且σ∝（1/√n）。這與中心極限定理所得到的的結果一致。

利用python計算定積分實例（代碼篇）

如下圖所示： 

利用蒙特卡洛計算

python代碼如下所示：

import numpy as npimport 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; N= 10000X = np.random.uniform(low=a, high=b, size=N) Y =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))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)   Y =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()
實驗成圖如下：
      從上圖可以看出隨著實驗次數的不斷增加，計算的誤差也越來越小。由此想要提高精確度有兩個方法：
增加實驗次數: 
隨著實驗次數的不斷增加，計算機計算的時間也越來越長，不是很符合實際情況；
重要採樣方法 
重要採樣方法的特點在於，它不是從給定的過程的概率抽樣，而是修改的概率分布抽樣，使對模擬結果有重要作用的事件更多多次出現，從而提高抽樣概率，減少花費在對模擬結果無關緊要的事件上的計算時間。比如在區間[a,b]上求g(x)的積分，若採用均勻抽樣，在函數值g(x)比較小的區間內產生的抽樣點跟函數值較大處區間內產生的抽樣點的數目接近，顯然抽樣效率不高，可以將抽樣概率密度函數改為f(x)，使f(x)與g(x)的形狀相近，就可以保證對積分計算貢獻較大的抽樣值出現的機會大於貢獻小的抽樣值，即可以將積分改為：
x是按照抽樣概率密度f(x)抽樣獲得的隨機變量，顯然在區間[a,b]上應該有: 
       因此，可以容易將積分值I看成是隨機變量Y=g(x)/f(x)的期望值，式子中服從概率密度為f(x)的採樣點： 
      下面的例子採用正太分布f(x)來近似g(x)=xsinx，並根據正太分布來選取採樣值區間為[0，π]上的積分: 
python代碼如下所示：
 from scipy import statsfrom scipy.stats import normimport numpy as npimport matplotlib.pyplot as pltmu = 2;sig =.7;f = lambda x: np.sin(x)*xinfun = 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)) xmax =np.pi xmin =0N =1000x=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()Iexact = infun(xmax)-infun(xmin)print(Iexact)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)) Iis = np.zeros(1000)for k in np.arange(0,1000):    xis = mu + sig*np.random.randn(N,1);  xis = xis[ (xis<xmax) & (xis>xmin)] ;    normal = normfun(np.pi)-normfun(0)     Iis[k] =np.mean(f(xis)/p(xis))*normal  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()
實驗效果圖如下：
      從圖中可以看出xsin⁡x的圖像與正太分布圖像很接近，因此在曲線峰值採樣點的數量要比曲線上其他低的位置多。精確計算的結果為π，從上面的右圖可以看出：兩種計算方法計算定積分1000次，精確靠近π=3.1415處的結果最多，離精確值越遠數目越少，顯然符合常規。但傳統的方法（紅色直方圖）計算出的積分值的差明顯比採用重要採樣法（藍色直方圖）要大。因此採用重要採樣法計算可以降低方差，提高精度。另外需要主要的是：關於f(x)的選擇會對計算結果的精度產生影響，當選擇的函數f(x)與g(x)相差較大時，計算結果的方差也會增大。
到此為止，蒙特卡洛算法從原理到實驗講述完畢，不知道我講明白沒有呢？
我們一直都在努力前行，將科學與現實聯繫在一起，將藝術與生活緊密結合