python實現蒙特卡羅模擬法的實踐

Posted on by

1.簡介

        蒙特卡洛又稱隨機抽樣或統計試驗，就是產生隨機變量，帶入模型算的結果，尋優方面，隻要模擬次數夠多，最終是可以找到最優解或接近最優的解。

         通常蒙特卡羅方法可以粗略地分成兩類：一類是所求解的問題本身具有內在的隨機性，借助計算機的運算能力可以直接模擬這種隨機的過程。例如在核物理研究中，分析中子在反應堆中的傳輸過程。中子與原子核作用受到量子力學規律的制約，人們隻能知道它們相互作用發生的概率，卻無法準確獲得中子與原子核作用時的位置以及裂變產生的新中子的行進速率和方向。科學傢依據其概率進行隨機抽樣得到裂變位置、速度和方向，這樣模擬大量中子的行為後，經過統計就能獲得中子傳輸的范圍，作為反應堆設計的依據。
另一種類型是所求解問題可以轉化為某種隨機分佈的特征數，比如隨機事件出現的概率，或者隨機變量的期望值。通過隨機抽樣的方法，以隨機事件出現的頻率估計其概率，或者以抽樣的數字特征估算隨機變量的數字

2.實例分析

2.1 模擬求近似圓周率

        繪制單位圓和外接正方形，正方形ABCD的面積為：2*2=4，圓的面積為：S=Π*1*1=Π，現在模擬產生在正方形ABCD中均勻分佈的點n個，如果這n個點中有m個點在該圓內，則圓的面積與正方形ABCD的面積之比可近似為m/n

程序如下：

#模擬求近似圓周率
import random
import numpy as np
import matplotlib.pyplot as plt
#解決圖標題中文亂碼問題
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默認字體
mpl.rcParams['axes.unicode_minus'] = False  # 解決保存圖像是負號'-'顯示為方塊的問題
#進入正題
r=random.random()
num=range(0,200000,10)
mypi=np.ones((1,len(num)))
for j in range(len(num)):
    # 投點次數
    n = 10000
    # 圓的半徑、圓心
    r = 1.0
    a,b = (0.,0.)
    # 正方形區域
    x_min, x_max = a-r, a+r
    y_min, y_max = b-r, b+r
    # 在正方形區域內隨機投點
    x = np.random.uniform(x_min, x_max, n) #均勻分佈
    y = np.random.uniform(y_min, y_max, n)
    #計算點到圓心的距離
    d = np.sqrt((x-a)**2 + (y-b)**2)
    #統計落在圓內點的數目
    res = sum(np.where(d < r, 1, 0))
    #計算pi的近似值（Monte Carlo:用統計值去近似真實值）
    mypi[0,j] = 4 * res / n
plt.plot(range(1,len(mypi[0])+1),mypi[0],'.-') 
plt.grid(ls=":",c='b',)#打開坐標網格
plt.axhline(y=np.pi,ls=":",c="yellow")#添加水平直線
# plt.axvline(x=4,ls="-",c="green")#添加垂直直線
plt.legend(['模擬', '實際'], loc='upper right', scatterpoints=1)
plt.title("近似圓周率")
plt.show()

返回：

2.2 估算定積分

程序如下：

#估算定積分
import random
import numpy as np
import matplotlib.pyplot as plt
#解決圖標題中文亂碼問題
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默認字體
mpl.rcParams['axes.unicode_minus'] = False  # 解決保存圖像是負號'-'顯示為方塊的問題
#進入正題
r=random.random()
num=range(1,10**6,500)
s = np.ones((1,len(num)))
for j in range(len(num)):
    n = 30000
    #矩形邊界區域
    x_min, x_max = 0.0, 1.0
    y_min, y_max = 0.0, 1.0
    #在矩形區域內隨機投點x
    x = np.random.uniform(x_min, x_max, n)#均勻分佈
    y = np.random.uniform(y_min, y_max, n)
    #統計落在函數y=x^2下方的點
    res = sum(np.where(y < x**2, 1 ,0))
    #計算定積分的近似值
    s[0,j] = res / n
plt.plot(range(1,len(s[0])+1),s[0],'.-') 
plt.grid(ls=":",c='b',)#打開坐標網格
plt.axhline(y=1/3,ls=":",c="red")#添加水平直線
# plt.axvline(x=4,ls="-",c="green")#添加垂直直線
plt.legend(['模擬', '實際1/3'], loc='upper right', scatterpoints=1)
plt.title("估算定積分")
plt.show()

返回：

2.3 求解整數規劃

要解的方程為：

條件如下：

程序如下：

# 求解整數規劃
import random
import numpy as np
import time
import matplotlib.pyplot as plt
#解決圖標題中文亂碼問題
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默認字體
mpl.rcParams['axes.unicode_minus'] = False  # 解決保存圖像是負號'-'顯示為方塊的問題
#進入正題
time_start=time.time() #計時開始
p0=0
for i in range(10**7):
    x=np.random.randint(0,99,(1,3))
    f=2*x[0,0]+3*x[0,0]**2+3*x[0,1]+x[0,1]**2+x[0,2]
    g=[
        x[0,0]+2*x[0,0]**2+x[0,1]+2*x[0,1]**2+x[0,2],
        x[0,0]+x[0,0]**2+x[0,1]+x[0,1]**2-x[0,2],
        2*x[0,0]+x[0,0]**2+2*x[0,1]+x[0,2],
        x[0,0]+2*x[0,1]**2
    ]
    if g[0]<=100 and g[1]<=500 and g[2]<=400 and g[3]>=10:
        if p0<f:
            x0=x
            p0=f
print('最優解：',x0)
print('最優值：',p0)
time_end=time.time() #計時結束
print('用時：',time_end-time_start)

返回：

 到此這篇關於python實現蒙特卡羅模擬法的實踐的文章就介紹到這瞭,更多相關python 蒙特卡羅模擬法內容請搜索WalkonNet以前的文章或繼續瀏覽下面的相關文章希望大傢以後多多支持WalkonNet！

推薦閱讀：