python实现蒙特卡洛模拟_随机模拟——蒙特卡洛算法及Python实现_weixin_39808143的博客-程序员秘密

随机模拟——蒙特卡洛算法及Python实现

irony • 2020 年 04 月 24 日

蒙特卡洛算法介绍蒙特卡洛算法(Monte Carlo method)也称统计模拟算法，是一种以概率统计理论为指导的数值计算方法方法。是指使用随机数来解决很多计算问题的方法。

蒙特卡洛方法基本思想

通常蒙特卡洛算法可以分为两种：

1，一种是所求解的问题本身具有内在随机性，借助计算机的运算能力可以直接模拟这种随机的过程。

2，另一种是所求解的问题可以转化为某种随机分布的特征数，比如随机事件出现的概率，或者随机变量的期望值。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的多维积分的问题。

蒙特卡洛方法的理论依据是大数定理，即重复试验多次，随机事件的频率近似于它的概率。

论文中利用蒙特卡洛法求定积分

蒙特卡洛一个重要的应用就是求定积分，例子如下：

假设曲线函数为$f(x)$,则易得阴影部分面积:

$$ \int_{a}^{b} f(x)\, dx $$

然后通过牛顿莱布尼兹公式进行积分求解，而蒙特卡洛算法则运用不同的思想。当我们在$[a,b]$之间随机取一点$x$时,它对应的函数值就是$f(x)$。接下来我们就可以用$f(x) * (b-a)$来估计曲线下方的面积，也就是积分值，但是很显然这种估计方式是非常粗略的。

在上图中，做了四次随机采样，得到了四个随机样本$x_1,x_2,x_3,x_4$，并且得到了这四个样本的$f(x_i)$的值分别为$f(x_1),f(x_2),f(x_3),f(x_4)$。对于这四个样本，每个样本都可求取一个近似的面积值，大小为$f(x_i)*(b-a)$。然后求取所有面积的平均值就完成了蒙特卡洛求积分的过程。

用以下公式表述上述过程：

\begin{array}{lcr}

S = \frac{1}{4}(b-a)(f(x_1)+f(x_2)+f(x_3)+f(x_4))\\

\quad = \frac{1}{4}(b-a)\sum_{i=1}^4 f(x_i)

\end{array}

蒙特卡洛Python实现

1,求取圆周率$\pi$

代码如下：import random

def cal_pai():

n = 1000000

r = 1

a, b = (0, 0)

# 求投点的上下界

x_low, x_up = a - r, a + r

y_low, y_up = b - r, b + r

# 记录落于圆内的点数

count = 0

for i in range(0, n):

x = random.uniform(x_low, x_up)

y = random.uniform(y_low, y_up)

if x**2 + y**2 <= r**2:

count += 1

return (count / n) * 4

if __name__ == '__main__':

print(cal_pai())

简单介绍一下算法思路：利用投点法，在边长为2的正方形内部找到内接圆，内接圆的半径为1，所以易得内接圆于外接正方形的面积比为$\frac{\pi}{4}$。现在向正方形内部投点，并计算点与圆心的距离。其中$n$为投点的总数，$count$为投在圆内的投点个数，则$\pi$可以由以下格式求出：

$$ \pi = \frac{count}{n} * 4 $$

代码运行结果:$3.140748$，每次运行结果不同。依据大数定理，n越大,$\pi$值越准确。

2,求定积分$\int_{0}^{1} x^2\, dx$

代码如下：import random

def integral():

n = 1000000

x_min, x_max = 0, 1

S = 0

for i in range(0, n):

x = random.uniform(x_min, x_max)

y = x ** 2

S += (x_max - x_min) * y

return 1 / n * S

if __name__ == '__main__':

print(integral())

代码运行结果:$0.33371372617627443$

本文链接：https://blog.csdn.net/weixin_39808143/article/details/111404153

原作者删帖不实内容删帖广告或垃圾文章投诉

智能推荐

任意输入三个数，判断能否组成直角三角形_php任意给定三个正数,能否构成直角三角形,_海芝叭的博客-程序员秘密

#include <stdio.h>int main(){int a,b,c;printf(“输人第一条: “);scanf(”%d”,&a);printf(“老铁输入第二条: “);scanf(”%d”,&b);printf(“小铁输入第三条: “);scanf(”%d”,&c);int temp,result;temp=a+b>c&&a+c>b&&b+c>a;result=aa==bb+cc||

在ARM-Linux下实现车牌识别（一）------车牌提取_风筝丶的博客-程序员秘密

你好！这里是风筝的博客，欢迎和我一起交流。最近在弄车牌识别这个项目，对于机器视觉有些了解的人都知道，这个东西算是比较成熟了，在书里也有代码。网上能找到的资料也比较多，所及借着这个机会在ARM开发板上实现以下车牌识别。反正对于神经网络这些什么的我是不知道了，所以代码也是网上借鉴了的，我稍微整理注释了下。先放下移植opencv的步骤：移植opencv到嵌入式arm详细过程 ...

git学习笔记（三）—— 远程仓库_weixin_30832405的博客-程序员秘密

一、gitHub&&SSH Key　　为git仓库提供托管服务的，所以注册一个github账号，就可以免费获得git远程仓库。本地Git仓库和GitHub仓库之间的传输是通过SSH加密的，所以进行下面设置1、创建SSH Key　　在用户主目录下，看看有没有.ssh目录，如果有，再看看这个目录下有没有id_rsa和id_rsa.pub这两个文件，如果已经有...

在IDEA中添加h2数据库_idea连接h2数据库_爱吃面包的小胖墩的博客-程序员秘密

H2数据库简介H2 是一个用 Java 开发的嵌入式数据库，它本身只是一个类库，即只有一个 jar 文件，可以直接嵌入到应用项目中，占用空间小：大约2 MB，常被用作单元测试。步骤在IDEA中添加h2数据库的方法有两种，第一种：通过下载jar包，然后加载到IDEA中；第二种：通过添加依赖加载到IDEA中，在此仅介绍第二种方法。1.找h2数据库的依赖，一般都是通过Maven仓库来进行搜索，Maven仓库的网址： link；2.将依赖添加在IDEA项目的pom.xml中，并更新Maven仓库。

python自动修图_3 行 Python 代码 5 秒抠图的 AI 神器，根本无需 PS_weixin_39757212的博客-程序员秘密

原标题：3 行 Python 代码 5 秒抠图的 AI 神器，根本无需 PS 文 | 苏克 [email protected]高级农名工曾几何时，「抠图」是一个难度系数想当高的活儿，但今天要介绍的这款神工具，只要 3 行代码 5 秒钟就可以完成高精度抠图，甚至都不用会代码，点两下鼠标就完成了。感受下这款抠图工具抠地有多精细：是不是很赞？什么 PS 、PPT 修图都 low 爆了，你软件还没打开，我这边都修完了。...

kd树简介在matlab下VLFeat中的kd-tree使用_matlab中kd是什么_damant的博客-程序员秘密

先简要介绍knn——K近邻算法和kd-tree——kd树，然后介绍matlab环境中有关使用kd树的函数。k-d树（k-dimensional树的简称），是一种分割k维数据空间的数据结构。主要应用于多维空间关键数据的搜索（如：范围搜索和最近邻搜索）。这里用的是在VLFeat 实现的kd树，其是一个开源的计算机视觉库，实现了 SIFT,MSER, k-means, hierarchical k-means, aggl

随便推点

如果结束进程拒绝访问，可以尝试以下_margot无法结束进程_咚踏取的博客-程序员秘密

如果taskkill /f /pid 123 出现拒绝访问时，可使用以下方式删除进程：wmic process where name=‘qq.exe’ delete11如果这样还杀不死，恐怕就要进[安全模式]删除了。杀了进程，想删除文件，可以这样删除cmd下运行：DEL /F /A /Q 文件名抄自-吉吉教主...

计算机网络第三章课后答案第七版（谢希仁著）_计算机网络第三章课后题答案_敢敢J的憨憨L的博客-程序员秘密

计算机网络第三章课后答案简介：平时老师布置的作业自己做的答案，仅是做参考，希望有帮助还可以把全部的答案分享给大家：&gt;完整答案&lt;第三章课后作业答案：第一页第二页第三页第四页...

StringBuilder和StringBuffer的区别_dnbug Blog的博客-程序员秘密

StringBuilder和StringBuffer区别这里从两个角度(性能、线程安全)去分析区别。StringBuilder与StringBuffer性能比较走DEMO0StringBuilder与StringBuffer线程安全分析走DEMO1

UVA 11987 Almost Union-Find （并查集，可持久化）_guognib的博客-程序员秘密

大白书：p246 UVA 11987Almost Union-Find加一个id[]，不修改，而是insert，然后改idconst int MAXN = 100010 * 2;const int MOD = 1000000; int fa[MAXN], num[MAXN];int sum[MAXN]; int idx[MAXN]; void init(int n){

Command /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/clan_LANGZI7758521的博客-程序员秘密

ProcessPCH++ "/Users/xxxx/Library/Developer/Xcode/DerivedData/English_Reader-alybxowggrudboejoylghzohavko/Build/Intermediates/PrecompiledHeaders/English Reader-Prefix-fldqdedtnerfrubpymgklsjvzccz/Engl

a随便写的nodejs学习笔记(4)_在知识的海洋里奋力狗刨的博客-程序员秘密

nodejs解析请求参数01-nodejs接收get请求参数前端发送get请求的参数是直接在url后面拼接，格式是： url?key1=value1&key2=value2…//1.导入模块const http = require('http');const url = require('url');//解析get请求参数//2.创建服务器let server = http.createServer((req, res) => { //(1)将url中的中文进行转码