蒙特卡罗分析法也叫随机模拟法,python用蒙特卡洛方法计算定积分的值
蒙特卡罗就是生成随机变量,带入模型计算的结果。在优化方面,本文主要介绍python蒙特卡罗模拟方法的实现,通过实例代码详细介绍,具有一定的参考价值。感兴趣的朋友可以参考一下。
00-1010 1.导言2。例题分析2.1近似圆周率的模拟2.2定积分的估计2.3求解整数规划
目录
蒙特卡罗,又称随机抽样或统计检验,是产生随机变量并带入模型计算的结果。在优化方面,只要仿真次数足够多,最终可以找到最优解或近优解。
一般来说,蒙特卡罗方法大致可以分为两类:一类是所要解决的问题本身具有随机性,这种随机过程可以借助计算机的计算能力直接模拟出来。比如在核物理的研究中,分析中子在反应堆中的传输过程。中子和原子核之间的相互作用受到量子力学定律的限制。人们只能知道它们相互作用的概率,却无法准确得到中子和原子核的位置,以及裂变产生的新中子的行进速度和方向。科学家根据其概率随机取样,得到裂变的位置、速度和方向,从而在模拟大量中子的行为后,通过统计得到中子的透射范围,可以作为反应堆设计的依据。
另一类是要解决的问题可以转化为一些随机分布的特征数,比如随机事件的概率或者随机变量的期望值。通过随机抽样,随机事件的概率由它们发生的频率来估计,或者随机变量的数量由样本的数字特征来估计。
1.简介
2.实例分析
画单位圆和外接圆。正方形ABCD的面积为:2*2=4,圆的面积为:s= * 1 * 1=。现在,模拟产生N个均匀分布在正方形ABCD中的点。如果这N个点中的M个在圆内,圆的面积与平方ABCD的面积之比可以近似为m/n。
该过程如下:
#通过模拟计算近似圆周率
随机导入
将numpy作为np导入
将matplotlib.pyplot作为plt导入
#解决图片标题中的乱码问题。
将matplotlib作为mpl导入
mpl . RC params[ font . sans-serif ]=[ sim hei ]#指定默认字体
mpl . RC params[ axes . unicode _ MINUS ]=false #解决保存的图像显示为带负号-的正方形的问题
#进入正题
r=random.random()
num=范围(0,200000,10)
mypi=np.ones((1,len(num)))
对于范围内的j(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))
#计算圆周率的近似值(Monte Carlo:使用统计值来逼近真实值)
mypi[0,j]=4 * res/n
plt.plot(range(1,len(mypi[0]) 1),mypi[0],。-)
Plt.grid (ls= 3360 ,c= b ,)#打开坐标网格
PLT。AXHLINE (Y=NP。PI,LS= 3360 ,C= Yellow) #添加一条水平线
# plt.axvline (x=4,ls=-,c= green) #添加一条垂直线
Plt.legend([模拟,实际],loc=右上,散点=1)
Plt.title(“近似圆周率”)
plt.show()
e>
返回:
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)
返回:
到此这篇关于python实现蒙特卡罗模拟法的实践的文章就介绍到这了,更多相关python 蒙特卡罗模拟法内容请搜索盛行IT软件开发工作室以前的文章或继续浏览下面的相关文章希望大家以后多多支持盛行IT软件开发工作室!
郑重声明:本文由网友发布,不代表盛行IT的观点,版权归原作者所有,仅为传播更多信息之目的,如有侵权请联系,我们将第一时间修改或删除,多谢。