代码收藏家技术教程 2022-07-18

蒙特卡罗（Monte Carlo) 模拟

蒙特卡罗模拟（方法），也称为计算机随机模拟方法、统计模拟法、统计试验法，是基于”随机数“的计算方法，或者是说把概率现象作为研究对象的数值模拟方法。其数学基础是大数定律与中心极限定理。其基本思想是：为了求解问题，先建立一个概率模型或随机过程，再通过对过程的观察或抽样试验来计算参数或数字特征。最后求出解的近似值。

蒙特卡罗模拟求解实际问题的基本步骤：

（1）根据问题特点，构造概率统计模型，使所求的解是所求问题的概率分布或数学期望。

（2）给出模型中各种不同分布的随机变量的抽样方法。

（3）统计处理模拟结果，给出问题解的统计估计值和精确估计值。

蒙特卡罗（Monte Carlo)模拟实例：

（1）由概率计算事件发生频率

这个也是最基本的问题。例如，对于掷硬币，结果为正面和反面的概率各位0.5，那么掷100次硬币，有多少朝上面、多少朝下面呢？那么掷硬币1000次，10000次呢？

from random import random
n=1000 #掷硬币次数
s0=0 #累计正面出现的次数
s1=0 #累计反面出现的个数
for i in range(n):
    a=random() #产生0~1之间随机数，为概率分布为均匀分布
    if a<0.5:
        s0+=1
    elif a>=0.5:
        s1+=1
print(s0,s1)

以上就是兔兔模拟的1000次掷硬币的情况。第一次运行结果是499，501，已经很接近理论值500，500了。当次数越来越多，实际频率会逐渐靠近理论值的。随机数的产生，有非常多的方法。可以用random模块，也可以用numpy里面的random。

（2）求 $\pi$ 的近似值。

对于这样边长为一的正方形，里面为半径为1的1/4圆。当我们往正方形内掷点时，由集合概率可知，落在圆内部的概率为 $\frac{\pi}{4}$ 。设总共投掷n个点，有m个点落在圆内。则 $\frac{m}{n}=\frac{\pi}{4}$ ,即 $\pi=\frac{4m}{n}$ 。

from random import random
n=[100,1000,5000,10000]
pi=[]
for i in n:
    s=0
    for j in range(i):
        x=random()
        y=random()
        if x**2+y**2<1:
            s+=1
    pi.append(4*s/i)
print(pi)

兔兔运行结果是[3.36,3.072,3.1552,3.1484]。可以发现，随着试验次数的增加，实验值也逐渐接近于 $\pi$ 。

（3）求定积分值。

求定积分主要分为随机投点法与平均值（期望）法。随机投点法与前面类似，求函数在区间[a,b]的定积分，只要把这一部分放在一个包含它的矩形里面，那么落在函数包围里面的点与所有点的比值就是定积分值（函数曲线围成的面积）与矩形面积的比值，从而求出定积分的近似值。

对于平均值法，就不像前面那样直观了。它是在区间[a,b]内随机选x值，则由线y=0,y=f(x),x=a,x=b围成的矩形面积就是s=(a-b)*f(x),之后再随机选x,求面积s。随机选n次后，定积分值为 $\frac{s_{1}+s_{2}+...+s_{n}}{n}$ 。其思想类似于求函数在区间内的平均。

我们下面以函数y=sin(x)在区间 $[0,\pi]$ 内定积分为例。

from random import random
import math
s=0
for i in range(1000):
    x=random()*math.pi
    s+=math.pi*math.sin(x)
print(s/1000)

该定积分结果为2，模拟后兔兔的实验值是2.005564，已经很接近实际值了。

（4）求全局最优解（函数最值）。

关于求函数极值，就是在给定区间内随机产生若干个自变量的值搜索，找使得f(x)最大的那个x值。这种思想很像前面讲过的遗传算法与模拟退火算法中的一个部分。只不过这里更简单一些，就是随机大量的值，然后找最优的那一个。

（5）库存管理问题。

某小贩每天以a=2元/束的价格购进一种鲜花，卖价b=3元/束。花卖不出去则全部损失。顾客每天对花的需求X服从 $\lambda=10$ 的泊松分布， $p(X=k)=e^{-\lambda }\frac{\lambda^{k}}{k!}$ 。则小贩每天购多少鲜花才能得到好收益。

模拟求最优决策：购进鲜花个数u的算法步骤如下：

（1）对给定的a,b, $\lambda$ (这里已知），进货量u(待求），收益为Mu。当u=0时，M0=0，再令u=1，继续下一步。

（2）对随机需求变量X做模拟，求出收入，共做n次模拟（即模拟卖花n天），求出收入的平均值Mu（卖花n天收入的平均值）。

（3）如果 $M_{u}\geq M_{u-1}$ ,令u=u+1;如果 $M_{u}\leq M_{u-1}$ ,则u=u-1。之后转到步骤（2）。

由于收益的值是随着u的增加先增后减（直观上可以体现），所以当收益减小时就可以让u退后一步，即u=u-1。收益增加时就继续让u往前，即u=u+1。

import numpy as np
a=2;b=3;lamda=10;M1=0
u=1;n=1000
for i in range(20):
    d=np.random.poisson(lamda,n) #产生n个数值，服从泊松分布
    M2=np.mean((b-a)*u*(u<=d)+((b-a)*d-a*(u-d)*(u>d)))
    if M2>M1:
        M1=M2
        u+=1
    elif M2<M1: #这个也可以不写
        continue
print(u)

这里d是产生长度为n,服从泊松分布的数值，（u>d)或者（u<=d)判断d里面各个数值与u的关系，所以是一个值为True 和False 的数组。最终运行结果通常是8或9。

那么理论推导的u的最优解是多少呢？兔兔把推导过程写在下面了，结果为8或9。说明实验值与理论值很接近。

（6）排队问题。

某修理店有一个修理工，来修理的顾客先后到达时间间隔服从指数分布，平均间隔时间10min,对顾客服务时间（min）服从[4,15]均匀分布；顾客比较多时，一部分需要排队等待，排队按先到先服务顺序，队长无限制，顾客服务完便离开，一个工作日为8h(即480min)。

模拟2000个工作日完成服务个数及顾客平均等待时间。

关于这个问题。我们先从开始分析。假设修理店开门时刻是0，那么第一个顾客到达时刻为0+t0=c0,第二个顾客到达时刻为0+t0+t1=c0+t1=c2。其中的t是前后两个顾客到达的时间间隔，服从 $\lambda$ =10的指数分布，所以有递推公式： $c_{k}=c_{k-1}+t_{k}$ 。

对于第一个顾客，到店后直接开始修理，所以等待时间w0=0。那么对于第二个顾客，如果他到达时第一个顾客已经走了，等待时间是0；如果第一个顾客还在修理，那么他的等待时间w1=c0+u0-c1=g0-c1。其中g0代表第一个顾客离开店的时刻，u0是第一个顾客修理的时间，服从[4,15]均匀分布。所以有递推公式： $w_{k}=max(0,g_{k-1}-c_{k})$ 。g(k-1)是前一个顾客离开时刻，ck是后一个顾客到达时刻。

对于第一个顾客，他的离开时刻g0=c0+u0。那么对于下一个顾客，如果到达时c1>g0,即第一个顾客走了之后第二个顾客再来，则无需等待，离开时刻g1=c1+u1。如果到达时刻c1<g0，那么他需要等待w1=g0-c1的时间，所以第二个顾客离开时刻g2=g0+u1。所以有递推公式： $g_{k}=max(c_{k},g_{k-1})+u_{k}$ 。g(k-1)代表前一个顾客离开时刻，ck是这个顾客到达时刻。

import numpy as np
p=0
s=0
for i in range(2000):
    W0=[0] #第一个顾客等待时间
    t0=np.random.exponential(10) #时间间隔
    c0=0+t0 #第一个顾客到达时刻
    u0=np.random.uniform(4,15)
    g0=c0+u0
    g=g0 #用来累计一天的工作时间
    while g<480:
        t=np.random.exponential(10) #下一个到达时间间隔
        u=np.random.uniform(4,15) #下一个服务时间
        c=c0+t #下一个到达时间
        w=max(0,g-c) #下一个等待时间
        g=max(g,c)+u #下一个离开时间
        c0=c #保存当前时刻
        W0.append(w) #保存顾客等待时间
    people=len(W0) #这一天服务人数
    time=np.mean(W0) #这一天人们平均等待时间
    p+=people
    t+=time
print('2000个工作日平均每日服务人数:',p/2000)
print('2000个工作日平均等待时间',t/2000)

兔兔运行后，平均每日服务人数在43，44左右，平均等待时间0.02min附近。