数学建模十大算法之——蒙特卡罗

杨利霞

数学建模十大算法之——蒙特卡罗
, u* P9 r) A# T7 i& @5 b% r3 v+ |* n& ^
( v3 [' T& ^& i2 S: T

专栏文章汇总

: M, P$ {# }4 ?' g0 a) t1 V; e

文章结构如下：

• 蒙特卡罗方法及基本原理
• 蒙特卡罗方法的步骤
• 蒙特卡罗方法在数学建模中的应用举例
• 参考文献
  

---

注：以前写过一篇马尔可夫链蒙特卡罗算法（MCMC），最近在看数学建模，看到有十大算法，因此在这里重新整理一个系列。这篇文章是卫艳荣，黄瑞芳[2012]的一篇学报论文，写的不错（作为博客论文来看）就搬来了，将文中matlab程序改成了python程序。
% y* O$ M* d# w% d. U
蒙特卡罗方法是一个以概率模型为基础，利用计算机通过多次反复模拟实验完成问题求解的一种数值计算方法。它特别适用于用传统的解析法难以解决甚至是无法解决的问题。

---

1： 蒙特卡罗方法及基本原理

蒙特卡罗方法源于美国在第二次世界大战中研制原子弹的“曼哈顿计划”。该计划的主持人之一数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法，其基本思想源于法国数学家蒲丰提出著名的蒲丰投针实验，并以该方法求圆周率$\pi$。

蒲丰投针实验

如图，在平面上画出一组平行线，他们之间的距离是 ，并找一根粗细均匀长度为 的细针，然后将小针一次次地任意投掷在平面上。这样反复投掷多次，记录下针与这组平行线中任意一条相交的次数，就可以计算出 的近似值。

设针的中点与最近一条平行线间的距离为 、针与直线的夹角为 ，则样本空间为 ，则 的边长为 和 ，其面积是 。

当且仅当 成立时，针与平行线相交。满足此条件的是上图阴影区域 ，其面积为： 。

由几何概率定义，得细针与平行线相交的概率为：

记 为投针中的试验次数， 为针与平行线相交的次数，则 为细针与平行线相交的频率。由大数定理，当试验次数足够多时，可将事件发生的频率看作事件发生的概率的近似值，即： ，于是得到：

蒲丰投针试验揭示了蒙特卡罗方法的基本思想：在大数定理的保证下，利用事件发生的“频率”作为事件发生的“概率”的近似值。只要设计一个随机试验，使一个事件的概率与某未知数有关，然后通过重复试验，以频率近似值表示概率，即可求得该未知数的近似值。显然利用随机试验求近似解，试验次数要相当多才行，且次数越多近似效果也越好。此种方法可以求解微分方程，求多重积分，求特征值等。

---

2： 蒙特卡罗方法的步骤

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

• 描述或构造概率过程：对于已有的随机性质问题可描述和模拟这个概率过程，对于不具有随机性质的确定性问题，需要人为地构造一个概率过程。
• 利用概率分布抽样：通过计算机产生已知概率分布的随机变量，常用的概率分布有均匀分布，正态分布、指数分布、泊松分布等。
• 建立各种估计量：构造了随机概率模型，并从中抽样后，就要确定一个随机变量，作为所要求问题的解。一般是把 次随机抽样结果的算术平均值作为解的近似值。
  

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

---

3： 蒙特卡罗方法在数学建模中的应用举例

3.1 问题的提出

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

5 M* E7 i0 Z. l% \+ \

3.2 问题的分析与假设及参数定义

解决问题的基本思路：利用蒙特卡罗方法随机模拟市场对该产品需求量，统计计算出按照两种不同方案 天后工厂的经济值，比较不同方案经济效益的大小，选出一个较好的方案。

假设市场对该产品的每天需求数量是一个随机变量，从统计学的角度分析得知，该随机变量服从正态分布 。

& v4 O% Q, i( d6 }, x

3.3 模型的建立与求解

根据上面的分析，利用蒙特卡罗方法编程实现，主要随机模拟前一天和前两天的各种不同的销售量，来确定当天的生产与库存量，依据可能的实际销售量，计算出当天的销售利润，选择使连续几天利润尽可能大的方案。

import numpy as npdef mcun(T, S1, S21, S22):    '''    params:    T: 模拟天数    S1: 方案(1)前一天的销售量    S21: S21表示方案(2)前一天的销售量    S22: S22表示方案(2)前二天的销售量    return:    方案（1）利润；方案（2）利润    '''    LS1 = 0    # 方案(1)实际累计总利润    LS2 = 0    # 方案(2)实际累计总利润    k = 1    while k < T:        KC1 = S1                           # 方案(1)当天的生产与库存量        KC2 = (S21 + S22) / 2              # 方案(2)当天的生产与库存量        C = np.random.normal(1500, 30^2)   # 每天需求量        if C > KC1:            ST1 = KC1                      # 方案(1)当天的实际销售量        else:            ST1 = C        if C > KC2:            ST2 = KC2                      # 方案(2)当天的实际销售量        else:            ST2 = C        L1 = 4 * ST1 - 1.5 * KC1           # 方案(1)当天的实际利润           L2 = 4 * ST2 - 1.5 * KC2           # 方案(2)当天的实际利润        LS1 = LS1 + L1        LS2 = LS2 + L2        k = k + 1        S1 = ST1        S22 = S21        S21 = ST2    return LS1, LS2mcun(20, 1489, 1498, 1502)(72791.25097491541, 72674.80873238193)mcun(20, 1499, 1488, 1505)(74334.30474556636, 73467.12908653304)

模拟几次可以看到方案（1）利润更高。