蒙特卡洛模擬是一種統計學方法,基本原理是通過大量的隨機樣本對系統進行模擬,從而求得所需計算的參量。使用蒙特卡洛模擬方法的基本要素包括:構建或描述概率模型、從已知概率分布采樣、建立各種估計量。
使用“簡書-朱煥”的"定量分析項目總持續時間"例子:比如說我們現在有個項目,該項目共有三個WBS要素分別是設計、建造和測試,為了簡單起見我們假設這三個WBS要素的預估的工期概率分布都呈標準正態分布,而且三者之間都是完成到開始的邏輯關系,這樣整個項目工期就是這三個WBS要素工期之和。例子的詳細情況參見[1]。
使用蒙特卡洛模擬方法,首先構建概率模型(這個應該是解決問題的關鍵),這個例子中假設了設計、建造和測試三個要素呈正態分布;然后根據正態分布對這三個要素進行采樣;估計量就是工期時間,工期時間是三個要素之和,根據采樣結果計算工期的模擬值,并計算工期模擬值的出現頻率(概率),最后根據出現頻率計算累積概率。下面是Python代碼:
import numpy as np
import matplotlib.pyplot as plt
import math
# 參數
mu = [14, 23, 22]
sigma = [2, 3, 4]
tips = ['design', 'build', 'test']
figureIndex = 0
fig = plt.figure(figureIndex, figsize=(10,8))
# 顯示分布圖
color = ['r', 'g', 'b']
ax = fig.add_subplot(111)
#ax = plt.subplot(1,1,1)
for i in range(3):
# 參考https://www.jb51.net/article/146073.htm.[2]
x = np.linspace(mu[i] - 3 * sigma[i], mu[i] + 3 * sigma[i], 100)
y_sig = np.exp(-(x - mu[i]) ** 2 / (2 * sigma[i] ** 2)) / (math.sqrt(2 * math.pi) * sigma[i])
ax.plot(x, y_sig, color[i]+'-', linewidth=2, alpha=0.6, label=tips[i])
#
ax.legend(loc='best', frameon=False)
ax.set_xlabel('# of days')
ax.set_ylabel('probability')
plt.grid(True)
# 蒙特卡洛采樣
# 三個WBS要素
size = 10000
samples = [np.random.normal(mu[i], sigma[i], size) for i in range(3)]
# 計算工期
data = np.zeros(len(samples[1]))
for i in range(len(samples[1])):
for j in range(3):
data[i] += samples[j][i]
data[i] = int(data[i])
# 統計一個列表中每個元素出現的次數
# 參考https://blog.csdn.net/qq_42467563/article/details/86182266.[3]
def count(lis):
lis=np.array(lis)
key=np.unique(lis)
x = []
y = []
for k in key:
mask =(lis == k)
list_new=lis[mask]
v=list_new.size
x.append(k)
y.append(v)
return x,y
#
# 計算工期出現頻率與累積概率
a,b = count(data)
pdf = [x/size for x in b]
cdf = np.zeros(len(a))
for i in range(len(a)):
if i > 0:
cdf[i] += cdf[i-1]
cdf[i] += b[i]
cdf = cdf/size
figureIndex += 1
fig = plt.figure(figureIndex, figsize=(10,8))
ax = fig.add_subplot(211)
ax.bar(a, height=pdf, color = 'blue',edgecolor = 'white', label='MC PDF')
ax.plot(a, pdf)
ax.legend(loc='best', frameon=False)
ax.set_xlabel('# of days for project')
ax.set_ylabel('probability')
ax.set_title('Monte Carlo Simulation')
ax = fig.add_subplot(212)
ax.plot(a, cdf, 'r-', marker='o', mfc='b', ms=4, lw=2, alpha=0.6, label='MC CDF')
ax.legend(loc='best', frameon=False)
ax.set_xlabel('# of days for project')
ax.set_ylabel('probability')
ax.grid(True)
plt.show()
模擬結果:
參考資料:
[1] 朱煥, “蒙特卡洛模擬(Monte Carlo Simulation)淺析”, https://www.jianshu.com/p/cb44f4b457c3.
[2] 腳本之家, “Python使用numpy產生正態分布隨機數的向量或矩陣操作示例”, https://www.jb51.net/article/146073.htm.
[3] 三尺秋水一點飛鴻, “Python統計一個列表中每個元素出現的次數。四種方法,總有一款適合你”, https://blog.csdn.net/qq_42467563/article/details/86182266.
更多文章、技術交流、商務合作、聯系博主
微信掃碼或搜索:z360901061

微信掃一掃加我為好友
QQ號聯系: 360901061
您的支持是博主寫作最大的動力,如果您喜歡我的文章,感覺我的文章對您有幫助,請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧,狠狠點擊下面給點支持吧,站長非常感激您!手機微信長按不能支付解決辦法:請將微信支付二維碼保存到相冊,切換到微信,然后點擊微信右上角掃一掃功能,選擇支付二維碼完成支付。
【本文對您有幫助就好】元
