蒙特卡罗（Monte Carlo) 模拟

Posted 2023-04-09 生信小兔

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了蒙特卡罗（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）求的近似值。

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

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]。可以发现，随着试验次数的增加，实验值也逐渐接近于。

（3）求定积分值。

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

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

我们下面以函数y=sin(x)在区间内定积分为例。

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服从的泊松分布，。则小贩每天购多少鲜花才能得到好收益。

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

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

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

（3）如果,令u=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是前后两个顾客到达的时间间隔，服从=10的指数分布，所以有递推公式：。

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

对于第一个顾客，他的离开时刻g0=c0+u0。那么对于下一个顾客，如果到达时c1>g0,即第一个顾客走了之后第二个顾客再来，则无需等待，离开时刻g1=c1+u1。如果到达时刻c1<g0，那么他需要等待w1=g0-c1的时间，所以第二个顾客离开时刻g2=g0+u1。所以有递推公式：。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附近。

总结

蒙特卡罗模拟本质还是利用随机数进行进行模拟计算。但同时也需要注意计算机产生的随机数其实也并不是真正的随机，相关实例兔兔在这里就不作展开了，在比较复杂的模拟的情况下是需要注意到这一点的。

蒙特卡罗算法（Monte Carlo method）

蒙特卡罗方法概述

蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。

蒙特卡罗方法的基本思想

用事件发生的“频率”来决定事件的“概率”。高速电子计算机使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。

蒙特卡罗方法的基本原理

设有统计独立的随机变量Xi(i=1，2，3，…，k)，其对应的概率密度函数分别为fx1，fx2，…，fxk，功能函数式为Z=g(x1，x2，…，xk)。

首先根据各随机变量的相应分布，产生N组随机数x1，x2，…，xk值，计算功能函数值 Zi=g(x1，x2，…，xk)(i=1，2，…，N)，若其中有L组随机数对应的功能函数值Zi≤0，则当N→∞时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。

从蒙特卡罗方法的思路可看出，该方法回避了结构可靠度分析中的数学困难，不管状态函数是否非线性、随机变量是否非正态，只要模拟的次数足够多，就可得到一个比较精确的失效概率和可靠度指标。

蒙特卡罗方法分子模拟计算的步骤

蒙特卡罗方法实施步骤：

1、通过敏感性分析，确定随机变量；

2、构造随机变量的概率分布模型；

3、为各输入随机变量抽取随机数；

4、将抽得的随机数转化为各输入随机变量的抽样值；

5、将抽样值组成一组项目评价基础数据；

6、根据基础数据计算出评价指标值；

7、整理模拟结果所得评价指标的期望值、方差、标准差和它的概率分布及累计概率，绘制累计概率分布图，计算项目可行或不可行的概率。

蒙特卡罗方法应用

1.求π

#include <bits/stdc++.h>

#define MAX_ITERS 10000000

using namespace std;

double Rand(double L, double R)
{
    return L + (R - L) * rand() * 1.0 / RAND_MAX;
}

double GetPi()
{
    srand(time(NULL));
    int cnt = 0;
    for(int i = 0; i < MAX_ITERS; i++)
    {
        double x = Rand(-1, 1);
        double y = Rand(-1, 1);
        if(x * x + y * y <= 1)
            cnt++;
    }
    return cnt * 4.0 / MAX_ITERS;
}

int main()
{
    for(int i = 0; i < 10; i++)
        cout <<  fixed << setprecision(10)<<GetPi() << endl;
    return 0;
}

2.求e

#include <bits/stdc++.h>

#define MAX_ITERS 10000000

using namespace std;

struct Point
{
    double x, y;
};

double Rand(double L, double R)
{
    return L + (R - L) * rand() * 1.0 / RAND_MAX;
}

Point getPoint()
{
    Point t;
    t.x = Rand(1.0, 2.0);
    t.y = Rand(0.0, 1.0);
    return t;
}

double getResult()
{
    int m = 0;
    int n = MAX_ITERS;
    srand(time(NULL));
    for(int i = 0; i < n; i++)
    {
        Point t = getPoint();
        double res = t.x * t.y;
        if(res <= 1.0)
            m++;
    }
    return pow(2.0, 1.0 * n / m);
}

int main()
{
    for(int i = 0; i < 20; i++)
        cout << fixed << setprecision(10) << getResult() << endl;
    return 0;
}
//precision() 返回当前的浮点数精度值
//precision(val) 设置val为新的浮点数精度值, 并返回原值
//setf(flags) 添加格式标志flags, 返回所有标志的原本状态.
//showpos    正数前面加上+号
//fixed 使用小数计数法
//scientific 使用科学计数法
//uppercase 使用大写字符
//showbase  显示数字的进制
//boolalpha    bool值使用字符表示 , true或者false
//noboolalpha     bool使用0和1表示
//left   靠左对齐
//right  靠右对齐
//internal 字符靠左对齐, 数字卡右对齐