Python & Matplotlib: Monte Carlos Method

　　Hey! 這裏是Lindy:) Hope you guys are doing well!  

　　今天想記錄的概念叫作 蒙特·卡羅 方法，是今年在cs課上老師作的擴展延伸。其實我在初次接觸這個概念時以爲很新奇，由於做爲一個編程菜鳥，在python裏試圖計算時（這裏指數學運算，也就是說output是以float,或integer的形式來表示），通常依賴於python的math module來作出肯定的計算。可是蒙特卡羅方法（如下簡稱爲MC method）卻帶給我了徹底不一樣的思路。

　　話很少說，先上從萬能的Wikipedia複製來的概念：

Monte Carlo methods, or Monte Carlo experiments, are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be deterministic in principle.

　　簡單的來講就是使用純粹機率統計的方法，在有必定數量的隨機數的統計基礎上計算出有高精確度的數據結果。

　　emmm...懂了麼？沒懂。

　　不要緊，咱們來看一個例子。

---

 

　　若是咱們想要得到 π 的精確數值（I know engineers use 3 as approximation for pi :) we are not doing that here)，最直接的方法是引入 python內置的math module 裏儲存的π，以下：　

1 import math
2 print(math.pi)

 　　這個方法的output是 3.141592653589793。可是須要說起的是這個數字是提早被儲存爲一個常數，因此這個方法嚴格的來講不算是計算，而是調用。

　　可是若是咱們不使用math module呢？你可能會說咱們也能夠用 Taylor Polynomial 來近似π啊。泰勒公式使用以下

$\mathrm{\pi}=4\cdot \mathrm{arc}\tan \left( 1 \right) =4\cdot \sum_{\mathrm{i}=0}^{\mathrm{N}}{\frac{\left( -1 \right) ^{\mathrm{i}-1}}{\left( 2\mathrm{i}-1 \right) ^{-1}}}$

　　以上公式能夠用python表達以下：

1 sum = 0
2 for i in range(1,100000):
3     sum += ((-1)**(i-1))*((2*i-1)**-1)
4 sum *= 4
5 print(sum)

　　這個方法的output是 3.1416026536897204，能夠看到，在循環了 99999次後這個結果開始接近了math module 儲存的π，但在小數點四位後便開始偏離了π的真實數值。

　　那還有什麼方法呢？好了我要點題了，咱們要用 Monte Carlos Method 了。

---

 

　　根據定義，MC method是根據在機率統計學上，因此須要必定的隨機數字做爲基礎量來計算機率密度。因此第一步咱們要製做不少隨機數，這一步須要用到python的random module：

1 import random
2 x , y = random.random() , random.random()

　　*random.random() 會隨機產生一個在 0-1之間的浮點數。

　　在進行下一步以前，咱們須要簡要複習一下π的定義。提起π，咱們都會很天然的把這個數字跟圓的周長和麪積進行聯繫。在這裏，爲了把機率密度視覺化，咱們要使用面積的概念：A = πR²。下圖是一個邊長爲1的正方形，內切1/4的半徑爲1的圓。

　　那若是咱們在以上圖片平面內隨機選一點的話，這一點落入圓的面積範圍的機率是多少呢？答案是$\mathrm{P}\left( \mathrm{x} \right) =\frac{\mathrm{A}_{\mathrm{circle}}}{\mathrm{A}_{\mathrm{square}}}=\frac{\frac{\mathrm{\pi}}{4}}{1}=\frac{\mathrm{\pi}}{4}$。利用這個概念，若是咱們可以計算點落入圓區域的機率(並乘以4)，咱們就能估算 π 的值。爲了實現以上概念，咱們要使用相似於以上 Talyor Polynomial 的循環寫法（這裏爲了模仿這個循環寫法做精度對比，我用了for loop，可是這裏while loop也能夠實現一樣的效果):

 1 import random
 2 
 3 total, total_random = 0, 99999 # 初始化
 4 
 5 for i in range(total_random):
 6     # 生成隨機數
 7     x , y = random.random() , random.random()
 8     # 計算點是否落入圓的範圍
 9     if (x**2 + y**2)**(1/2)<1: 
10         total += 1
11 # 計算pi的值
12 pi = 4*(total/total_random)
13 print(pi)

　　由於機率統計計算的不穩定性，我重複跑了5遍一樣的代碼，但生出的結果卻大不同。在平均了100次代碼產出結果後，輸出結果爲 3.141013010130101。可見，這個計算方法在小數點後3位後便偏離了 π的精確數值。在把隨機數量提高到100000後，精確數字增長到了小數點後4位。可是比起泰勒公式的精確計算，由於統計學的不肯定性，MC method的精確度確有欠缺。MC Method 的顯著優勢是實現度高，沒有針對數學理解的要求，並且可視度高，理解上更容易一些。

---

　　用一樣的方法，咱們能夠用相似的機率面積估算來計算函數的積分。這是由於 Integral = Area under the Curve, 因此咱們能夠用相似的面積計算方式來解決這個問題。由於面積是一個有限的量，因此在這裏咱們計算的是definite integral。

模仿以上π的MC算法，咱們能夠試着來用python實現$\int_0^1{\mathrm{x}^2}$ 注意這裏咱們能夠繼續使用random.random()，由於爲了代碼的簡潔性，我假設了這個積分的函數的local max在1之內

 1 total, total_random = 0, 999999 #初始化
 2 
 3 def f(x): #定義須要計算的積分函數，在這個例子裏是x^2
 4     return x**2
 5 
 6 for i in range(total_random): #產生大量隨機數
 7     x , y = random.random(), random.random()
 8     if y < f(x): #若是隨機點在函數下方，則計入總數
 9         total += 1
10 
11 integral = (total/total_random) #計算機率密度
12 print(integral)

　　代碼運行結果爲0.3331783331783332，和理論結果（1/3）仍是很類似的。這也證實了MC Method在積分計算上的可行性。以上代碼模仿了初始的π的MC算法，咱們能夠看到，這個例子的變形主要是在函數的定義上；咱們從新寫了一個函數，定義了須要返回的函數的具體值。如此，這個函數才能在以後的隨機數計算是被用於歸類（這裏的歸類是指看點是否落在函數下方）。

　　可是以上代碼只適用於$\int_0^1{\mathrm{x}^2}$這個特定的例子，若是咱們要計算更復雜的函數，或者不一樣界限的函數，這些代碼將再也不適用（應該說怎麼會可以適用呢：）連函數可能都不同了。）另外代碼的計算過程並無被實際的表達出來，也就是說即便經過代碼咱們能夠成功的得到積分的最終計算值，可僅僅經過這個值咱們也沒有辦法得到關於這個函數的其餘信息，例如形狀，最大，最小值等等。

　　因此爲了提高以上代碼的實用性，咱們能夠利用python經常使用的畫圖模塊，matplotlib,以及數據處理模塊，numpy進行如下處理 (具體代碼使用請看註釋，文中將再也不對代碼做用進行解釋)：

 1 import random
 2 import math
 3 import numpy as np
 4 import matplotlib.pyplot as plt
 5 
 6 def f(x): 
 7     return np.sin(x)
 8 
 9 def definite_integral_calculation(f,a,b,N):
10     '''
11     Approximate the definite integral value of the function, f, with the left bound a
12     and right bound b, using a total number of N points
13 
14     Parameters:
15         f: a defined function that is being approximated
16         a: left bound of the function
17         b: right bound of the function
18         N: the total number of random points
19     '''
20     # 首先，咱們爲了確保MC Method的正確性，要如今這裏用Numpy寫出準確函數的圖像表達。
21     act_x = np.arange(a, b, 0.01) #以0.01個單位爲間隔，產出array
22     act_y = f(act_x) #計算以上array相應的函數值
23 
24     #咱們注意到，在計算積分時，函數不必定老是落在1的下方。並且在函數同時有
25     #正負值時，咱們也要對此進行鍼對性的處理。因此咱們要對極值作出如下操做：
26     f_max = max(act_y)
27     f_min = min(act_y)
28 
29     if f_min <= 0:
30         negative = True
31     else:
32         negative = False
33 
34     #如今咱們要生成相應的隨機數並進行分類，看是否落進函數造成的區域
35     x_rand = a + np.random.random(N)*(b-a) #生成符合函數範圍的隨機數
36     y_rand = 0 + np.random.random(N)*f_max
37 
38     ind_in_pos = np.where(y_rand < f(x_rand)) #進行條件判斷，看點是否落入函數範圍內
39     ind_out_pos = np.where(y_rand >= f(x_rand))
40 
41     #咱們能夠把以上的代碼用matplotlib來具現化
42     plt.plot(act_x, act_y, color='red') #準確的函數用紅色的光滑曲線來表示
43     plt.scatter(x_rand[ind_in_pos], y_rand[ind_in_pos],color ='green') #落入範圍的隨機點用綠色點來表示
44     plt.scatter(x_rand[ind_out_pos], y_rand[ind_out_pos],color ='yellow') #未落入範圍的隨機點用藍色來表示
45 
46     #同時咱們也要考慮到同時有正負值的狀況
47     if negative:
48         y_rand_neg = 0 + np.random.random(N)*f_min #同上，生成負的隨機數
49 
50         ind_in_neg = np.where(y_rand_neg > f(x_rand)) #同上，條件判斷
51         ind_out_neg = np.where(y_rand_neg <= f(x_rand))
52 
53         plt.scatter(x_rand[ind_in_neg], y_rand_neg[ind_in_neg],color ='green') #同上，根據分類畫出隨機點
54         plt.scatter(x_rand[ind_out_neg], y_rand_neg[ind_out_neg],color ='yellow')
55 
56     #最後，咱們要計算並輸出函數相關的計算結果
57     positive_area = f_max*(b-a) #計算正的總面積
58     if negative:
59         negative_area = f_min*(b-a) #計算負的總面積
60         area = (len(ind_in_pos[0])/N)*positive_area+(len(ind_in_neg[0])/N)*negative_area #利用機率計算複合面積
61     else:
62         area = (len(ind_in_pos[0])/N)*positive_area #若是隻有正數，則計算正的面積
63         
64     #輸出結果
65     print(f'The value of definite integral calculation is {area}')
66     print(f'The approximated maximum is {f_max}')
67     print(f'The approximated minimum is {f_min}')
68 
69 
70 if __name__ == '__main__':
71     definite_integral_calculation(f,-math.pi,math.pi*0.5,99999)

　　以上代碼使用MC Method計算了$\int_{-\mathrm{\pi}}^{\frac{\mathrm{\pi}}{2}}{\sin \left( \mathrm{x} \right) \,\,\mathrm{dx}}$的具體數值，並用matplotlib具現化。首先咱們先來看看打印出來的圖像是否符合咱們預期：

　　顯然，sine 波形的形狀已經被綠色和黃色勾勒了出來，符合紅色波形的預計。而代碼的具體輸出結果

 

• The value of definite integral calculation is -0.9978134302969512
• The approximated maximum is 0.9999971463877178
• The approximated minimum is -0.9999996829318346

　　也在數學的角度上符合了咱們的預判（sine的最大值爲1，最小值爲-1，$\int_{-\mathrm{\pi}}^{\frac{\mathrm{\pi}}{2}}{\sin \left( \mathrm{x} \right) \,\,\mathrm{dx}}$的結果爲-1）。到此，在保留必定結果精度的狀況下，咱們能夠判斷在積分計算方面，Monte Carlos Method是足夠準確的。

　　最後，咱們仍是要指出最後這個方程的幾個使用時的注意事項：

1. 在計算積分時，咱們假設了這個函數是連續的，也就說若是在某一點方程會出現不平滑過渡的情況，最終呈現的函數圖像（這裏指紅色的預判線）將會再也不準確。
2. 在使用以上方程時，每次f函數須要被從新定義，而且每次的返回值必須是一個整數或浮點數。
3. 返回值的精準度會隨着適用的隨機數的數量而變化。使用的數量越大，返回的值將會越精確，反而亦然。在以上的例子裏，99999個隨機數能夠產生三位小數的精確度。

　　完成！

　　Thanks for reading! Have a good day!

 

 

 

 

　　

　　[部分延伸例題選自 University of Toronto, ESC180 課程，見：https://github.com/guerzh/esc180lec]