蒙特卡洛算法

什么是蒙特卡洛算法？

蒙特·卡罗（ Monte Carlo method），又称统计模拟方法，是二十世纪四十年代中叶由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随

基本思想：

当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。

 

样例1:圆周率的计算

基本思想：在图中区域产生足够多的随机数点，然后计算落在圆内的点的个数与总个数的比值再乘以4，就是圆周率。

import math
import random
m = 10000
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))

运行结果：

 样例2
利用python计算函数X的平方在闭区间0到1上的积分

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))

运行结果：

 样例3：

 先绘制函数图像：

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,'-')

 计算积分：

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
print(ans)

运行结果：

样例4

现在有个项目，共三个WBS要素，分别是设计、建造和测试。假设这三个WBS要素预估工期的概率分布呈标准正态分布，且三者之间都是完成到开始的逻辑关系，于是整个项目工期就是三个WBS要素工期之和。

 

首先，先画出这三个要素的概率密度函数

从表中可以看出 时间最少为8天 最长为34天 所以我们设置时间范围为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()

运行结果：

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 = np.mean(y)
b = np.var(y)
print(a,b)

 运行结果：

2、三门问题 

三门问题（Monty Hall probelm）亦称为蒙提霍尔问题，出自美国电视游戏节目Let’s Make a Deal。
参赛者会看见三扇关闭了的门，其中一扇的后面有一辆汽车，选中后面有车的那扇门可赢得该汽车，另外两扇门则各藏有一只山羊。
当参赛者选定了一扇门，但未去开启它的时候，节目主持人开启剩下两扇门的其中一扇，露出其中一只山羊。主持人其后问参赛者要不要换另一扇仍然关上的门。
问题是：换另一扇门是否会增加参赛者赢得汽车的几率？如果严格按照上述条件，即主持人清楚地知道，自己打开的那扇门后面是羊，那么答案是会。不换门的话，赢得汽车的几率是1/3,，换门的话，赢得汽车的几率是2/3。
 

import math
import  random
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)

运行结果：

 3、M*M豆M&M豆是有各种颜色的糖果巧克力豆。制造M&M豆的Mars公司会不时变更不同颜色巧克力豆之间的混合比例。
1995年，他们推出了蓝色的M&M豆。在此前一袋普通的M&M豆中，颜色的搭配为：30%褐色，20%黄色，20%红色，10%绿色，10%橙色，10%黄褐色。这之后变成了：24%蓝色，20%绿色，16%橙色，14%黄色，13%红色，13%褐色。
假设我的一个朋友有两袋M&M豆，他告诉我一袋是1994年，一带是1996年。
但他没告诉我具体那个袋子是哪一年的，他从每个袋子里各取了一个M&M豆给我。一个是黄色，一个是绿色的。那么黄色豆来自1994年的袋子的概率是多少？
 

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: {'蓝色': 24, '绿色': 20, '橙色': 16, '黄色': 14, '红色': 13, '褐色': 13}}
    num = 10000
    list_1994 = ['褐色'] * 30 * num + ['黄色'] * 20 * num + ['红色'] * 20 * num + ['绿色'] * 10 * num + [
        '橙色'] * 10 * num + ['黄褐'] * 10 * num
    list_1996 = ['蓝色'] * 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()))

运行结果：