這篇文章介紹蒙地卡羅法(Monte carlo)求圓周率及橢圓面積，使用Python來實現並將其結果用matlibplot視覺化。
蒙地卡羅法：
不用數值計算方式來解題，而以機率(亂數)來解題的方式，稱為
蒙地卡羅法。
摩洛哥王國是位於法國與義大利之間的小國，該國首都即為蒙地卡羅，以賭博聞名。
而這種以亂數解題所求得的數值不如以數值分析計算所得的結果精確，這種解題方式類似賭博，故稱
蒙地卡羅法。
A.使用蒙地卡羅法求圓周率π：
如下圖Python 程式所畫出的結果來解釋，假設正方形邊長為1，圓半徑也為1，
那麼1/4 圓的面積為 ：
而正方形面積為:1x1=1
那麼產生分布於0~1 之間的亂數，應當會均勻的分布於正方形之內，如下圖綠點所示，
而分布於1/4圓內的數量假設為a ，分布於圓外的數量為b，N則是所產生的亂數總數為
N=a+b。
那麼其亂數分布數量a與N的比值應與1/4圓面積及正方形面積成正比，於是：
即可算出圓周率，如圖一，求得圓周率為3.1268，此種算法沒那麼精確，但是如果將
亂數總數N放大，理當亂數分布會越均勻，則應會越趨近於精確值，如圖二，N=100000的結果，求得圓周率為3.14708
<圖一>使用蒙地卡羅法求圓周率π，當亂數N=10000 ，圓周率π求得3.1268
<圖二>使用蒙地卡羅法求圓周率π，當亂數N=100000 ，圓周率π求得3.14708
B.使用蒙地卡羅法求橢圓面積：
假設一橢圓方程式：
如下圖三，產生亂數0~2之間的實數對應於x軸， 產生亂數0~1之間的實數對應於y軸，其產生的亂數點應分布於2x1的長方形內。如同上述求圓周率π，假設分布於1/4橢圓內的亂數數量為a，1/4橢圓面積為s，全部分布於長方形內的亂數總數為N，那麼根據面積比關係：
<圖三>使用蒙地卡羅法求橢圓面積，當亂數N=100000，橢圓面積求得6.26872
<蒙地卡羅法求圓周率範例程式>
import math
import numpy as np
from matplotlib import pyplot as plt
from decimal import *
def circle(x):
if x <=1:
y = math.sqrt(1.0-math.pow(x,2))
return y
N = 100000
pei_count=0
px = np.random.rand(N)
py = np.random.rand(N)
sx=np.linspace(0.0, 1, 100)
sg=np.array([circle(t) for t in sx])
for i in range(N):
if math.pow(px[i],2)+math.pow(py[i],2)<=1:
pei_count+=1
print("在圓內的點數有",pei_count)
pei=4*pei_count/N
print("圓周率=",pei)
Pei_msg="Pei="+str(pei)
#畫出circle函數
plt.subplot(1,1,1)
plt.plot(sx, sg, 'b-',linewidth=2.0)
plt.axis([0,1.5, 0,1.5]) #分別設定X,Y軸的最小最大值
plt.title(' circle ',fontsize=10)
plt.ylabel('Y')
plt.xlabel('x')
plt.grid(True)
plt.scatter(px, py, c="green", alpha=0.9)
plt.text(0.5, 1.2,Pei_msg)
plt.show()
<蒙地卡羅法求橢圓面積範例程式>
import math
import numpy as np
from matplotlib import pyplot as plt
from decimal import *
def circle(x):
y = math.sqrt(1.0-math.pow(x,2)/4)
return y
N = 100000
Es_count=0
px = np.random.randint(2000, size=N)
pd=np.mat([1000.0]*N)
px=np.divide(px, pd)
px=px.reshape(N,-1)
py = np.random.rand(N)
sx=np.linspace(-2, 2, 100)
sg=np.array([circle(t) for t in sx])
sg2=-sg
for i in range(N):
if (math.pow(px[i],2)/4)+math.pow(py[i],2)<=1:
Es_count+=1
print("在圓內的點數有",Es_count)
Es_area=4*(2.0*Es_count/N)
print("面積=",Es_area)
Area_msg="Ellipse Area="+str(Es_area)
#畫出ellipse函數
plt.subplot(1,1,1)
plt.plot(sx, sg, 'b-',linewidth=2.0)
plt.plot(sx, sg2, 'b-',linewidth=2.0)
plt.axis([-2.0,2.0, -1.5,1.5]) #分別設定X,Y軸的最小最大值
plt.title(' ellipse ',fontsize=10)
plt.ylabel('Y')
plt.xlabel('x')
plt.grid(True)
plt.scatter(px, py, c="green", alpha=0.9)
plt.text(0.5, 1.2,Area_msg)
plt.show()
加入阿布拉機的3D列印與機器人的FB專頁
https://www.facebook.com/arbu00/
<其他有關OPENCV 文章>
機器學習(2)--使用OPENCV SVM實作手寫辨識
機器學習(1)--使用OPENCV KNN實作手寫辨識
OPENCV(11)--contours and convex hull(輪廓與凸包)
OPENCV(10)--Canny Edge Detection(Canny邊緣檢測)
OPENCV(9)--Image Gradients(圖像梯度)
OPENCV(8)--Histogram & Histograms Equalization(長條圖與長條圖均衡化)
OPENCV(7)--2D Convolution ，Image Filtering and Blurring (旋積，濾波與模糊)
OPENCV(6)--Trackbar(軌道桿)
OPENCV(5)--Drawing
OPENCV(4)--Grayscale,Binarization,Threshole(灰階化,二值化,閥值)
OPENCV(3)--Matplotlib pyplot bassic function
OPENCV(2)--Capture Video from Camera
OPENCV(1 )--How to install OPENCV in Python
<其他文章>
人工神經網路(2)--使用Python實作後向傳遞神經網路演算法(Backprogation artificial neature network)
人工神經網路(1)--使用Python實作perceptron(感知器)
