python算积分蒙特卡罗_python编程通过蒙特卡洛法计算定积分详解
想当初,考研的时候要是知道有这么个好东西,计算定积分。。。开玩笑,那时候计算定积分根本没有这么简单的。但这确实给我打开了一种思路,用编程语言去解决更多更复杂的数学问题。下面进入正题。
如上图所示,计算区间[a b]上f(x)的积分即求曲线与X轴围成红色区域的面积。下面使用蒙特卡洛法计算区间[2 3]上的定积分:∫(x2+4*x*sin(x))dx
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return x**2 + 4*x*np.sin(x)
def intf(x):
return x**3/3.0+4.0*np.sin(x) - 4.0*x*np.cos(x)
a = 2;
b = 3;
# use N draws
N= 10000
X = np.random.uniform(low=a,high=b,size=N) # N values uniformly drawn from a to b
Y =f(X) # CALCULATE THE f(x)
# 蒙特卡洛法计算定积分:面积=宽度*平均高度
Imc= (b-a) * np.sum(Y)/ N;
exactval=intf(b)-intf(a)
print "Monte Carlo estimation=",Imc,"Exact number=",intf(b)-intf(a)
# --How does the accuracy depends on the number of points(samples)? Lets try the same 1-D integral
# The Monte Carlo methods yield approximate answers whose accuracy depends on the number of draws.
Imc=np.zeros(1000)
Na = np.linspace(0,1000,1000)
exactval= intf(b)-intf(a)
for N in np.arange(0,1000):
X = np.random.uniform(low=a,size=N) # N values uniformly drawn from a to b
Y =f(X) # CALCULATE THE f(x)
Imc[N]= (b-a) * np.sum(Y)/ N;
plt.plot(Na[10:],np.sqrt((Imc[10:]-exactval)**2),alpha=0.7)
plt.plot(Na[10:],1/np.sqrt(Na[10:]),'r')
plt.xlabel("N")
plt.ylabel("sqrt((Imc-ExactValue)$^2$)")
plt.show()
>>>
Monte Carlo estimation= 11.8181144118 Exact number= 11.8113589251
从上图可以看出,随着采样点数的增加,计算误差逐渐减小。想要提高模拟结果的精确度有两个途径:其一是增加试验次数N;其二是降低方差σ2. 增加试验次数势必使解题所用计算机的总时间增加,要想以此来达到提高精度之目的显然是不合适的。下面来介绍重要抽样法来减小方差,提高积分计算的精度。
重要性抽样法的特点在于,它不是从给定的过程的概率分布抽样,而是从修改的概率分布抽样,使对模拟结果有重要作用的事件更多出现,从而提高抽样效率,减少花费在对模拟结果无关紧要的事件上的计算时间。比如在区间[a b]上求g(x)的积分,若采用均匀抽样,在函数值g(x)比较小的区间内产生的抽样点跟函数值较大处区间内产生的抽样点的数目接近,显然抽样效率不高,可以将抽样概率密度函数改为f(x),使f(x)与g(x)的形状相近,就可以保证对积分计算贡献较大的抽样值出现的机会大于贡献小的抽样值,即可以将积分运算改写为:
x是按照概率密度f(x)抽样获得的随机变量,显然在区间[a b]内应该有:
因此,可容易将积分值I看成是随机变量 Y = g(x)/f(x)的期望,式子中xi是服从概率密度f(x)的采样点
下面的例子采用一个正态分布函数f(x)来近似g(x)=sin(x)*x,并依据正态分布选取采样值计算区间[0 pi]上的积分个∫g(x)dx
# -*- coding: utf-8 -*-
# Example: Calculate ∫sin(x)xdx
# The function has a shape that is similar to Gaussian and therefore
# we choose here a Gaussian as importance sampling distribution.
from scipy import stats
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
mu = 2;
sig =.7;
f = lambda x: np.sin(x)*x
infun = lambda x: np.sin(x)-x*np.cos(x)
p = lambda x: (1/np.sqrt(2*np.pi*sig**2))*np.exp(-(x-mu)**2/(2.0*sig**2))
normfun = lambda x: norm.cdf(x-mu,scale=sig)
plt.figure(figsize=(18,8)) # set the figure size
# range of integration
xmax =np.pi
xmin =0
# Number of draws
N =1000
# Just want to plot the function
x=np.linspace(xmin,xmax,1000)
plt.subplot(1,2,1)
plt.plot(x,f(x),'b',label=u'Original $x\sin(x)$')
plt.plot(x,p(x),'r',label=u'Importance Sampling Function: Normal')
plt.xlabel('x')
plt.legend()
# =============================================
# EXACT SOLUTION
# =============================================
Iexact = infun(xmax)-infun(xmin)
print Iexact
# ============================================
# VANILLA MONTE CARLO
# ============================================
Ivmc = np.zeros(1000)
for k in np.arange(0,1000):
x = np.random.uniform(low=xmin,high=xmax,size=N)
Ivmc[k] = (xmax-xmin)*np.mean(f(x))
# ============================================
# IMPORTANCE SAMPLING
# ============================================
# CHOOSE Gaussian so it similar to the original functions
# Importance sampling: choose the random points so that
# more points are chosen around the peak,less where the integrand is small.
Iis = np.zeros(1000)
for k in np.arange(0,1000):
# DRAW FROM THE GAUSSIAN: xis~N(mu,sig^2)
xis = mu + sig*np.random.randn(N,1);
xis = xis[ (xisxmin)] ;
# normalization for gaussian from 0..pi
normal = normfun(np.pi)-normfun(0) # 注意:概率密度函数在采样区间[0 pi]上的积分需要等于1
Iis[k] =np.mean(f(xis)/p(xis))*normal # 因此,此处需要乘一个系数即p(x)在[0 pi]上的积分
plt.subplot(1,2)
plt.hist(Iis,30,histtype='step',label=u'Importance Sampling');
plt.hist(Ivmc,color='r',label=u'Vanilla MC');
plt.vlines(np.pi,100,color='g',linestyle='dashed')
plt.legend()
plt.show()
从图中可以看出曲线sin(x)*x的形状和正态分布曲线的形状相近,因此在曲线峰值处的采样点数目会比曲线上位置低的地方要多。精确计算的结果为pi,从上面的右图中可以看出:两种方法均计算定积分1000次,靠近精确值pi=3.1415处的结果最多,离精确值越远数目越少,显然这符合常规。但是采用传统方法(红色直方图)计算出的积分值方的差明显比采用重要抽样法(蓝色直方图)要大。因此,采用重要抽样法计算可以降低方差,提高精度。另外需要注意的是:关于函数f(x)的选择会对计算结果的精度产生影响,当我们选择的函数f(x)与g(x)相差较大时,计算结果的方差也会加大。
总结
以上就是本文关于python编程通过蒙特卡洛法计算定积分详解的全部内容,希望对大家有所帮助。感兴趣的朋友可以继续参阅本站:
如有不足之处,欢迎留言指出。感谢朋友们对本站的支持!
python算积分蒙特卡罗_python编程通过蒙特卡洛法计算定积分详解相关推荐
- python计算定积分_python编程通过蒙特卡洛法计算定积分详解
这篇文章主要介绍了python编程通过蒙特卡洛法计算定积分详解,具有一定借鉴价值,需要的朋友可以参考下. 想当初,考研的时候要是知道有这么个好东西,计算定积分...开玩笑,那时候计算定积分根本没有这么 ...
- python用蒙特卡洛法区间_python编程通过蒙特卡洛法计算定积分详解
想当初,考研的时候要是知道有这么个好东西,计算定积分...开玩笑,那时候计算定积分根本没有这么简单的.但这确实给我打开了一种思路,用编程语言去解决更多更复杂的数学问题.下面进入正题. 如上图所示,计算 ...
- python多态的例子_Python编程之多态用法实例详解
本文实例讲述了Python编程之多态用法.分享给大家供大家参考.具体分析如下: 什么是多态?顾名思义,多态就是多种表现形态的意思.它是一种机制.一种能力,而非某个关键字.它在类的继承中得以实现,在类的 ...
- python布尔类型运算_Python对象类型及其运算方法(详解)
基本要点: 程序中储存的所有数据都是对象(可变对象:值可以修改 不可变对象:值不可修改) 每个对象都有一个身份.一个类型.一个值 例: >>> a1 = 'abc' >> ...
- python 自动化办公 案例_python自动化工具之pywinauto实例详解
python自动化工具之pywinauto实例详解 来源:中文源码网 浏览: 次 日期:2019年11月5日 [下载文档: python自动化工具之pywinauto实例详解.txt ] (友情提示: ...
- python属性使用教程_Python对象的属性访问过程详解
只想回答一个问题: 当编译器要读取obj.field时, 发生了什么? 看似简单的属性访问, 其过程还蛮曲折的. 总共有以下几个step: 1. 如果obj 本身(一个instance )有这个属性, ...
- python dataframe loc函数_python pandas.DataFrame.loc函数使用详解
官方函数 DataFrame.loc Access a group of rows and columns by label(s) or a boolean array. .loc[] is prim ...
- python迭代器创建序列_Python 中迭代器与生成器实例详解
Python 中迭代器与生成器实例详解 本文通过针对不同应用场景及其解决方案的方式,总结了Python中迭代器与生成器的一些相关知识,具体如下: 1.手动遍历迭代器 应用场景:想遍历一个可迭代对象中的 ...
- python装饰器性能_python装饰器的特性原理详解
这篇文章主要介绍了python装饰器的特性原理详解,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 今天发现了装饰器的另一种用法,下面就先上代码: d ...
- python迭代器与生成器_python的迭代器与生成器实例详解
本文以实例详解了python的迭代器与生成器,具体如下所示: 1. 迭代器概述: 迭代器是访问集合元素的一种方式.迭代器对象从集合的第一个元素开始访问,直到所有的元素被访问完结束.迭代器只能往前不会后 ...
最新文章
- React navtive
- 【唠叨两句】如何将一张树型结构的Excel表格中的数据导入到多张数据库表中...
- 33.Linux系统介绍
- 解决Mac Pro上IDEA卡顿的问题
- python话雷达图-Python简单雷达图绘制
- mysql数据库的增删改查命令_MySQL 初识别语句,数据库、表、行的增删改查
- MFC子线程访问主线程对话框程序的控件对象
- java 反射泛型方法_java基础之反射和泛型以及注解
- cfb为什么不需要填充_学日语为什么不需要准备,现在就可以学?
- cin gt gt a用c语言怎么写写,cin、cin.get()、cin.getline()、getline()、gets()等函数的用法...
- Flex接受任意拖拽
- oracle中序号生成器,Oracle序列生成器
- php header 无法跳转,PHP利用header跳转失效解决方法
- GCN实践——可视化cora-network
- window10运行不了1stopt_1stopt点击运行没有反应求大佬指点
- 《缠论》的精髓是什么?
- 计算机机械硬盘系统安装,电脑硬盘安装图解,机械硬盘安装-
- 软件生命周期-SDLC-的六个阶段简单介绍
- Mega2560串口通信实现
- MapStruct系列(6)-映射集合、映射Stream流、映射枚举
热门文章
- p图软件pⅰc_P图教程|教你做超火的iMessage图 所需软件:Picsart QQ_修图软件_滤镜_picsart怎么样_纯白色_相册_我超会p图der_摄影_摄影技巧_修图技巧...
- Corrupted STDOUT by directly writing to native stream in forked JVM 1. See FAQ web page and the dump
- Android反编译工具合集
- (三)Detecting Spacecraft Anomalies Using LSTMs and Nonparametric Dynamic Thresholding
- 网页设计的步骤和标准都有哪些?
- NUC搭建Centos8服务器
- 超搜索引擎BBMAO
- 利用Matplotlib绘制各类图表
- 天津大学仁爱学院c语言期末考试题,天津大学仁爱学院2014-2015学年第1学期期末C语言复习...
- Sql Server——Sql Server中进行查询操作时提示“对象名无效”