Python數(shù)學(xué)建模系列（六）：蒙特卡洛算法
1、蒙特卡洛算法
什么是蒙特卡洛算法？
?蒙特·卡羅（ Monte Carlo method），又稱統(tǒng)計模擬方法，是二十世紀(jì)四十年代中葉由于科學(xué)技術(shù)的發(fā)展和電子計算機(jī)的發(fā)明，而被提出的一種以概率統(tǒng)計理論為指導(dǎo)的一類非常重要的數(shù)值計算方法。是指使用隨機(jī)數(shù)（或更常見的偽隨機(jī)數(shù)）來解決很多計算問題的方法。
?
起源：
?蒙特卡洛方法于20世紀(jì)40年代美國在第二次世界大戰(zhàn)中研制原子彈的“曼哈頓計劃”計劃成員S.M.烏拉姆和J.馮·諾依曼首先提出。數(shù)學(xué)家馮·諾依曼用馳名世界的賭城——摩納哥的Monte Carlo——來命名這種方法，為它蒙上一層神秘色彩。公認(rèn)的蒙特卡洛方法的起源是1777年法國數(shù)學(xué)家布豐提出用投針實驗的方法求圓周率π。
?
基本思想：
?當(dāng)所求解問題是某種隨機(jī)事件出現(xiàn)的概率，或者是某個隨機(jī)變量的期望值時，通過某種“實驗”的方法，以這種事件出現(xiàn)的頻率估計這一隨機(jī)事件的概率，或者得到這個隨機(jī)變量的某些數(shù)字特征，并將其作為問題的解。
?
工作步驟：
構(gòu)造或描述概率過程 實現(xiàn)從已知概率分布抽樣 建立各種估計量
樣例1
題目：經(jīng)典的蒙特卡洛方法求圓周率
?基本思想：在圖中區(qū)域產(chǎn)生足夠多的隨機(jī)數(shù)點，然后計算落在圓內(nèi)的點的個數(shù)與總個數(shù)的比值再乘以4，就是圓周率。
?
Demo代碼 「（模擬次數(shù)為10000次時）」
import math
import random
m = 10000
n = 0
for i in range(m):
# x、y為0-1之間的隨機(jī)數(shù)
x = random.random()
y = random.random()
# 若點(x,y) 屬于圖中1/4圓內(nèi) 則有效個數(shù)+1
if math.sqrt(x**2 + y**2) < 1:
n += 1
# 計算pi
pi = 4 * n / m
print("pi = {}".format(pi))
# pi = 3.1508（結(jié)果具有隨機(jī)性 不一定完全一樣）
Demo代碼 「（模擬次數(shù)為1000000次時）」
import math
import random
m = 1000000
n = 0
for i in range(m):
x = random.random()
y = random.random()
if math.sqrt(x**2 + y**2) < 1:
n += 1
pi = 4 * n / m
print("pi = {}".format(pi))
# pi = 3.141416
樣例2
利用python計算函數(shù)在[0,1]區(qū)間的定積分
Demo代碼
import math
import random
m = 1000000
n = 0
for i in range(m):
x = random.random()
y = random.random()
if y >= x**2:
n += 1
r = n / m
print("r = {}".format(r))
# r = 0.666817
變式：計算積分
?題目來源：https://blog.csdn.net/FightingBob/article/details/109694130
理論答案：
?
先繪制函數(shù)圖像：
import numpy as np
import matplotlib.pylab as plt
x = np.linspace(0,1,num=50)
y = np.log(1 + x) / (1 + x**2)
plt.plot(x,y,'-')
函數(shù)圖像如下：
計算積分
import numpy as np
import random
m = 100000
n = 0
for i in range(m):
x = random.random()
y = random.random()
if np.log(1 + x) / (1 + x**2) > y:
n += 1
ans = n / m
ans
# 0.27293
# 理論答案是 pi/8*log(2) = 0.2721983
樣例3
現(xiàn)在有個項目，共三個WBS要素，分別是設(shè)計、建造和測試。假設(shè)這三個WBS要素預(yù)估工期的概率分布呈標(biāo)準(zhǔn)正態(tài)分布，且三者之間都是完成到開始的邏輯關(guān)系，于是整個項目工期就是三個WBS要素工期之和。
首先，先畫出這三個要素的概率密度函數(shù)
從表中可以看出 時間最少為8天 最長為34天 所以我們設(shè)置時間范圍為7-35 天
import numpy as np
import matplotlib.pyplot as plt
import math
def normal_distribution(x, mean, sigma):
return np.exp(-1*((x-mean)**2)/(2*(sigma**2)))/(math.sqrt(2*np.pi) * sigma)
mean1, sigma1 = 14, 2
x1 = np.linspace(7, 35,num=100)
mean2, sigma2 = 23, 3
x2 = np.linspace(7, 35, num=100)
mean3, sigma3 =22, 4
x3 = np.linspace(7, 35, num=100)
y1 = normal_distribution(x1, mean1, sigma1)
y2 = normal_distribution(x2, mean2, sigma2)
y3 = normal_distribution(x3, mean3, sigma3)
plt.plot(x1, y1, 'r', label='m=14,sig=2')
plt.plot(x2, y2, 'g', label='m=23,sig=3')
plt.plot(x3, y3, 'b', label='m=22,sig=4')
plt.legend()
plt.grid()
plt.show()
圖像如下：
Demo代碼
import numpy as np
import matplotlib.pylab as plt
import random
import math
m = 10000
a = 0
b = 0
y1 = np.random.normal(loc=14,scale=2,size=m)
y2 = np.random.normal(loc=23,scale=3,size=m)
y3 = np.random.normal(loc=22,scale=4,size=m)
y =y1 + y2 + y3
a,b = np.mean(y),np.var(y)
a,b
# (59.08396232409535, 29.120136564111085)
# 理論均值為59 = 14 + 23 + 22 方差為29 = 2*2 + 3*3 + 4*4
2、三門問題
三門問題
?三門問題（Monty Hall probelm）亦稱為蒙提霍爾問題，出自美國電視游戲節(jié)目Let’s Make a Deal。
參賽者會看見三扇關(guān)閉了的門，其中一扇的后面有一輛汽車，選中后面有車的那扇門可贏得該汽車，另外兩扇門則各藏有一只山羊。
當(dāng)參賽者選定了一扇門，但未去開啟它的時候，節(jié)目主持人開啟剩下兩扇門的其中一扇，露出其中一只山羊。主持人其后問參賽者要不要換另一扇仍然關(guān)上的門。
問題是：換另一扇門是否會增加參賽者贏得汽車的幾率？如果嚴(yán)格按照上述條件，即主持人清楚地知道，自己打開的那扇門后面是羊，那么答案是會。不換門的話，贏得汽車的幾率是1/3,，換門的話，贏得汽車的幾率是2/3。
?
蒙特卡洛思想的應(yīng)用
?應(yīng)用蒙特卡洛重點在使用隨機(jī)數(shù)來模擬類似于賭博問題的贏率問題，并且通過多次模擬得到所要計算值的模擬值。
在三門問題中，用0、1、2分代表三扇門的編號，在[0,2]之間隨機(jī)生成一個整數(shù)代表獎品所在門的編號prize，再次在[0,2]之間隨機(jī)生成一個整數(shù)代表參賽者所選擇的門的編號guess。用變量change代表游戲中的換門（true）與不換門（false）。
?
蒙特卡洛過程
Demo代碼
import math
def play(change):
prize = random.randint(0,2)
guess = random.randint(0,2)
if guess == prize:
if change:
return False
else:
return True
else:
if change:
return True
else:
return False
def winRate(change, N):
win = 0
for i in range(N):
if(play(change)):
win += 1
print("中獎率為{}".format(win / N))
N = 1000000
print("每次換門的中獎概率：")
winRate(True,N)
print("每次都不換門的中獎概率：")
winRate(False,N)
# 理論換門2/3 不換門1/3
#每次換門的中獎概率：
#中獎率為0.665769
#每次都不換門的中獎概率：
#中獎率為0.33292
運行結(jié)果
3、M*M豆問題
M*M豆貝葉斯統(tǒng)計問題
?M&M豆是有各種顏色的糖果巧克力豆。制造M&M豆的Mars公司會不時變更不同顏色巧克力豆之間的混合比例。
1995年，他們推出了藍(lán)色的M&M豆。在此前一袋普通的M&M豆中，顏色的搭配為：30%褐色，20%黃色，20%紅色，10%綠色，10%橙色，10%黃褐色。這之后變成了：24%藍(lán)色，20%綠色，16%橙色，14%黃色，13%紅色，13%褐色。
假設(shè)我的一個朋友有兩袋M&M豆，他告訴我一袋是1994年，一帶是1996年。
但他沒告訴我具體那個袋子是哪一年的，他從每個袋子里各取了一個M&M豆給我。一個是黃色，一個是綠色的。那么黃色豆來自1994年的袋子的概率是多少？
?
Demo代碼
import time
import random
for i in range(10):
print(time.strftime("%Y-%m-%d %X",time.localtime()))
dou = {1994:{'褐色':30,'黃色':20,'紅色':20,'綠色':10,'橙色':10,'黃褐':30},
1996:{'藍(lán)色':24,'綠色':20,'橙色':16,'黃色':14,'紅色':13,'褐色':13}}
num = 10000
list_1994 = ['褐色']*30*num+['黃色']*20*num+['紅色']*20*num+['綠色']*10*num+['橙色']*10*num+['黃褐']*10*num
list_1996 = ['藍(lán)色']*24*num+['綠色']*20*num+['橙色']*16*num+['黃色']*14*num+['紅色']*13*num+['褐色']*13*num
random.shuffle(list_1994)
random.shuffle(list_1996)
count_all = 0
count_key = 0
for key in range(100 * num):
if list_1994[key] == '黃色' and list_1996[key] == '綠色':
count_all += 1
count_key += 1
if list_1994[key] == '綠色' and list_1996[key] == '黃色':
count_all += 1
print(count_key / count_all,20/27)
print(time.strftime("%Y-%m-%d %X",time.localtime()))
# ...
# 0.7407064573459715 0.7407407407407407
# 20/27是理論答案
運行結(jié)果
結(jié)語
參考：
https://www.bilibili.com/video/BV12h411d7Dm https://blog.csdn.net/wangyuhang124/article/details/114241080
學(xué)習(xí)來源：B站及其課堂PPT，對其中代碼進(jìn)行了復(fù)現(xiàn)
「文章僅作為學(xué)習(xí)筆記，記錄從0到1的一個過程」
希望對您有所幫助，如有錯誤歡迎小伙伴指正～