机器学习回顾篇（16）：蒙特卡洛方法

蒙特卡罗（MC，Monte Carlo）方法是一种随机采样模拟求解的方法，又被称统计试验方法或者统计模拟方法。起初，蒙特卡罗方法的提出是20世纪40年代冯·诺伊曼，斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯等人为推进研制原子弹的“曼哈顿”计划而提出，但大概是因为蒙特卡罗方法是一种随机模拟的方法，与赌博场里面的扔骰子的过程十分相似而以赌城的名字命名这一方法。现如今，这一方法已被广泛应用到科学计算的众多应用领域中。

在这里，之所以没有说是蒙特卡罗算法，而是称蒙特卡罗方法，是因为任何通过生成合适的随机数，并使用这些随机树建立概率模型，以找到难以通过其他方法解决的数值问题的近似解的方法归类为蒙特卡罗方法。蒙特卡罗这一家族中还包括蒙特卡罗算法、蒙特卡罗模拟、蒙特卡罗过程、蒙特卡罗搜索树（AlphoGo中使用的方法）等众多分支。

所以，在理解蒙特卡罗方法时，需要将其当做一个大类的算法去对待，它是一种思想，只要符合了这种思想，就属于蒙特卡罗方法，不是一个固定化的数学模型，更没有特定的数学公式。只需要理解：任何通过生成合适的随机数，并使用这些随机树建立概率模型，以找到数值问题的近似解的方法都可以称为蒙特卡罗方法

蒙特卡罗方法的一个最著名应用就是求定积分。函数f(x)$f(x)$图像如下，现在要求函数f(x)$f(x)$在[a,b]$[a,b]$之间的积分，也就是阴影部分面积，即F=∫baf(x)dx$F={\int }_a^{b}f(x)dx$。

蒙特卡罗方法的解法如下。如下图所示，我们在[a,b]$[a,b]$之间取一个数x$x$，那么我们可以粗略地将f(x)⋅(b−a)$f(x)\cdot (b-a)$来估计阴影部分的面积。当然，只用一个值进行估算，结果可能过于粗糙，所以我们可以多取几个值，然后用多个结果估计值的平均值作为最终的结果。如下图所示，取4个值进行估计：

S=14(f(x1)(b−a)+f(x2)(b−a)+f(x3)(b−a)+f(x4)(b−a))=14(b−a)(f(x1)+f(x2)+f(x3)+f(x4))=14(b−a)∑i=14f(xi)$$\begin{array}{rl}S& =\frac{1}{4}(f(x_1)(b-a)+f(x_2)(b-a)+f(x_3)(b-a)+f(x_4)(b-a))\\ & =\frac{1}{4}(b-a)(f(x_1)+f(x_2)+f(x_3)+f(x_4))\\ & =\frac{1}{4}(b-a)\sum _{i=1}^{4}f(x_i)\end{array}$$

可以预见，取4个值进行估计时，结果将更加准确。可以认为，当进行估计的值个数越多时，最终结果将越接近真实值。

这个例子就是蒙特卡罗方法的典型应用：使用随机数（x$x$）建立概率模型，对最终结果（定积分）进行估计，当估计次数越多时，估计值越接近真实值。

进一步地，我们列举一些例子来说明什么是蒙特卡罗方法。

例子1：黑箱子里装有100个球，包含数量不等的红、绿、蓝三种颜色，求红色球数量。如果值抽取一次，为红球，那么我们可以认为，黑箱子里全是红球；重复抽取一万次，去一万次抽取到红球的数量的平均值最为估计结果。

例子2：求下图所示图片中白色图案占整张图片面积的比例。随机从图像中抽取一个像素，抽取n次，看n次中白色像素所在比例。

例子3：已知圆面积公式S=π⋅r2$S=\pi \cdot r^2$，求圆周率π$\pi$的值。可以这么做，构造如下图所示的面积为1的正方形，并以边长为半径画四分之一圆。之后，随机网正方形内添加n$n$个点，落在四分之一圆内点所在比例记为四分之一圆的面积，然后通过圆面积公式可估计π$\pi$的值。

例子4：报童问题。有报童购买报纸进行售卖，利润成本计算如下：

当进货量大于需求量时，利润 = 销售报纸收入 + 废纸收入 - 成本

当进货量小于需求量时，利润 = 销售保值收入 - 成本

成本 = 每份报纸进货价格（1.3元） ×$\times$ 进货量

废纸收入 = （进货量 - 需求量） - 每份废纸价格（0.2元）

每份报纸售出价格为2元。

每天报纸中可能出现重大、一般、平淡三种类型的新闻，概率分别为0.35、0.45和0.2，每种新闻类型下的销量和可能概率如下表所示。每天的进货量相等情况下，怎么实现利润最大化。

对于报童问题，难点在于就算对于同一进货量，也会因为不同新闻类型以不同概率出现不同需求量。对于这一问题，运用蒙特卡罗方法进行模拟求解，以相应的概率构造新闻类型随机事件和不同新闻事件下的需求量，然后在不同进货量下分别进行n$n$次模拟，求n$n$次模拟的平均利润。详细解法参看下方Python实现。

下面对上述的第2、3、4个例子使用Python实现蒙特卡罗方法求解。

估算π$\pi$的值

In [1]:

import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

In [13]:

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 中文字体支持
plt.style.use('ggplot')

In [24]:

def cal_pi(total):
    inner = 0
    for i in range(total):
        x = random.random()
        y = random.random()

        dis = (x**2 + y ** 2)**0.5

        if dis <= 1:
            inner += 1
    pi = 4 * inner / total
    return pi

In [37]:

pi_lst = []
for i in range(20):
    total = 2 ** i
    pi = cal_pi(total)
    pi_lst.append(pi)
fig = plt.figure(figsize=(15,3))
axes = fig.add_axes((0,0,0.8,1))
line1, = axes.plot([str(2**i) for i in range(20)], pi_lst)
line2 = axes.plot([str(2**i) for i in range(20)], [3.14 for i in range(20)])
axes.set_xlabel('蒙特卡罗试验次数')
axes.set_ylabel('π值')
plt.show()

 

从图像可以看出，随着试验次数的增多，估计值趋于平稳，越来越靠近真实的π$\pi$值。

估计白色图案所在比例

In [41]:

from PIL import Image

In [49]:

img = Image.open('data/图片.jpg')

In [93]:

def count_white(times):
    in_count = 0
    width = img.width - 1
    height = img.height -1
    for i in range(times):
        x = random.randint(0, width - 1)
        y = random.randint(0, height - 1)
        a = img.getpixel((x, y))
        if sum(img.getpixel((x, y))) > 200:  # 白色是（255， 255， 255），三通道像素值之和大于200则认为是白色，这不严谨，但有效就好
            in_count += 1
    return in_count / times

In [97]:

rat_lst = []
for i in range(20):
    times = 2 ** i
    rat = count_white(times)
    rat_lst.append(rat)
fig = plt.figure(figsize=(15,3))
axes = fig.add_axes((0,0,0.8,1))
line1, = axes.plot([i for i in range(20)], rat_lst)
axes.set_xticklabels(['2的%s次方次'%j for j in range(20)])
axes.set_xlabel('蒙特卡罗试验次数（2的x次方）')
axes.set_ylabel('白色图案占比')
plt.show()

 

报童问题求解

In [3]:

def random_envent(env_pro):
    """
    根据时间和概率，生成指定概率分布的事件组合
    """
    env_list = []
    for env, pro in env_pro:
        [env_list.append(env) for _ in range(int(pro*100))]
    return env_list

In [4]:

env_pro = [('重大', 0.35), ('一般', 0.45), ('平淡', 0.2)]
env_list = random_envent(env_pro)

In [5]:

# 重大新闻类型随机需求量
important_env_pro = [(40, 0.03), (50, 0.05), (60, 0.15), (70, 0.2), (80, 0.35), (90, 0.15), (100, 0.07)]
important_env_list = random_envent(important_env_pro)
# 一般新闻类型随机需求量
common_env_pro = [(40, 0.1), (50, 0.18), (60, 0.4), (70, 0.2), (80, 0.08), (90, 0.04), (100, 0.0)]
common_env_list = random_envent(common_env_pro)
# 平淡新闻类型随机需求量
dull_env_pro = [(40, 0.44), (50, 0.22), (60, 0.16), (70, 0.12), (80, 0.06), (90, 0.0), (100, 0.0)]
dull_env_list = random_envent(dull_env_pro)

In [6]:

def mk_requirements(news_type):
    """
    根据指定新闻类型，随机生成销量
    """
    if news_type == '重大':
        # 重大新闻类型随机需求量
        important_env_pro = [(40, 0.03), (50, 0.05), (60, 0.15), (70, 0.2), (80, 0.35), (90, 0.15), (100, 0.07)]
        important_env_list = random_envent(important_env_pro)
        return random.choice(important_env_list)
    elif news_type == '一般':
        # 重大新闻类型随机需求量
        common_env_pro = [(40, 0.1), (50, 0.18), (60, 0.4), (70, 0.2), (80, 0.08), (90, 0.04), (100, 0.0)]
        common_env_list = random_envent(common_env_pro)
        return random.choice(common_env_list)
    else:
        dull_env_pro = [(40, 0.44), (50, 0.22), (60, 0.16), (70, 0.12), (80, 0.06), (90, 0.0), (100, 0.0)]
        dull_env_list = random_envent(dull_env_pro)
        return random.choice(dull_env_list)

In [7]:

def cal_profit(purchase, requirements):
    """
    根据进货量和需求量计算利润
    """
    cost = purchase * 1.3  # 成本
    if purchase > requirements:  # 当进货量大于需求量时
        return requirements * 2 + (purchase - requirements) * 0.2 -cost
    else:
        return purchase * 2 -cost

In [8]:

def cal_profit_10000(purchase):
    """
    对指定进货量进行10000的实验，返回平均利润
    """
    total_profit = 0
    for _ in range(10000):
        news_type = random.choice(env_list) # 按照指定概率随机生产新闻事件
        requirements = mk_requirements(news_type)  # 根据新闻类型，按照指定概率生成需求量
        total_profit += cal_profit(purchase, requirements)
    return total_profit / 10000

In [38]:

def mk_pic():
    x = [i for i in range(100)]
    y = [cal_profit_10000(i) for i in range(100)]
    fig = plt.figure(figsize=(15,3))
    axes = fig.add_axes((0,0,0.8,1))
    axes.plot(x, y)
    axes.set_xlabel('进货量')
    axes.set_ylabel('1w次平均利润')
    plt.show()

In [39]:

mk_pic()

 

从图像可知，60次左右的进货量可获得最大利润。