掷镖游戏和约算Pi

你和你的几个朋友在一起玩飞镖。你们对这个游戏都非常不会玩。你们丢的每一镖一定会击中以下的正方形板，但是除此之外每次丢的镖会随机落到正方形中的任何位置。为了在展现你垃圾技巧的同时自娱自乐，你决定这是一个约算无理数 \(\pi \approx 3.14159\) 的好机会。

因为你的每一镖都会随机落到正方形中，你发现镖落到圆的几率等于圆面积相比正方形面积的比例：

\begin{equation} P_{circle} = \frac{Area_{circle}}{Area_{square}} = \frac{\pi r^2}{(2r)^2} \end{equation}

更进一步，我们可以通过镖落入圆形的比例来约算 \(P_{circle}\)。因此，我们写出以下公式：

\begin{equation} \frac{N_{circle}}{N_{total}} \approx \frac{\pi r^2}{(2r)^2} = \frac{\pi}{4} \end{equation}

其中 \(N_{total}\) 是总共丢的飞镖数，\(N_{circle}\) 是落入圆中的飞镖数。因此，通过计算飞镖落到的不同位置的次数，你可以开始约算 \(\pi\) 的值！

问题1

编写代码来模拟丢掷飞镖和数其是否落入圆中的代码，并“在飞镖被丢掷”的同时计算实时的 \(\pi\) 值。为了简单性，你可以假设板的中心在 \((0, 0)\)，且 \(r = 1\)（圆的半径）。使用 numpy.random.rand（说明文档链接）来随机生成飞镖落到板上的位置。为总共 \(N = 10,000\) 个镖进行以上计算。为每个丢出的镖决定其是否落入了圆中，并根据以下公式更新你对 \(\pi\) 的估测值：\(N_{circle} / N_{total} \approx \pi / 4\)

请记住，每个镖可以落入的范围是 \((x \in [-1, 1], y \in [-1, 1])\)，而一个落到 \((x, y)\) 的镖只有满足以下条件才算是落入到了圆中：

\begin{equation} \sqrt{x^2 + y^2} < 1 \end{equation}

你可以首先通过显式for循环来编写一个直观的解。虽然如此，你应该努力编写一个完全矢量化的解（也就是说，不使用任何显式for循环来计算 \(10,000\) 次丢镖中每丢一镖后的 \(\pi\) 约算值）。

相关阅读

你会需要熟悉NumPy的矢量化操作和根据轴求和。你也很可能会发现布尔索引很有用。

这里强烈建议你使用matplotlib来可视化你对 \(\pi\) 每一轮后的约算值。你可以参阅PLYMI的这一小节来学习如何创建一个简单的线图。

提示

了解NumPy的累计和函数 numpy.cumsum 会很有用。这对计算每一轮后的总数很有用——也就是说，这可以帮助你计算每一轮后总共落入圆中的镖的数量，而不仅仅是这一轮是否在圆中。

解（未矢量化）

首先，我们想要生成 \(N = 10,000\) 镖的 \((x, y)\) 坐标。使用 numpy.random.rand(N, 2)，我们可以生成有 \(N\) 行的2维数组——每一行包含着一镖的 \((x, y)\) 坐标。

我们想要每一镖的 \(x\) 和 \(y\) 坐标落入 \([-1, 1]\)。numpy.random.rand 生成在范围 \([0, 1)\) 中的数。我们可以对生成的数乘以2并减去1来使得生成树随机落入 \([-1, 1)\)。

import numpy as np
N = 10000
dart_positions = 2 * np.random.rand(N, 2) - 1  # 生成 [-1, 1] 中的数字

现在，我们可以循环迭代所有的镖位置，决定其是否在圆中，并根据结果更新 \(N_{circle}\)。

Ncircle = [0] # 将首位值设为0来简化循环中的逻辑

for x,y in dart_positions:
    if np.sqrt(x**2 + y**2) < 1:
        Ncircle.append(Ncircle[-1] + 1) # 镖落入了圆中
    else:
        Ncircle.append(Ncircle[-1])  # 镖落入了圆外——Ncircle没有变化

现在，让我们使用我们的列表 Ncircle 来计算每一轮后对 \(\pi\) 的约算值。

running_estimate = []

for n_total, n_circle in enumerate(Ncircle[1:]): # 跳过第一个的零值
    # n_total将从0开始，所以我们需要加1
    running_estimate.append(4 * n_circle / (n_total + 1))

让我们打印我们对前10镖后和最后10镖后的约算值。我们应该能注意在前10镖后的约值非常不精确且受噪声影响，而在最后10镖后明显更加紧密。

running_estimate[:10]

[4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0]

running_estimate[-10:]

[3.1380242217996197,
 3.1381104883907125,
 3.1381967377164015,
 3.1378827296377825,
 3.1379689844922463,
 3.1380552220888354,
 3.137741322396719,
 3.1378275655131027,
 3.137913791379138,
 3.138]

解（矢量化）

首先，我们想要生成 \(N = 10,000\) 镖的 \((x, y)\) 坐标。使用 numpy.random.rand(N, 2)，我们可以生成有 \(N\) 行的2维数组——每一行包含着一镖的 \((x, y)\) 坐标。

我们想要每一镖的 \(x\) 和 \(y\) 坐标落入 \([-1, 1]\)。numpy.random.rand 生成在范围 \([0, 1)\) 中的数。我们可以对生成的数乘以2并减去1来使得生成树随机落入 \([-1, 1)\)。

import numpy as np
N = 10000
dart_positions = 2 * np.random.rand(N, 2) - 1  # 生成在 [-1, 1] 中的数字

现在，我们将要为 dart_positions 中的每一镖计算和原点的距离 \(\sqrt{x^2 + y^2}\)。我们可以平方 dart_positions 中的每一个值并顺着它的列（轴1）求和，然后平方根其结果。这将产生一个形状为 \((N,)\) 的储存每一镖和原点距离的数组。

dist_from_origin = np.sqrt((dart_positions**2).sum(axis=1))  # 形状为 (N,) 的数组

你也可以使用内置的 np.linalg.norm 来编写更加简短的等值代码。

现在，我们想要确认这些镖是否落入圆中。也就是说，我们想要求 \(\sqrt{x^2 + y^2} < 1\) 是否为真。我们可以直接使用 < 来对每个成员进行对比。这将返回一个布尔值数组，其值在镖落入圆中时为 True 而不落入圆中时为 False。

is_in_circle = dist_from_origin < 1  # 形状为 (N,) 的布尔数组

最后，我们想要计算每一轮后落入圆中的总镖数。请回忆，True 的行为类似 1，而 False 的行为类似 0。因此，我们可以对 in_circle 进行累计和来计算这个。

# 累计和：num_in_circle[i] = sum(is_in_circle[:i])
num_in_circle = np.cumsum(is_in_circle)

num_thrown = np.arange(1, N+1) # 1, 2, ..., N

最后，我们可以通过 \(N_{circle} / N_{total} \approx \pi / 4\) 计算每一镖后对 \(\pi\) 的预估值。

running_estimate = 4 * num_in_circle / num_thrown

让我们通过绘制结果来查看它们。我们将会创建一个简单的线图并将pi的正确值作为虚的横线画出。因为我们丢了那么多镖，将镖数绘制在一个对数标尺（log scale）中会更加直观。这将允许我们查看我们的预估是如何在丢几十次和几百次和几千次等等后提升的。

import matplotlib.pyplot as plt
%matplotlib notebook

fig, ax = plt.subplots()

ax.plot(num_thrown, running_estimate);
ax.hlines(y=np.pi, xmin=1, xmax=N+1, linestyles="--")  # 代表pi真实值的横线
ax.set_xscale("log")

# 译者注：无法翻译轴标注，因为matplotlib不支持中文字体
ax.set_ylabel(r"Estimated value of $\pi$")
ax.set_xlabel("Number of randomly-thrown darts (log-scale)")
ax.grid(True)

问题2

尝试多次重新运行并绘制你的解。你会发现你的约算值的曲线在每次运行之间有着很大的区别，但因为投镖是随机的所以也是可以理解的。虽然如此，这些曲线应该每次都应该在更多镖被丢出后接近 \(\pi\) 的正确值。

让我们来学习一下这个过程中的统计学知识。让我们模拟在 \(M = 100\) 次独立测试中丢出 \(N = 10,000\) 个镖的结果。为每一镖，计算在 \(M\) 次测试中 \(\pi\) 预估的平均值。绘制这个平均值曲线以及 平均值 + 标准差 的上限和 平均值 - 标准差 的下限曲线。你可以使用ax.fill_between来涂色这两个极限中的面积。

如果你没有觉得上一题的矢量化解对比未矢量化解更加优雅，那么你或许将在这一题中体会到矢量化的强大。

解

与其生成 \((N, 2)\) 个镖位置，我们将生成 \((M, N, 2)\) 个位置。也就是说，轴0代表着 \(N=10,000\) 镖的各次测试，轴1代表着每一镖，且轴2代表着镖落到的 \((x, y)\) 坐标。在如此整理我们的镖坐标后，将之前的矢量化解拓展到 \(M\) 次测试就很简单了。

import numpy as np

M = 100
N = 10000

dart_positions = np.random.rand(M, N, 2) * 2 - 1               # 形状为 (M, N, 2) 的位置数组
dist_from_origin = np.sqrt((dart_positions**2).sum(axis=2))    # 形状为 (M, N) 的距离数组
is_in_circle = dist_from_origin <= 1                           # 形状为 (M, N) 的布尔数组

num_thrown = np.arange(1, N+1)  # 1, 2, ..., N, 形状为 (N,)
num_in_circle = np.cumsum(is_in_circle, axis=1)  # 形状为 (M, N)

# 广播相除来产生M次测试中的pi约值
running_estimate = 4 * num_in_circle / num_thrown

当我们有了 \((M, N)\) 个 \(\pi\) 约值——\(N\) 个镖在 \(M\) 次计算中的约算值——后，我们可以简单地计算 \(M\) 次测试的平均值和标准差。

mean_in_circled = running_estimate.mean(axis=0)  # 不同测试的平均值
std_in_circle = running_estimate.std(axis=0)     # 不同测试的标准差

想象一下使用for循环会有多麻烦。矢量化使得这一额外维度的分析仅仅成为第一题解的简单拓展。

我们现在可以可视化我们每一轮后的平均值和标准差。ax.fill_between 提供了涂色平均值上下标准差的面积。请注意，我们提供了参数 alpha=0.2 来使得这个涂色区域为半透明。

import matplotlib.pyplot as plt
%matplotlib notebook

fig, ax = plt.subplots()

ax.plot(num_thrown, mean_in_circled, label="mean");
ax.fill_between(num_thrown, y1=mean_in_circled-std_in_circle, y2=mean_in_circled+std_in_circle,
                alpha=0.2, label="standard deviation")
ax.hlines(y=np.pi, xmin=1, xmax=N+1, linestyles="--")

ax.set_xscale("log")
ax.grid(True)
ax.set_ylabel("Estimated value of pi")
ax.set_xlabel("Number of darts thrown")
ax.legend();

和之前猜测的一样，我们的约算值在投的镖越多后约接近正确值。而标准差为我们提供了一个约算需要投掷多少镖才能达到某个精度的方法。

我希望你喜欢这个有趣的虚拟实验并为在NumPy中运行了完全矢量化的数字模拟而自豪！