蒙特卡洛模擬是一種統(tǒng)計(jì)學(xué)方法，基本原理是通過大量的隨機(jī)樣本對系統(tǒng)進(jìn)行模擬，從而求得所需計(jì)算的參量。使用蒙特卡洛模擬方法的基本要素包括：構(gòu)建或描述概率模型、從已知概率分布采樣、建立各種估計(jì)量。

使用“簡書-朱煥”的"定量分析項(xiàng)目總持續(xù)時間"例子：比如說我們現(xiàn)在有個項(xiàng)目，該項(xiàng)目共有三個WBS要素分別是設(shè)計(jì)、建造和測試，為了簡單起見我們假設(shè)這三個WBS要素的預(yù)估的工期概率分布都呈標(biāo)準(zhǔn)正態(tài)分布，而且三者之間都是完成到開始的邏輯關(guān)系，這樣整個項(xiàng)目工期就是這三個WBS要素工期之和。例子的詳細(xì)情況參見[1]。

使用蒙特卡洛模擬方法，首先構(gòu)建概率模型(這個應(yīng)該是解決問題的關(guān)鍵)，這個例子中假設(shè)了設(shè)計(jì)、建造和測試三個要素呈正態(tài)分布；然后根據(jù)正態(tài)分布對這三個要素進(jìn)行采樣；估計(jì)量就是工期時間，工期時間是三個要素之和，根據(jù)采樣結(jié)果計(jì)算工期的模擬值，并計(jì)算工期模擬值的出現(xiàn)頻率(概率)，最后根據(jù)出現(xiàn)頻率計(jì)算累積概率。下面是Python代碼：

            
              import numpy as np
import matplotlib.pyplot as plt
import math

# 參數(shù)
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)]
# 計(jì)算工期
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])

# 統(tǒng)計(jì)一個列表中每個元素出現(xiàn)的次數(shù)
# 參考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
#

# 計(jì)算工期出現(xiàn)頻率與累積概率
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()
            
          

模擬結(jié)果：

參考資料：

[1] 朱煥, “蒙特卡洛模擬(Monte Carlo Simulation)淺析”, https://www.jianshu.com/p/cb44f4b457c3.

[2] 腳本之家, “Python使用numpy產(chǎn)生正態(tài)分布隨機(jī)數(shù)的向量或矩陣操作示例”, https://www.jb51.net/article/146073.htm.

[3] 三尺秋水一點(diǎn)飛鴻, “Python統(tǒng)計(jì)一個列表中每個元素出現(xiàn)的次數(shù)。四種方法，總有一款適合你”, https://blog.csdn.net/qq_42467563/article/details/86182266.