Monte Carlo Algorithms. 蒙特卡洛算法是一大类随机算法，又称为随机抽样或统计试验方法，通过随机样本估计真实值。

下面用几个实例来理解蒙特卡洛算法。

6. 蒙特卡洛算法

6.1 计算 \(\pi\)

a. 原理

如果我们不知道 \(\pi\) 的值，我们能不能用随机数 来近似 \(\pi\) 呢？

假设我们用一个随机数生成器，每次生成两个范围在 \([-1,+1]\) 的随机数，一个作为 x，另一个作为 y，即生成了一个二维随机点：

假如生成 1亿 个随机样本，会有多少落在 半径=1 的圆内？这个概率就是圆的面积除以正方形的面积。

即：\(P = \frac{\pi{r^2}}{2^2}=\frac{\pi}{4}\)

假设从正方形区域中随机抽样 n 个点，那么落在圆内点个数的期望为：\(P_n=\frac{\pi{n}}{4}\)，

下面我们去求落在圆内的点的个数，只需满足\(x^2+y^2\leqslant1\) 即为圆内。

如果生成的随机点的个数足够多，落在圆内的实际观测值 \(m\approx \frac{\pi{n}}{4}\)；

我们已知了m 与 n，所以\(\pi \approx \frac{4m}{n}\).

事实上，根据概率论大数定律：

\(\frac{4m}{n}\rightarrow \pi\)，as n → ∞

这保证了蒙特卡洛的正确性。

伯恩斯坦概率不等式还能确定 观测值和真实值之间误差的上界。

\(|\frac{4m}{n}-\pi|=O(\frac{1}{\sqrt{n}})\)

说明 这个误差与样本n的根号成反比。

b. 代码

下面放一个Python代码

#coding=utf-8
#蒙特卡罗方法计算 pi
import random,math,time
start_time = time.perf_counter()
s = 1000*1000
hits = 0
for i in range(s):
    x = random.random()
    y = random.random()
    z = math.sqrt(x**2+y**2)
    if z<=1:
        hits +=1

PI = 4*(hits/s)
print(PI)
end_time = time.perf_counter()
print("{:.2f}S".format(end_time-start_time))

# 输出
3.141212
0.89S

另外可还有一个可视化程序，可以模拟点落在方块区域圆内外：http://www.anders.wang/monte-carlo/

6.2 Buffon's Needle Problem

a. 原理

布封投针，也是用蒙特卡洛来近似 \(\pi\) 值。这是一个可以动手做的实验。

用一张纸，画若干等距平行线（距离为 d），撒上一把等长的针（长度为l），通过与平行线相交的针的数量，就可以推算出 \(\pi\)。

通过微积分可以算出：相交的概率为：\(P = \frac{2l}{\pi{d}}\)

微积分推导过程：

课程里并没有讲解推导，这里我参考的是一下两篇博客的推导过程：

1. https://zhuanlan.zhihu.com/p/479953215
2. https://cosx.org/2009/11/a-brief-talk-on-buffon-throwing-needle-problems/

主流做法是通过对针的斜率进行积分：

这里我后续补充。

跟 6.1 类似，我们随机扔 n 根针，这样相交个数的期望为 \(Pn = \frac{2ln}{\pi{d}}\) 。我们可以观察到（如果是电脑模拟即为通过公式判断出）有 m 跟针实际与线相交，如果n足够大，则 \(m\approx \frac{2ln}{\pi{d}}\)。

求 \(\pi\) 公式即为： \(\pi\approx \frac{2ln}{md}\)

b. 代码

有了公式 \(\pi\approx \frac{2ln}{md}\)，代码实现其实很简单了，仅列出一种实现思路：

import numpy as np

def buffon(a,l,n):
  xl = np.pi*np.random.random(n)
  yl = 0.5*a*np.random.random(n)
  m = 0
  for x,y in zip(xl,yl):
    if y < 0.5*l*np.sin(x):
      m+=1
  result = 2*l/a*n/m
  print(f'pi的估计值是{result}')
  
buffon(2,1,1000000)

# 输出为：
pi的估计值是3.153977165205324

当然，也有可视化的代码：

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

NUMBER_OF_NEEDLES = 5000


class DefineNeedle:
    def __init__(self, x=None, y=None, theta=None, length=0.5):
        if x is None:
            x = random.uniform(0, 1)
        if y is None:
            y = random.uniform(0, 1)
        if theta is None:
            theta = random.uniform(0, math.pi)

        self.needle_coordinates = np.array([x, y])
        self.complex_representation = np.array(
            [length/2 * math.cos(theta), length/2*math.sin(theta)])
        self.end_points = np.array([np.add(self.needle_coordinates, -1*np.array(
            self.complex_representation)), np.add(self.needle_coordinates, self.complex_representation)])

    def intersects_with_y(self, y):
        return self.end_points[0][1] < y and self.end_points[1][1] > y


class BuffonSimulation:
    def __init__(self):
        self.floor = []
        self.boards = 2
        self.list_of_needle_objects = []
        self.number_of_intersections = 0

        fig = plt.figure(figsize=(10, 10))
        self.buffon = plt.subplot()
        self.results_text = fig.text(
            0, 0, self.estimate_pi(), size=15)
        self.buffon.set_xlim(-0.1, 1.1)
        self.buffon.set_ylim(-0.1, 1.1)

    def plot_floor_boards(self):
        for j in range(self.boards):
            self.floor.append(0+j)
            self.buffon.hlines(
                y=self.floor[j], xmin=0, xmax=1, color='black', linestyle='--', linewidth=2.0)

    def toss_needles(self):
        needle_object = DefineNeedle()
        self.list_of_needle_objects.append(needle_object)
        x_coordinates = [needle_object.end_points[0]
                         [0], needle_object.end_points[1][0]]
        y_coordinates = [needle_object.end_points[0]
                         [1], needle_object.end_points[1][1]]

        for board in range(self.boards):
            if needle_object.intersects_with_y(self.floor[board]):
                self.number_of_intersections += 1
                self.buffon.plot(x_coordinates, y_coordinates,
                                 color='green', linewidth=1)
                return
        self.buffon.plot(x_coordinates, y_coordinates,
                         color='red', linewidth=1)

    def estimate_pi(self, needles_tossed=0):
        if self.number_of_intersections == 0:
            estimated_pi = 0
        else:
            estimated_pi = (needles_tossed) / \
                (1 * self.number_of_intersections)
        error = abs(((math.pi - estimated_pi)/math.pi)*100)
        return (" Intersections:" + str(self.number_of_intersections) +
                "\n Total Needles: " + str(needles_tossed) +
                "\n Approximation of Pi: " + str(estimated_pi) +
                "\n Error: " + str(error) + "%")

    def plot_needles(self):
        for needle in range(NUMBER_OF_NEEDLES):
            self.toss_needles()
            self.results_text.set_text(self.estimate_pi(needle+1))
            if (needle+1) % 200 == 0:
                plt.pause(1/200)
        plt.title("Estimation of Pi using Probability")

    def plot(self):
        self.plot_floor_boards()
        self.plot_needles()
        plt.show()


simulation = BuffonSimulation()
simulation.plot()

效果如图：

以上内容参考：

1. 课程视频
2. https://www.section.io/engineering-education/buffon-needle/
3. https://github.com/topics/buffon-needle
4. https://github.com/GunnarDahm/buffon_monte_carlo_sim/blob/master/buffon_monte_carlo.py
5. https://blog.csdn.net/qq_45757739/article/details/108387567
6. https://blog.csdn.net/TSzero/article/details/111604960

理解思想即可，如果后续有机会，可能单出一篇介绍介绍，也有可能将这部分丰富一下。

6.3 估计阴影部分的面积

我们稍微推广一下，试着用蒙特卡洛解决一个阴影部分面积的求解。比如下图：

我们如何使用蒙特卡洛的思路解决这个阴影部分面积的求解呢？

类似于上面的思路，在正方形内做随机均匀抽样，得到很多点，怎么确定点在阴影部分呢？

可知，阴影部分的点满足：

• 易知，正方形面积 \(A_1=4\)；设阴影部分面积为 \(A_2\)
• 随机抽样的点落在阴影部分的概率为：\(P=\frac{A_2}{A_1}=\frac{A_2}{4}\)
• 从正方形区域抽样 n 个点，n尽可能大，则来自阴影部分点的期望为：\(nP=\frac{nA_2}{4}\)；
• 如果实际上满足上述条件的点 有 m 个，则令 \(m\approx nP\)
• 得到：\(A_2\approx \frac{4m}{n}\)

代码与 6.1 相近。

6.4 求不规则积分

近似求积分是蒙特卡洛在工程和科学问题中最重要的应用。很多积分是没有解析的积分（即可以计算出来的积分），特别是多元积分，而只能用数值方法求一个近似值，蒙特卡洛就是最常用的数值方法。

一元函数步骤如下：

我们要计算一个一元函数的定积分 \(I = \int_a^bf(x)dx\);

• 从区间 \([a,b]\) 上随机均匀抽样 \(x_1,x_2,...,x_n\);

• 计算 \(Q_n = (b-a)\frac{1}{n}\sum_{i=1}^nf(x_i)\)，即均值乘以区间长度；

  这里均值乘以区间长度是 实际值，而 I 是期望值

• 用 \(Q_n\) 近似 \(I\)

大数定律保证了 当\(n\rightarrow\infty,Q_n\rightarrow I\)

多元函数步骤如下：

我们要计算一个多元函数的定积分 \(I = \int_a^bf(\vec{x})d\vec{x}\)，积分区域为 \(\Omega\);

• 从区间 \(\Omega\) 上随机均匀抽样 \(\vec{x_1},\vec{x_2},...,\vec{x_n}\);

• 计算 \(\Omega\) 的体积V（高于三维同样）：\(V=\int_\Omega{d\vec{x}}\)；

  hh值得注意的是，这一步仍要计算定积分，如果形状过于复杂，无法求得 V，那么无法继续进行，则无法使用蒙特卡洛算法。所以只能适用于比较规则的区域，比如圆形，长方体等。

• 计算 \(Q_n =V \frac{1}{n}\sum_{i=1}^nf(\vec{x_i})\)，即均值乘以区间长度；

  这里均值乘以区间长度是 实际值，而 I 是期望值

• 用 \(Q_n\) 近似 \(I\)

下面我们从积分的角度再来看看 蒙特卡洛近似求 pi

• 定义一个二元函数 \(f(x,y)=\begin{cases} 1 \ \ if点在圆内\\ 0 \ \ if 点在圆外\end{cases}\);
• 定义一个区间 \(\Omega=[-1,1]×[-1,1]\)
• \(I =\pi {r^2}=\pi\)
• 接下来用蒙特卡洛近似 I，得到关于 \(\pi\)的算式即可得到近似的\(\pi\);
  • 随机抽样 n 个点，记为\((x_1,y_1),...,(x_n,y_n)\)
  • 计算 区域面积 \(V = \int_\Omega{dxdy}=4\)；
  • 计算 \(Q_n =V \frac{1}{n}\sum_{i=1}^nf(x_i,y_i)\)
  • 蒙特卡洛近似 Q 与 I 近似相等：\(\pi=Q_n=\int_\Omega{f(x,y)}{dxdy}\)

这是从蒙特卡洛积分的角度得到的pi，6.1 中则是从蒙特卡洛概率和期望的角度得到的。

6.5 用蒙特卡洛近似期望

这个方法对于统计学和机器学习很有用。

• 定义 X 是 d 维的随机变量，函数 p(x) 是一个PDF，概率密度函数；
• 函数 \(f(x)\) 的期望：\(\mathbb{E}_{x\sim{p}}[f(X)]=\int_{R^d}f(X)\cdotp(x)dx\)
• 直接以上面的方式求期望可能并不容易，所以通常使用蒙特卡洛近似求期望：
  1. 随机抽样：根据概率密度函数 \(p(x)\) 进行随机抽样，记为\(X_1,X_2,...,X_n\)；
  2. 计算 \(Q_n =\frac{1}{n}\sum_{i=1}^nf(x_i)\)
  3. 用 Q 近似 期望\(\mathbb{E}_{x\sim{p}}[f(X)]\)

6.6 总结 | 蒙特卡洛算法的思想

我的想法是尽量精简，即：

模拟---抽样---估值，通过模拟出来的大量样本集或者随机过程，以随机抽样的方式，去近似我们想要研究的实际问题对象。

补充蒙特卡洛相关：

• 蒙特卡洛是摩洛哥的赌场；

• 蒙特卡洛算法得到的结果通常是错误的，但很接近真实值，对于对精度要求不高的机器学习已经足够。

  随机梯度下降就是一种蒙特卡洛算法，用随机的梯度近似真实的梯度，不准确但是降低了计算量。

• 蒙特卡洛是一类随机算法，除此以外还有很多随机算法，比如拉斯维加斯算法（结果总是正确的算法）

x. 参考教程

• 视频课程：深度强化学习（全）_哔哩哔哩_bilibili
• 视频原地址：https://www.youtube.com/user/wsszju
• 课件地址：https://github.com/wangshusen/DeepLearning