python计算定积分_python编程通过蒙特卡洛法计算定积分详解

这篇文章主要介绍了python编程通过蒙特卡洛法计算定积分详解,具有一定借鉴价值,需要的朋友可以参考下。

想当初,考研的时候要是知道有这么个好东西,计算定积分。。。开玩笑,那时候计算定积分根本没有这么简单的。但这确实给我打开了一种思路,用编程语言去解决更多更复杂的数学问题。下面进入正题。

0b83b92ea128fc84bfef313422314c2e-0.png

如上图所示,计算区间[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, high=b, 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

0b83b92ea128fc84bfef313422314c2e-1.png

从上图可以看出,随着采样点数的增加,计算误差逐渐减小。想要提高模拟结果的精确度有两个途径:其一是增加试验次数N;其二是降低方差σ2. 增加试验次数势必使解题所用计算机的总时间增加,要想以此来达到提高精度之目的显然是不合适的。下面来介绍重要抽样法来减小方差,提高积分计算的精度。

重要性抽样法的特点在于,它不是从给定的过程的概率分布抽样,而是从修改的概率分布抽样,使对模拟结果有重要作用的事件更多出现,从而提高抽样效率,减少花费在对模拟结果无关紧要的事件上的计算时间。比如在区间[a b]上求g(x)的积分,若采用均匀抽样,在函数值g(x)比较小的区间内产生的抽样点跟函数值较大处区间内产生的抽样点的数目接近,显然抽样效率不高,可以将抽样概率密度函数改为f(x),使f(x)与g(x)的形状相近,就可以保证对积分计算贡献较大的抽样值出现的机会大于贡献小的抽样值,即可以将积分运算改写为:

aef36d58cb08c36504324c157240d903-2.png

x是按照概率密度f(x)抽样获得的随机变量,显然在区间[a b]内应该有:

aef36d58cb08c36504324c157240d903-3.png

因此,可容易将积分值I看成是随机变量 Y = g(x)/f(x)的期望,式子中xi是服从概率密度f(x)的采样点

aef36d58cb08c36504324c157240d903-4.png

下面的例子采用一个正态分布函数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,2)

plt.hist(Iis,30, histtype='step', label=u'Importance Sampling');

plt.hist(Ivmc, 30, color='r',histtype='step', label=u'Vanilla MC');

plt.vlines(np.pi, 0, 100, color='g', linestyle='dashed')

plt.legend()

plt.show()

bb29e27d4fa60419f01c2b01cb50ca61-5.png

从图中可以看出曲线sin(x)*x的形状和正态分布曲线的形状相近,因此在曲线峰值处的采样点数目会比曲线上位置低的地方要多。精确计算的结果为pi,从上面的右图中可以看出:两种方法均计算定积分1000次,靠近精确值pi=3.1415处的结果最多,离精确值越远数目越少,显然这符合常规。但是采用传统方法(红色直方图)计算出的积分值方的差明显比采用重要抽样法(蓝色直方图)要大。因此,采用重要抽样法计算可以降低方差,提高精度。另外需要注意的是:关于函数f(x)的选择会对计算结果的精度产生影响,当我们选择的函数f(x)与g(x)相差较大时,计算结果的方差也会加大。

相关推荐:

以上就是python编程通过蒙特卡洛法计算定积分详解的详细内容,更多请关注php中文网其它相关文章!

article_wechat2021.jpg?1111

本文原创发布php中文网,转载请注明出处,感谢您的尊重!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/566971.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

java容易掉发吗_容易被忽略的面试题—Java高并发

常见实现如CAS等。部分乐观锁削弱了一致性,但中低并发程度下的效率大大提高。并发编程Java中如何创建一个线程?从面相接口的角度上讲,实际上只有一种方法实现Runable接口;但Thread类为线程操作提供了更多的支持,所以通…

java spring mvc api_SpringMVC实现REST API

JSON使用Jackson jar包、RequestBody、ResponseBody注解,达到:1. 请求JSON消息体映射为JAVA对象2. 返回JAVA对象映射为JSON消息体Step 1. 导入Jackson jar包:Step 2. 在Spring MVC配置中加入annotation-driven,该配置可以确保Requ…

java 定义构造器_java的构造器定义以及使用

构造器,是面向对象所特有的概念,是一种特殊的方法,与对象创建有关1、构造器没有返回值类型2、构造器方法名与类名相同,而且可以重载构造器3、构造器不能手动调用,只能在创建对象时自动调用一次4、如果没有在类中定义构…

gateway动态路由_无语!SpringCloud Gateway动态路由之Nacos,我已经讲得很清楚了

前言当我们的网关Gateway程序开发完成之后,需要部署到生产环境,这个时候你的程序不能是单点运行的,肯定是多节点启动(独立部署或者docker等容器部署),防止单节点故障导致整个服务不能访问,网关是对客户端的入口与出口&…

json里面的list数据取不出来_sql盲注的困局:利用DNSlog快速导出数据

对于一个sql注入点来说最幸运的就是支持堆叠注入&#xff0c;最蛋疼的就是盲注&#xff0c;盲注里面难搞的就是基于时间的盲注。我们在本地利用这段代码进行演示<?php error_reporting(0); $link mysqli_connect(localhost,root,root); mysqli_set_charset($link,utf8); m…

python判断文件是否打开_python如何判断一个文件是否处于打开状态?

请问你的是linux吗&#xff1f;如果是linux可以借助/proc来获取。 import os class File(object): def __init__(self, file_path): if not os.path.exists(file_path): raise OSError({file_path} not exist.format(file_path file_path)) self.file_path os.path.abspath(f…

python怎么测试程序_python如何测试程序

测试函数是用于自动化测试&#xff0c;使用python模块中的unittest中的工具来测试 附上书中摘抄来的代码&#xff1a;#codingutf-8 import unittest from name_function import get_formatted_name class NamesTestCase(unittest.TestCase): def test_first_last_name(self): f…

判定覆盖白盒测试java_白盒测试系列(四)条件判定覆盖

条件判定覆盖一、定义&#xff1a;程序中每个判定至少有一次为真值&#xff0c;有一次为假值,使得程序中每个分支至少执行一次&#xff0c;且使得各判定中的每个条件获得各种可能的取值至少满足一次。二、特点&#xff1a;1、综合了条件覆盖和判定覆盖的特点2、满足条件判定覆盖…

discard python_Netty入门教程(一) 实现DISCARD服务

官方那个给出的介绍是&#xff1a;Netty是由JBOSS提供的一个java开源框架。Netty提供异步的、事件驱动的网络应用程序框架和工具&#xff0c;用以快速开发高性能、高可靠性的网络服务器和客户端程序。然后我们简单理解一下&#xff0c;这玩意就是个程序&#xff0c;干什么的&am…

python 每天执行一次_python 定时器每天就执行一次的实现代码

1.实现功能编写python脚本一直运行&#xff0c;判断当下是否是新的一天&#xff0c;如果是就执行一次任务代码2.具体实现代码#-*-coding:utf-8 -*-__author__ Administratorimport os,threading,timecurTimetime.strftime("%Y-%M-%D",time.localtime())#记录当前时间…

python 类的内置方法_【转】[python] 类常用的内置方法

原文&#xff1a;http://xukaizijian.blog.163.com/blog/static/170433119201111894228877/ 内置方法 说明 __init__(self,...) 初始化对象&#xff0c;在创建新对象时调用 __del__(self) 释放对象&#xff0c;在对象被删除之前调用 __new__(cls,*args,**kwd) 实例的生成操作 _…

java 文件与base64_java之文件与base64字符之间的相互转换

package cn.xuanyuan.util;import java.io.File;import java.io.FileInputStream;import java.io.FileOutputStream;import sun.misc.BASE64Decoder;import sun.misc.BASE64Encoder;public class FileUitl {/*** 将文件转成base64 字符串* param path文件路径* return ** thro…

python向量化编程技巧_神经网络基础之Python与向量化

Vectorization 深度学习算法中&#xff0c;数据量很大&#xff0c;在程序中尽量减少使用loop循环语句&#xff0c;而可以使用向量运算来提高程序运行速度。 向量化(Vectorization)就是利用矩阵运算的思想&#xff0c;大大提高运算速度。例如下面所示在Python中使用向量化要比使…

java 父子线程 调用链_ZipKin原理学习--Zipkin多线程及线程池中追踪一致性问题解决...

在学习Zipkin分布式追踪系统中我们了解到Trace在整个调用链是一致的&#xff0c;在web服务中可以通过在header设置Trace值在不同的服务中进行传递&#xff0c;那样在一个服务内部不同的线程&#xff0c;甚至是线程池中Zipkin是如何处理的&#xff0c;接下来我们来了解学习一下。…

SQL server 数据库面试题及答案(实操2)

使用你的名字创建一个数据库 创建表&#xff1a; 数据库中有三张表&#xff0c;分别为student,course,SC&#xff08;即学生表&#xff0c;课程表&#xff0c;选课表&#xff09; 问题&#xff1a; --1.分别查询学生表和学生修课表中的全部数据。--2.查询成绩在70到80分之间…

python电子相册制作软件_电子相册怎么做

电子相册制作 本文来自#千兆网络有什么用#征稿活动&#xff0c;不断提速的网络给你的生活带来了什么变化&#xff1f;快来参与活动&#xff0c;聊聊你玩转互联网&#xff0c;高速网上冲浪的经历&#xff01;>点击这里查看活动详情< 现在手机的拍照功能日趋强大&#xff0…

java list 范围删除_JAVA中循环删除list中元素(移除list两时间范围外的元素)

印象中循环删除list中的元素使用for循环的方式是有问题的&#xff0c;但是可以使用增强的for循环&#xff0c;然后今天在使用时发现报错了&#xff0c;然后去科普了一下&#xff0c;再然后发现这是一个误区。下面就来讲一讲。。伸手党可直接跳至文末。看总结。。JAVA中循环遍历…

python标注审核_Python类型标注

机器学习越来越火&#xff0c;大量的机器学习包都支持Python&#xff0c;导致了Python近几年非常火爆&#xff0c;入手门槛低&#xff0c;编程简单&#xff0c;概念非常少。越来越多的新手小白加入到Python编程。 Python虽然简单&#xff0c;但也带来很多问题。尤其是弱类型一直…

php的在线问卷调查_基于php技术的问卷调查系统

本系统前台主要使用php作为开发语言&#xff0c;后台使用mysql作为数据库管理系统&#xff0c;开发环境是wamp&#xff0c;服务器采用apache。系统的主要功能包括&#xff1a;管理登陆、问卷调查题目及内容选项的添加、修改和查询&#xff0c;调查结果统计等。分为管理员用户、…

python reduce函数_Python reduce()函数的用法小结

reduce()函数也是Python内置的一个高阶函数。 reduce() 格式&#xff1a; reduce (func, seq[, init()]) reduce()函数即为化简函数&#xff0c;它的执行过程为&#xff1a;每一次迭代&#xff0c;都将上一次的迭代结果&#xff08;注&#xff1a;第一次为init元素&#xff0c;…