python實(shí)現(xiàn)蒙特卡羅模擬法的實(shí)踐
1.簡(jiǎn)介
蒙特卡洛又稱隨機(jī)抽樣或統(tǒng)計(jì)試驗(yàn),就是產(chǎn)生隨機(jī)變量,帶入模型算的結(jié)果,尋優(yōu)方面,只要模擬次數(shù)夠多,最終是可以找到最優(yōu)解或接近最優(yōu)的解。
通常蒙特卡羅方法可以粗略地分成兩類:一類是所求解的問題本身具有內(nèi)在的隨機(jī)性,借助計(jì)算機(jī)的運(yùn)算能力可以直接模擬這種隨機(jī)的過程。例如在核物理研究中,分析中子在反應(yīng)堆中的傳輸過程。中子與原子核作用受到量子力學(xué)規(guī)律的制約,人們只能知道它們相互作用發(fā)生的概率,卻無法準(zhǔn)確獲得中子與原子核作用時(shí)的位置以及裂變產(chǎn)生的新中子的行進(jìn)速率和方向??茖W(xué)家依據(jù)其概率進(jìn)行隨機(jī)抽樣得到裂變位置、速度和方向,這樣模擬大量中子的行為后,經(jīng)過統(tǒng)計(jì)就能獲得中子傳輸?shù)姆秶?,作為反?yīng)堆設(shè)計(jì)的依據(jù)。
另一種類型是所求解問題可以轉(zhuǎn)化為某種隨機(jī)分布的特征數(shù),比如隨機(jī)事件出現(xiàn)的概率,或者隨機(jī)變量的期望值。通過隨機(jī)抽樣的方法,以隨機(jī)事件出現(xiàn)的頻率估計(jì)其概率,或者以抽樣的數(shù)字特征估算隨機(jī)變量的數(shù)字
2.實(shí)例分析
2.1 模擬求近似圓周率
繪制單位圓和外接正方形,正方形ABCD的面積為:2*2=4,圓的面積為:S=Π*1*1=Π,現(xiàn)在模擬產(chǎn)生在正方形ABCD中均勻分布的點(diǎn)n個(gè),如果這n個(gè)點(diǎn)中有m個(gè)點(diǎn)在該圓內(nèi),則圓的面積與正方形ABCD的面積之比可近似為m/n
程序如下:
#模擬求近似圓周率 import random import numpy as np import matplotlib.pyplot as plt #解決圖標(biāo)題中文亂碼問題 import matplotlib as mpl mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認(rèn)字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負(fù)號(hào)'-'顯示為方塊的問題 #進(jìn)入正題 r=random.random() num=range(0,200000,10) mypi=np.ones((1,len(num))) for j in range(len(num)): # 投點(diǎn)次數(shù) n = 10000 # 圓的半徑、圓心 r = 1.0 a,b = (0.,0.) # 正方形區(qū)域 x_min, x_max = a-r, a+r y_min, y_max = b-r, b+r # 在正方形區(qū)域內(nèi)隨機(jī)投點(diǎn) x = np.random.uniform(x_min, x_max, n) #均勻分布 y = np.random.uniform(y_min, y_max, n) #計(jì)算點(diǎn)到圓心的距離 d = np.sqrt((x-a)**2 + (y-b)**2) #統(tǒng)計(jì)落在圓內(nèi)點(diǎn)的數(shù)目 res = sum(np.where(d < r, 1, 0)) #計(jì)算pi的近似值(Monte Carlo:用統(tǒng)計(jì)值去近似真實(shí)值) mypi[0,j] = 4 * res / n plt.plot(range(1,len(mypi[0])+1),mypi[0],'.-') plt.grid(ls=":",c='b',)#打開坐標(biāo)網(wǎng)格 plt.axhline(y=np.pi,ls=":",c="yellow")#添加水平直線 # plt.axvline(x=4,ls="-",c="green")#添加垂直直線 plt.legend(['模擬', '實(shí)際'], loc='upper right', scatterpoints=1) plt.title("近似圓周率") plt.show()
返回:
2.2 估算定積分
程序如下:
#估算定積分 import random import numpy as np import matplotlib.pyplot as plt #解決圖標(biāo)題中文亂碼問題 import matplotlib as mpl mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認(rèn)字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負(fù)號(hào)'-'顯示為方塊的問題 #進(jìn)入正題 r=random.random() num=range(1,10**6,500) s = np.ones((1,len(num))) for j in range(len(num)): n = 30000 #矩形邊界區(qū)域 x_min, x_max = 0.0, 1.0 y_min, y_max = 0.0, 1.0 #在矩形區(qū)域內(nèi)隨機(jī)投點(diǎn)x x = np.random.uniform(x_min, x_max, n)#均勻分布 y = np.random.uniform(y_min, y_max, n) #統(tǒng)計(jì)落在函數(shù)y=x^2下方的點(diǎn) res = sum(np.where(y < x**2, 1 ,0)) #計(jì)算定積分的近似值 s[0,j] = res / n plt.plot(range(1,len(s[0])+1),s[0],'.-') plt.grid(ls=":",c='b',)#打開坐標(biāo)網(wǎng)格 plt.axhline(y=1/3,ls=":",c="red")#添加水平直線 # plt.axvline(x=4,ls="-",c="green")#添加垂直直線 plt.legend(['模擬', '實(shí)際1/3'], loc='upper right', scatterpoints=1) plt.title("估算定積分") plt.show()
返回:
2.3 求解整數(shù)規(guī)劃
要解的方程為:
條件如下:
程序如下:
# 求解整數(shù)規(guī)劃 import random import numpy as np import time import matplotlib.pyplot as plt #解決圖標(biāo)題中文亂碼問題 import matplotlib as mpl mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默認(rèn)字體 mpl.rcParams['axes.unicode_minus'] = False # 解決保存圖像是負(fù)號(hào)'-'顯示為方塊的問題 #進(jìn)入正題 time_start=time.time() #計(jì)時(shí)開始 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('最優(yōu)解:',x0) print('最優(yōu)值:',p0) time_end=time.time() #計(jì)時(shí)結(jié)束 print('用時(shí):',time_end-time_start)
返回:
到此這篇關(guān)于python實(shí)現(xiàn)蒙特卡羅模擬法的實(shí)踐的文章就介紹到這了,更多相關(guān)python 蒙特卡羅模擬法內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家!
相關(guān)文章
向量化操作改進(jìn)數(shù)據(jù)分析工作流的Pandas?Numpy示例分析
這篇文章主要介紹了向量化操作改進(jìn)數(shù)據(jù)分析工作流的Pandas?Numpy示例分析,有需要的朋友可以借鑒參考下,希望能夠有所幫助,祝大家多多進(jìn)步,早日升職加薪2023-10-10Python繪制分段函數(shù)的實(shí)現(xiàn)示例
本文主要介紹了Python繪制分段函數(shù)的實(shí)現(xiàn)示例,文中通過示例代碼介紹的非常詳細(xì),對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2023-04-04詳解python的四種內(nèi)置數(shù)據(jù)結(jié)構(gòu)
這篇文章主要介紹了python的四種內(nèi)置數(shù)據(jù)結(jié)構(gòu),文中通過示例代碼介紹的非常詳細(xì),對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2019-03-03Python3之亂碼\xe6\x97\xa0\xe6\xb3\x95處理方式
這篇文章主要介紹了Python3之亂碼\xe6\x97\xa0\xe6\xb3\x95處理方式,具有很好的參考價(jià)值,希望對(duì)大家有所幫助。一起跟隨小編過來看看吧2020-05-05Python特征降維知識(shí)點(diǎn)總結(jié)
在本篇文章里小編給大家整理了一篇關(guān)于Python特征降維知識(shí)點(diǎn)總結(jié)內(nèi)容,有需要的朋友們可以學(xué)習(xí)參考下。2021-08-08Python實(shí)現(xiàn)簡(jiǎn)單石頭剪刀布游戲
這篇文章主要為大家詳細(xì)介紹了Python實(shí)現(xiàn)簡(jiǎn)單的石頭剪刀布的游戲,文中示例代碼介紹的非常詳細(xì),具有一定的參考價(jià)值,感興趣的小伙伴們可以參考一下2019-01-01python爬蟲beautifulsoup庫使用操作教程全解(python爬蟲基礎(chǔ)入門)
這篇文章主要介紹了python爬蟲beautifulsoup庫使用操作全解(python爬蟲基礎(chǔ)入門),本文給大家介紹的非常詳細(xì),對(duì)大家的學(xué)習(xí)或工作具有一定的參考借鑒價(jià)值,需要的朋友可以參考下2021-02-02