python实现蒙特卡罗模拟法的实践
作者:洋洋菜鸟 时间:2023-08-11 02:22:02
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
标签:python,蒙特卡罗模拟法
![](/images/zang.png)
![](/images/jiucuo.png)
猜你喜欢
git stash的正确用法详解
2022-11-06 19:10:07
Mysql将一个表中的某一列数据复制到另一个表中某一列里的方法
2024-01-14 17:36:29
![](https://img.aspxhome.com/file/2023/0/68050_0s.png)
Python迭代器定义与简单用法分析
2022-10-24 02:44:29
ASP JSON类源码
2011-04-30 16:38:00
Python logging模块handlers用法详解
2022-08-24 18:08:08
学习ASP和编程的28个观点
2008-06-27 12:57:00
Python操作word文档的示例详解
2021-11-10 02:54:51
![](https://img.aspxhome.com/file/2023/1/93461_0s.png)
pymysql 插入数据 转义处理方式
2024-01-23 08:43:29
将滚动条(scrollbar)保持在最底部的方法
2008-02-21 10:05:00
Python利用redis-py实现哈希数据类型的常用指令操作
2021-09-24 18:16:45
![](https://img.aspxhome.com/file/2023/7/134777_0s.jpg)
Python流程控制语句的深入讲解
2023-12-01 00:06:31
![](https://img.aspxhome.com/file/2023/6/101066_0s.png)
Python中的exec、eval使用实例
2022-07-05 21:01:41
CentOS下安装Jenkins的完整步骤
2022-12-13 16:43:34
![](https://img.aspxhome.com/file/2023/0/132460_0s.jpg)
php strstr查找字符串中是否包含某些字符的查找函数
2023-11-17 01:42:23
vue 代码压缩优化方式
2024-04-09 10:44:46
在 Django/Flask 开发服务器上使用 HTTPS
2023-02-21 22:15:33
深入理解python中的闭包和装饰器
2023-03-20 06:27:48
![](https://img.aspxhome.com/file/2023/6/79766_0s.jpg)
django模型中的字段和model名显示为中文小技巧分享
2021-04-01 15:03:54
用sql脚本创建sqlserver数据库触发器范例语句
2024-01-27 16:42:14
MySQL选错索引的原因以及解决方案
2024-01-19 21:34:34
![](https://img.aspxhome.com/file/2023/8/102388_0s.png)