基于python的数学建模—蒙特卡洛算法

import math
import random
m = input('请输入一个较大的整数')
n = 0
for i in range(int(m)):
    x = random.random()
    y = random.random()
    if math.sqrt(x**2 + y**2) < 1:
        n += 1
pi = 4 * n /int(m)
print("pi = {}".format(pi))

请输入一个较大的整数>? 10000000
pi = 3.1425488　

 计算积分

 绘制图像

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,'-')
plt.show()

 计算积分

import random
import numpy as np


m = 100000
n = 0
for i in range(m):
    x = random.random()   # random.random()用于生成一个0到1的随机符点数: 0 <= n < 1.0
    y = random.random()
    if np.log(1 + x) / (1 + x ** 2) > y:
        n += 1
ans = n / m
print(ans)



0.27331

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):
    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])
 
# 统计一个列表中每个元素出现的次数
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()

 三门问题

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)

每次换门的中奖概率：
中奖率为0.667476
每次都不换门的中奖概率：
中奖率为0.333089
为什么两次中将概率相加不等于1   两次不是同时发生的  没有联系　

 M*M豆问题

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