機器學習回顧篇（16）：蒙特卡洛演算法
蒙特卡羅（MC，Monte Carlo）方法是一種隨機取樣模擬求解的方法，又被稱統計試驗方法或者統計模擬方法。起初，蒙特卡羅方法的提出是20世紀40年代馮·諾伊曼，斯塔尼斯拉夫·烏拉姆和尼古拉斯·梅特羅波利斯等人為推進研製原子彈的「曼哈頓」計劃而提出，但大概是因為蒙特卡羅方法是一種隨機模擬的方法，與賭博場裡面的扔骰子的過程十分相似而以賭城的名字命名這一方法。現如今，這一方法已被廣泛應用到科學計算的眾多應用領域中。
在這裡，之所以沒有說是蒙特卡羅演算法，而是稱蒙特卡羅方法，是因為任何通過生成合適的隨機數，並使用這些隨機樹建立概率模型，以找到難以通過其他方法解決的數值問題的近似解的方法歸類為蒙特卡羅方法。蒙特卡羅這一家族中還包括蒙特卡羅演算法、蒙特卡羅模擬、蒙特卡羅過程、蒙特卡羅搜索樹（AlphoGo中使用的方法）等眾多分支。
所以，在理解蒙特卡羅方法時，需要將其當做一個大類的演算法去對待，它是一種思想，只要符合了這種思想，就屬於蒙特卡羅方法，不是一個固定化的數學模型，更沒有特定的數學公式。只需要理解：任何通過生成合適的隨機數，並使用這些隨機樹建立概率模型，以找到數值問題的近似解的方法都可以稱為蒙特卡羅方法
蒙特卡羅方法的一個最著名應用就是求定積分。函數$f(x)$影像如下，現在要求函數$f(x)$在$[a, b]$之間的積分，也就是陰影部分面積，即$F = \int_a^b {f(x)dx} $。
蒙特卡羅方法的解法如下。如下圖所示，我們在$[a, b]$之間取一個數$x$，那麼我們可以粗略地將$f(x) \cdot (b – a)$來估計陰影部分的面積。當然，只用一個值進行估算，結果可能過於粗糙，所以我們可以多取幾個值，然後用多個結果估計值的平均值作為最終的結果。如下圖所示，取4個值進行估計：
可以預見，取4個值進行估計時，結果將更加準確。可以認為，當進行估計的值個數越多時，最終結果將越接近真實值。
這個例子就是蒙特卡羅方法的典型應用：使用隨機數（$x$）建立概率模型，對最終結果（定積分）進行估計，當估計次數越多時，估計值越接近真實值。
進一步地，我們列舉一些例子來說明什麼是蒙特卡羅方法。
例子1：黑箱子里裝有100個球，包含數量不等的紅、綠、藍三種顏色，求紅色球數量。如果值抽取一次，為紅球，那麼我們可以認為，黑箱子里全是紅球；重複抽取一萬次，去一萬次抽取到紅球的數量的平均值最為估計結果。
例子2：求下圖所示圖片中白色圖案占整張圖片面積的比例。隨機從影像中抽取一個像素，抽取n次，看n次中白色像素所在比例。
例子3：已知圓面積公式$S = \pi \cdot {r^2}$，求圓周率$π$的值。可以這麼做，構造如下圖所示的面積為1的正方形，並以邊長為半徑畫四分之一圓。之後，隨機網正方形內添加$n$個點，落在四分之一圓內點所在比例記為四分之一圓的面積，然後通過圓面積公式可估計$π$的值。
例子4：報童問題。有報童購買報紙進行售賣，利潤成本計算如下：
當進貨量大於需求量時，利潤 = 銷售報紙收入 + 廢紙收入 – 成本
當進貨量小於需求量時，利潤 = 銷售保值收入 – 成本
成本 = 每份報紙進貨價格（1.3元） $\times $ 進貨量
廢紙收入 = （進貨量 – 需求量） – 每份廢紙價格（0.2元）
每份報紙售出價格為2元。
每天報紙中可能出現重大、一般、平淡三種類型的新聞，概率分別為0.35、0.45和0.2，每種新聞類型下的銷量和可能概率如下表所示。每天的進貨量相等情況下，怎麼實現利潤最大化。
對於報童問題，難點在於就算對於同一進貨量，也會因為不同新聞類型以不同概率出現不同需求量。對於這一問題，運用蒙特卡羅方法進行模擬求解，以相應的概率構造新聞類型隨機事件和不同新聞事件下的需求量，然後在不同進貨量下分別進行$n$次模擬，求$n$次模擬的平均利潤。詳細解法參看下方Python實現。
下面對上述的第2、3、4個例子使用Python實現蒙特卡羅方法求解。
估算$π$的值
import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 中文字體支援
plt.style.use('ggplot')
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
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()
從影像可以看出，隨著試驗次數的增多，估計值趨於平穩，越來越靠近真實的$π$值。
估計白色圖案所在比例
from PIL import Image
img = Image.open('data/圖片.jpg')
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
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()
報童問題求解
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
env_pro = [('重大', 0.35), ('一般', 0.45), ('平淡', 0.2)]
env_list = random_envent(env_pro)
# 重大新聞類型隨機需求量
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)
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)
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
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
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()
mk_pic()
從影像可知，60次左右的進貨量可獲得最大利潤。
