python实现蒙特卡罗模拟法的实践
作者:洋洋菜鸟 发布时间:2023-08-11 02:22:02
标签:python,蒙特卡罗模拟法
1.简介
蒙特卡洛又称随机抽样或统计试验,就是产生随机变量,带入模型算的结果,寻优方面,只要模拟次数够多,最终是可以找到最优解或接近最优的解。
通常蒙特卡罗方法可以粗略地分成两类:一类是所求解的问题本身具有内在的随机性,借助计算机的运算能力可以直接模拟这种随机的过程。例如在核物理研究中,分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约,人们只能知道它们相互作用发生的概率,却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。科学家依据其概率进行随机抽样得到裂变位置、速度和方向,这样模拟大量中子的行为后,经过统计就能获得中子传输的范围,作为反应堆设计的依据。
另一种类型是所求解问题可以转化为某种随机分布的特征数,比如随机事件出现的概率,或者随机变量的期望值。通过随机抽样的方法,以随机事件出现的频率估计其概率,或者以抽样的数字特征估算随机变量的数字
2.实例分析
2.1 模拟求近似圆周率
绘制单位圆和外接正方形,正方形ABCD的面积为:2*2=4,圆的面积为:S=Π*1*1=Π,现在模拟产生在正方形ABCD中均匀分布的点n个,如果这n个点中有m个点在该圆内,则圆的面积与正方形ABCD的面积之比可近似为m/n
程序如下:
#模拟求近似圆周率
import random
import numpy as np
import matplotlib.pyplot as plt
#解决图标题中文乱码问题
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#进入正题
r=random.random()
num=range(0,200000,10)
mypi=np.ones((1,len(num)))
for j in range(len(num)):
# 投点次数
n = 10000
# 圆的半径、圆心
r = 1.0
a,b = (0.,0.)
# 正方形区域
x_min, x_max = a-r, a+r
y_min, y_max = b-r, b+r
# 在正方形区域内随机投点
x = np.random.uniform(x_min, x_max, n) #均匀分布
y = np.random.uniform(y_min, y_max, n)
#计算点到圆心的距离
d = np.sqrt((x-a)**2 + (y-b)**2)
#统计落在圆内点的数目
res = sum(np.where(d < r, 1, 0))
#计算pi的近似值(Monte Carlo:用统计值去近似真实值)
mypi[0,j] = 4 * res / n
plt.plot(range(1,len(mypi[0])+1),mypi[0],'.-')
plt.grid(ls=":",c='b',)#打开坐标网格
plt.axhline(y=np.pi,ls=":",c="yellow")#添加水平直线
# plt.axvline(x=4,ls="-",c="green")#添加垂直直线
plt.legend(['模拟', '实际'], loc='upper right', scatterpoints=1)
plt.title("近似圆周率")
plt.show()
返回:
2.2 估算定积分
程序如下:
#估算定积分
import random
import numpy as np
import matplotlib.pyplot as plt
#解决图标题中文乱码问题
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#进入正题
r=random.random()
num=range(1,10**6,500)
s = np.ones((1,len(num)))
for j in range(len(num)):
n = 30000
#矩形边界区域
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0
#在矩形区域内随机投点x
x = np.random.uniform(x_min, x_max, n)#均匀分布
y = np.random.uniform(y_min, y_max, n)
#统计落在函数y=x^2下方的点
res = sum(np.where(y < x**2, 1 ,0))
#计算定积分的近似值
s[0,j] = res / n
plt.plot(range(1,len(s[0])+1),s[0],'.-')
plt.grid(ls=":",c='b',)#打开坐标网格
plt.axhline(y=1/3,ls=":",c="red")#添加水平直线
# plt.axvline(x=4,ls="-",c="green")#添加垂直直线
plt.legend(['模拟', '实际1/3'], loc='upper right', scatterpoints=1)
plt.title("估算定积分")
plt.show()
返回:
2.3 求解整数规划
要解的方程为:
条件如下:
程序如下:
# 求解整数规划
import random
import numpy as np
import time
import matplotlib.pyplot as plt
#解决图标题中文乱码问题
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#进入正题
time_start=time.time() #计时开始
p0=0
for i in range(10**7):
x=np.random.randint(0,99,(1,3))
f=2*x[0,0]+3*x[0,0]**2+3*x[0,1]+x[0,1]**2+x[0,2]
g=[
x[0,0]+2*x[0,0]**2+x[0,1]+2*x[0,1]**2+x[0,2],
x[0,0]+x[0,0]**2+x[0,1]+x[0,1]**2-x[0,2],
2*x[0,0]+x[0,0]**2+2*x[0,1]+x[0,2],
x[0,0]+2*x[0,1]**2
]
if g[0]<=100 and g[1]<=500 and g[2]<=400 and g[3]>=10:
if p0<f:
x0=x
p0=f
print('最优解:',x0)
print('最优值:',p0)
time_end=time.time() #计时结束
print('用时:',time_end-time_start)
返回:
来源:https://blog.csdn.net/qq_25990967/article/details/122878392
0
投稿
猜你喜欢
- phpMyAdmin 是一套可以通过WEB来管理 MySQL-server 以及单一数据库的 PHP
- 你和用户之间的网站堆栈(简化版)在TXJS大会的最后一天,一个开发者问我:面向对象的CSS没有给你留下一大堆基于表现的class名?网络堆栈
- 前言以前在浏览文章时,看到过一个Android的加载效果,觉得挺好看的,于是自己就模仿了一个。下面话不多说,我们直接来看看详细的介绍吧。运行
- 随着jQuery、Mootools、prototype等知名的JavaScript框架的应用变的越来越强大,浏览器对最新版本CSS属性的支持
- 从前有三只小猪,长大自立了分别造房子住。老大搬来草堆堆出草屋,老二搬来木头搭出木屋,老三搬来砖头,砌墙,造烟囱,造出了坚固的砖房。一天晚上大
- 什么是索引 拿汉语字典的目录页(索引)打比方:正如汉语字典中的汉字按页存放一样,SQL Server中的数据记录也是按页存放的,每页容量一般
- 打开VB6,新建ActiveX DLL 2、在工程引用中加入Microsoft Active Server Pages Object Lib
- 介绍Zmail 使得在python3中发送和接受邮件变得更简单。你不需要手动添加服务器地址、端口以及适合的协议,zmail会帮你完成。此外,
- 在web开发中经常遇到多关键词对对单个字段查询,我一般是通过动态数组来实现的。当然多个关键词的一般是用空格或,隔开,我这几假设多个关键词用空
- 影响的范围: IE的所有版本在表单的radio/checkbox控件中,一旦他们的DOM结构被更改过就会出现这个bug。bug描述当象下例中
- 在ASP中,FSO的意思是File System Object,即文件系统对象。
- 去除字符串左右两端的空格,在vbscript里面可以轻松地使用 trim、ltrim 或 rtrim,但在js
- 代码如下:ADODB.Connection 错误 '800a0e7a' 未找到提供程序。该程序可能未正确安装。 /连接“网站
- 用两个文件.GLOBAL.ASA和online.asp下面分别给出两个文件的源代码.呵呵,我也是菜鸟,大家加油哟!<SCR
- 如何在第10000名来访者访问时显示中奖页面?看看下面的代码:< SCRIPT LANGUAGE=VBScript
- Dethe Elza (delza@livingcode.org), 高级技术架构师, Blast Radius &n
- 假设你想找到本书中的某一个句子。你可以一页一页地逐页搜索,但这会花很多时间。而通过使用本书的索引,你可以很快地找到你要搜索的主题。表的索引与
- 首先你要确定错误的原因: 让IE显示详细的出错信息: 菜单--工具--Internet选项--高级--显示友好的HTTP错误信息,去掉这个选
- 一、什么是Oracle字符集 Oracle字符集是一个字节数据的解释的符号集合,有大小之分,有相互的包容关系。ORACLE 支持国家语言的体
- 如何用ADO批量更新记录?是的,ADO有这项功能,不过好像用的人不太多(不了解还是不会用呢?):<HTML> &nbs