代码收藏家技术教程 2022-08-10

蒙特卡洛方法（入门详解）

一、定义

蒙特卡洛又称统计试验法，是基于概率论的算法。其实质就是将问题转化为一个概率问题，并用计算机模拟产生一堆随机数，再对随机数进行统计工作。

蒙特卡洛模拟方法=建立概率模型+计算机模拟+数理统计。

二、原理

大数定理证明,在大样本的情况下，事件发生频率等于概率

三、起源(布丰投针试验）

法国数学家设计了投针试验

1.取一张白纸，在上面画出许多条间距为a的平行线。

2.取一根长度为l（l<=a)的针，随机地向画有平行线的纸上投掷n次，观察针与直线相交的次数记为m.

3.计算针与直线相交的概率

对于布丰试验，我的解读如下：

针长l,间距为a,角度为[0,π】

1.理论概率：落下的针有可能和直线相交，有可能与直线不相交。针是否与直线相交取决于针中心到较近直线的距离和针的偏斜角度。当 $x\leq (l/2)*sin(\Theta )$ 时，相交。

由上图可以观察到只需要考虑针在偏斜的角度问题，每对应一个角度，针中心落在相应白色区域则不会相交，而针中心点在属于均匀分布，可以将针是否相交问题转化为 $d=(l/2)*sin(\Theta )$ 该长度在a的几何分布问题。因此每一个角度 $\Theta$ 对应一个是否相交的概率值 $f$ .men

本文我假设l=a=2,计算出理论概率为曲线lsinx/a下的面积与(o,π）*（0,1）矩形的比，即为2/π

import numpy as np
import math
import matplotlib.pyplot as plt
n=100 
l=2
a=2
phi=np.linspace(0,1,n)*math.pi#产生n个从小到大的角度值，角度范围在[0,180]
B=list(range(0,100))#假设投100次
fx=[]#概率值集合
for i in B:
    f=((l/2)*math.sin(phi[i]))/(a/2)
    fx.append(f)
plt.plot(phi,fx)
plt.show()#经过计算，针相交的概率为2L/πa(2/π）

2.模拟频率：模拟随机数，产生相交的频数及频率 $p=m/k$

由理论概率=模拟频率 $\frac{2*l}{\pi*a}=\frac{m}{k}=p$ ,可以推出 $\pi=\frac{2*l}{a*p}$

k=10000#模拟投10000根针
m=0
x=np.random.rand(k)*a/2#随机取10000根针中心与较近直线的距离
y=np.random.rand(k)*math.pi#随机取10000根针与直线的偏斜角度
#判断相交的针的数量m，我输出结果为6398
C=list(range(0,10000))
for i in C:
    if x[i]<=l/2*math.sin(y[i]):
        m=m+1
p=m/k#模拟结果的频率
mypi=2*l/(a*p)#计算π的值，我输出结果为3.125976867771178

3.一次模拟具有偶然性，可以模拟多次，计算平均值

k=10000

num=1
mm=[]
while num<=1000:
    m=0
    num=num+1
    x=np.random.rand(k)*a/2
    y=np.random.rand(k)*math.pi
    C=list(range(0,10000))
    for i in C:
        if x[i]<=l/2*math.sin(y[i]):
             m=m+1
    mm.append(m)
mmean=np.mean(mm)
p=mmean/k
mypi=2*l/(a*p)#我输出结果为3.1412944834785193

小编有话说

模拟随机数越大，结果越准确，由于小编我计算机运算速度不好，就简单做了下哦，但结果接近，证明小编基本思路是正确的。

四、蒙特卡洛计算积分

蒙特卡洛计算积分是利用面积的几何概率原理，积分一般情况等于计算面积，积分等于正面积-负面积，可以为负。小编用一个简单的定积分来说明蒙特卡洛的应用。

用蒙特卡洛计算 $\int_{4}^{7} x^{3}\varrho ^{x}dx$ 的值

n=10000#模拟10000次随机数
x=np.random.uniform(4,7,n)#产生（4,7）上均匀分布随机数作为x坐标
y=np.random.uniform(0,(7**3)*math.exp(7),n)#产生0到7对应函数值上均匀分布的随机数作为y坐标
#判断落在曲线和x轴随机数的数量
k=[]
for i in list(range(0,n)):
    b=y[i]-(x[i]**3)*math.exp(x[i])
    k.append(b)
num=len(list(filter(lambda k:k<0,k)))
p=num/n#计算落在积分面积内的随机数频率
result=p*(7-4)*(7**3)*math.exp(7)#根据频率计算积分面积

来源：梦在码下