IV.python初探日記：python實現蒙特卡洛方法計算π值

完整的代碼在文末，本文適合小白看，完全白紙的都可以，也希望大神們不吝賜教。

一、最簡單的實現方法

下面是最簡單的實現方式，模擬試驗一千萬次，但模擬出來的π值並不精確。

import randomzongshu = 10000000jishu = 0for i in range(zongshu) : x = random.random() y = random.random() if (x ** 2 + y ** 2) < 1 : jishu+=1print(jishu * 4.0 / zongshu)

二、關於代碼的效率問題

上面這段代碼的計算速度和精度都相當粗糙，然而我還沒有對改進代碼的想法，所以下面就直接暴力模擬了。下面是粗略地試驗了一下各種實驗次數下的實現時間，此處引用的是我修改後的代碼，原代碼為——使用Python語言的蒙特卡洛方法計算圓周率π的一種實現。

from __future__ import divisionimport randomimport timenum = 1for j in range(1, 7): startT = time.clock() # 落入圓內計數 counter = 0 # 往正方形中扔了 10 的 j 次方次 for i in range(10 ** j): x = random.uniform(-1, 1) y = random.uniform(-1, 1) # 落入了圓內 if x**2 + y**2 < 1: counter = counter + 1 endT = time.clock() print ( num) print ( pi:{0}.format(4 * (counter / 10 ** j))) print ( Time:{0}.format(endT - startT)) print () num += 1

模擬試驗一千萬次以上的時間就相當漫長了。

三、matplotlib輸出

左圖輸出的是模擬的1/4圓試驗，右圖輸出的是試驗次數和計算π值的關係，這裡主要碰到了下面這幾個問題。

1、生成隨機點

主要的思想是生成兩組隨機數分別加入兩組數組。原本不知道matplotlib是可以直接輸出數組的，折騰了好久才發現這麼簡單......嘗試過的方法有：循環「生成一個隨機點」這個步驟；生成兩對隨機數組，然後一一匹配得到隨機點；還有一些奇葩的就不列舉了......總之這些步驟都沒有去搜索，純粹想鍛煉一下自己的思維。

2、1/4圓「那條線」

一開始是想用matplotlib直接畫一個圓，結果發現matplotlib不能直接畫圓（也有可能是我沒找到，知道的同學告訴我一聲）；隨後想到圓環、實心圓塗色等幾個方法，但實現起來都較為複雜，最後我找到了這樣一個方法numpy.meshgrid，直接把這個1/4圓描了出來：

x = y = np.arange(-1, 1, 0.001)x, y = np.meshgrid(x,y)p1.contour(x, y, x**2 + y**2, [1])

3、右圖代碼的取捨

右邊的關係圖，考慮到速度和精度，取的是1000。

C = [ (j+1) for j in range(1000)]D = [ M((j+1)) for j in range(1000)]

其實可以改進為下面這樣，但我原代碼中沒有修改。

C = [ 100 * (j+1) for j in range(100)]D = [ M(100 * (j+1)) for j in range(100)]

4、排版問題

用的是「121」的排版方式，出來的圖手工拉動了一下大小，不過似乎matplotlib是可以直接控制輸出的，這點我要再琢磨下。

5、遺留問題

左圖其實可以直接生成隨機點；「1/4圓」的問題沒有從根本上解決；輸出的排版沒有解決；代碼的優化完全沒有涉及；另外比較現實的是，如何把「1/4圓」一下的點抹成紅色，「1/4圓」以上的點保持藍色，這點我還沒有去嘗試解決。上面的遺留問題如果有進一步的探索的話，我會在這裡跟進的。

四、完整代碼

最後是完整的代碼。

import randomimport matplotlib.pyplot as pltimport numpy as npA = []B = []cishu = 100for i in range(cishu): X = random.random() A.append(X) Y = random.random() B.append(Y)def M(zongshu): jishu = 0.0 for j in range(zongshu): x = random.random() y = random.random() if ( x**2 + y**2 ) < 1: jishu += 1 return (jishu/zongshu * 4.0)fig=plt.figure()p1=fig.add_subplot(121)p1.axis([0,1,0,1])p1.scatter(A, B) x = y = np.arange(-1, 1, 0.001)x, y = np.meshgrid(x,y)p1.contour(x, y, x**2 + y**2, [1])p2=fig.add_subplot(122)C = [ (j+1) for j in range(1000)]D = [ M((j+1)) for j in range(1000)]p2.plot(C, D)plt.show()

推薦閱讀：