Table of Contents

1. 蒙特卡洛方法大致原理

蒙特卡洛起源于求面积，但受限于准确求解面积的方法，提出了利用概率的方法来粗略计算面积。

原理就是通过大量的数据采样，利用大数定律，由事件发生的频率得出相对应的概率。根据不同问题，可以抽象出类似需要进行概率计算的模型，从而利用该方法求解。

蒙特卡洛方法又称为统计模拟方法，与其概念完全一致。

蒙特卡洛方法的基本思想

  • 当问题可以抽象为某个确定的数学问题时，应当首先建立一个恰当的概率模型，也即确定某个随机事件A或者是随机变量X，使得待求的解等于随机事件出现的概率或随机变量的数学期望。

  • 然后进行模拟实验，即重复多次地模拟随机事件A或随机变量X。

  • 最后对随机实验结果进行统计平均，求出A出现地频数或X的平均值最为问题的近似解。这种方法也叫做间接蒙特卡洛模拟。

2. 蒙特卡洛方法实践步骤

蒙特卡洛方法的三个主要步骤：

  1. 描述或构造概率过程：对于实际存在概率随机的问题可以直接模拟这个概率过程，而对于不具有随机性质的确定性问题，需要人为来进行构造。

  2. 利用概率分布抽样：通过计算机产生已知概率分布的随机变量，常用的概率分布有均匀分布，正态分布、指数分布、泊松分布等。

  3. 计算所求解的估计量：通过构造的概率模型与抽样结果，进行估计量的测定。一般是把 次随机抽样结果的算术平均值作为解的近似值。

其中产生已知概率分布的随机变量是蒙特卡罗方法的重点步骤，当不知道随机变量的概率模型服从那个分布时，可以使用均匀分布来构造；各种测量的误差、射击命中率、人的身高与体重等服从正态分布；指数分布可用在排队论与可靠性分析中；泊松分布可用在产品检验、排队系统、物理等领域中。

3. 应用举例

3.1 问题提出

某食品加工厂主要生产即食产品，一般当天生产的产品必须当天售出，否则就会出现不能保质、或变质、造成一定的经济损失，如果市场需求量大而生产量不足，则也会影响工厂的销售收入，该产品的单位成本为1.5元，单位产品售价为4元。工厂为了避免产品滞销存货过多而造成的经济损失，提出了如何制定合理的生产与库存数量的方案问题，能够使得工厂能有尽可能多的收益，经初步考虑拟从以下两种生产与库存方案中选出一个较好的方案： 方案(1)：按前一天的销售量作为当天的生产库存量。 方案(2)：按前两天的平均销售量作为当天的生产库存量。

3.2 问题分析与假设

分析可知：

  • 问题本身存在随机概率性质，其每日的销售量X就是一个随机变量，并且满足正态分布。

  • 问题的解是方案的比较结果，可以从收益率来进行比较，收益率 Y 就是 X 的函数。

  • 通过模拟一定时间的销售过程就可以得到收益率随机变量的期望。因此很好进行模拟。

假设：

  • 模拟方案为T天后的收益情况。

  • 随机变量 X 满足 N(1500, 30^2) 的正态分布。

3.3 编码模拟

要点

  • 设置3个变量，维护 cur、pre1、pre2 的销售量，表示当前的销售量，前一天、前两天的销售量
  • 多次模拟求平均
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

import numpy as np

def mcun(T, S1, S2, S11):
    LS1 = 0
    LS2 = 0
    k = 1
    while k < T:
        KC1 = S11
        KC2 = (S1 + S2) / 2
        C = np.random.normal(1500, 30)
        if C > KC1:
            ST1 = KC1
        else:
            ST1 = C
        if C > KC2:
            ST2 = KC2
        else:
            ST2 = C
        L1 = 4 * ST1 - 1.5 * KC1
        L2 = 4 * ST2 - 1.5 * KC2
        LS1 = LS1 + L1
        LS2 = LS2 + L2
        k = k + 1
        S2 = S1
        S1 = ST2
        S11 = ST1
    return LS1, LS2

    
S1 = np.random.normal(1500, 30)
S21 = np.random.normal(1500, 30)
S22 = np.random.normal(1500, 30)

ls1, ls2 = mcun(6000, S1, S21, S22)
print("first day: ", ls1)
print("second day: ", ls2)

可以得到第二种方案较好