模拟

模拟是使用计算机模仿物理实验的过程。在本课程中,这些实验几乎总是涉及随机性。

我们已经看到了如何模拟抛硬币的结果。该模拟中的步骤是我们将在本课程中进行的每一次模拟的构成步骤的示例。在本节中,我们将阐述这些步骤并在示例中加以遵循。

过程

第 1 步:确定模拟什么

决定要模拟哪个量。例如,你可能决定要模拟抛硬币的结果。每个模拟值将是正面或反面。

第 2 步:模拟一个值

找出如何模拟你在第 1 步中指定的量的“一个”值。在我们的示例中,你必须找出如何模拟“一次”抛硬币的结果。如果你的量更复杂,可能需要多行代码来生成一个模拟值。通常,我们将定义一个返回模拟值的函数。

第 3 步:重复次数

决定你想模拟该量多少次。你需要重复第 2 步中的模拟那么多次。在我们之前的一个示例中,我们决定模拟 1000 次抛硬币的结果,因此我们需要重复生成单次抛掷结果 1000 次。

第 4 步:模拟多个值

最后,按以下方式将所有步骤整合起来。

  • 创建一个空数组用于收集所有模拟值。我们称其为收集数组。
  • 创建一个“重复序列”,即一个长度为你第 3 步中指定的重复次数的序列。对于 n 次重复,我们几乎总是使用序列 np.arange(n)
  • 创建一个 for 循环。对于重复序列中的每个元素:
    • 使用第 2 步中编写的函数模拟“一个”值。
    • 将该模拟值追加到收集数组中。

通常的 Python 约定是将数组放置在表格中以便于分析。

[In ]:
from datascience import *
path_data = '../../../assets/data/'
import matplotlib
matplotlib.use('Agg')
%matplotlib inline
import matplotlib.pyplot as plots
plots.style.use('fivethirtyeight')
import numpy as np

示例:100 次抛掷中的正面次数

很自然地预期,在 100 次抛硬币中,将会有 50 次正面,有增有减。

但“一些”是多少?恰好得到 50 次正面的概率是多少?诸如此类的问题在数据科学中很重要,不仅因为它们涉及随机性的有趣方面,还因为它们可以用于分析实验中治疗组和对照组的分配是由抛硬币决定的。

在这个示例中,我们将模拟 100 次抛硬币中的正面次数。结果的直方图将帮助我们了解可能的正面次数。

让我们开始模拟,按照上述步骤进行。

第 1 步:确定模拟什么

我们要模拟的量是 100 次抛掷中的正面次数。

第 2 步:模拟一个值

我们必须找出如何进行一组 100 次抛掷并统计正面次数。让我们从创建一枚硬币开始。

[In ]:
coin = make_array('Heads', 'Tails')

在我们之前的示例中,我们使用了 np.random.choice 和一个 for 循环来生成多次抛掷。但在数据科学中,抛硬币集经常需要用到,以至于如果我们包含第二个参数(即抛掷次数),np.random.choice 会直接为我们模拟。

以下是 10 次抛掷的结果。

[In ]:
ten_tosses = np.random.choice(coin, 10)
ten_tosses
array(['Heads', 'Tails', 'Heads', 'Tails', 'Tails', 'Heads', 'Tails',
       'Tails', 'Tails', 'Heads'], dtype='<U5')

我们可以像之前一样使用 np.count_nonzero 来统计正面次数:

[In ]:
np.count_nonzero(ten_tosses == 'Heads')
4

我们的目标是模拟 100 次抛掷中的正面次数,而不是 10 次。为此,我们可以重复相同的代码,将 10 替换为 100。

[In ]:
outcomes = np.random.choice(coin, 100)
num_heads = np.count_nonzero(outcomes == 'Heads')
num_heads
46

由于我们想多次执行此操作,让我们定义一个返回正面次数的模拟值的函数。我们可以使用上面单元格中开发的代码来实现。

[In ]:
def one_simulated_value():
    outcomes = np.random.choice(coin, 100)
    return np.count_nonzero(outcomes == 'Heads')

第 3 步:重复次数

我们将使用多少次重复由我们自己决定。使用越多,我们的模拟就越可靠,但运行代码所需的时间也越长。Python 抛硬币的速度相当快,所以让我们进行 20,000 次重复。这意味着我们将执行以下操作 20,000 次: - 抛硬币 100 次并统计正面次数。

这真是大量的抛掷!幸好我们有 Python 替我们完成。

第 4 步:模拟多个值

我们准备创建一个包含 20,000 个模拟值的数组,每个值代表 100 次抛硬币中的正面次数。

[In ]:
num_repetitions = 20000   # number of repetitions

heads = make_array() # empty collection array

for i in np.arange(num_repetitions):   # repeat the process num_repetitions times
    new_value = one_simulated_value()  # simulate one value using the function defined
    heads = np.append(heads, new_value) # augment the collection array with the simulated value

# That's it! The simulation is done.

检查数组 heads 是否包含 20,000 个条目,每个条目对应实验的一次重复。

[In ]:
len(heads)
20000

为了了解 100 次抛掷中正面次数的波动性,我们可以将结果收集到表格中并绘制直方图。

[In ]:
simulation_results = Table().with_columns(
    'Repetition', np.arange(1, num_repetitions + 1),
    'Number of Heads', heads
)
[In ]:
simulation_results.show(3)
<IPython.core.display.HTML object>
[In ]:
simulation_results.hist('Number of Heads', bins = np.arange(30.5, 69.6, 1))
Histogram with 'Number of Heads' on the x-axis and 'Percent per unit' on the y-axis. The x-axis is labeled from 30 to 70 and the data appears in a bell shape centered around 50. The tallest bar is a litle above 8. The bars have a width of 1.

每个箱宽度为 1,并以每個正面数值为中心。

毫不奇怪,直方图大致以 50 次正面为中心对称。50 处的条形高度约为每单位 8%。由于每个箱宽度为 1 个单位,这相当于说大约 8% 的重复恰好产生了 50 次正面。这个百分比不算很大,但它比其他每个正面数值的百分比都要大。

直方图还显示,在几乎所有重复中,100 次抛掷中的正面次数都在 35 到 65 之间。实际上,大部分重复产生的正面次数在 45 到 55 的范围内。

虽然在理论上“可能”正面次数在 0 到 100 之间,但模拟显示“可能的”值范围要小得多。

这是抛硬币波动性更一般现象的一个实例,我们将在后续课程中看到。

示例:大富翁中的移动

大富翁游戏中每次移动由两次掷骰子的总点数决定。如果你玩大富翁,你应当预期两次掷骰子得到什么结果?

我们可以通过模拟两次掷骰子的和来探索这个问题。我们将运行 10,000 次模拟。注意在本段中我们已经完成了模拟过程的第 1 步和第 3 步。

第 2 步是我们编写一个函数来模拟一对骰子的总点数。首先,让我们规划我们的代码。我们将创建一个包含数字 1 到 6 的数组,从数组中随机放回地抽取两次,然后将抽取的两个数字相加。

[In ]:
die = np.arange(1, 7)
sum(np.random.choice(die, 2))
7

我们可以使用数组 die 和上面的表达式来定义一个模拟大富翁游戏中一次移动的函数。

[In ]:
def one_simulated_move():
    return sum(np.random.choice(die, 2))

现在我们可以通过从一个空的收集数组开始,并用每个新的模拟移动来扩充它,从而创建一个包含 10,000 次模拟的大富翁移动的数组。

[In ]:
num_repetitions = 10000

moves = make_array()
for i in np.arange(num_repetitions):
    new_move = one_simulated_move()
    moves = np.append(moves, new_move)

以下是结果的直方图。

[In ]:
results = Table().with_columns(
    'Repetition', np.arange(1, num_repetitions + 1),
    'Sum of Two Rolls', moves
)

results.hist('Sum of Two Rolls', bins = np.arange(1.5, 12.6, 1))
Histogram with 'Sum of Two Rolls' on the x-axis and 'Percent per unit' on the y-axis. The bars have a width of 1 and range centered on 2 to 12. The graph is bell shaped with the peak centered at 7.

7 是最常见的值,频数在其两侧对称下降。