自助法

一位数据科学家正在使用随机样本中的数据来估计一个未知参数。她使用样本来计算一个统计量的值，并将其作为她的估计值。

一旦她计算出统计量的观测值，她本可以将其作为估计值呈现，然后继续她的工作。但她是一位数据科学家。她知道她的随机样本只是众多可能随机样本中的一个，因此她的估计值也只是众多合理估计值中的一个。

这些估计值可能有多大差异？要回答这个问题，看起来她需要从总体中抽取另一个样本，并基于新样本计算一个新的估计值。但她没有资源回到总体中再抽取一个样本。

看起来这位数据科学家陷入了困境。

幸运的是，一个名为“自助法”（bootstrap）的巧妙想法可以帮助她。由于从总体中生成新样本不可行，自助法通过一种称为“重抽样”（resampling）的方法生成新的随机样本：新样本是从原始样本中随机抽取的。

[In ]:

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

在本节中，我们将了解自助法的工作原理及其原因。在本章的其余部分，我们将使用自助法进行推断。

旧金山市员工薪酬

SF OpenData 是旧金山市和县公开部分数据的网站。其中一个数据集包含该市员工的薪酬数据。这些员工包括市立医院的医疗专业人员、警察、消防员、交通工人、民选官员以及该市的所有其他员工。

2019日历年的薪酬数据在表格 sf2019 中。

[In ]:

sf2019 = Table.read_table(path_data + 'san_francisco_2019.csv')

[In ]:

sf2019.show(3)

<IPython.core.display.HTML object>

超过44,500名员工每人对应一行。有众多列包含关于城市部门归属以及员工薪酬包不同部分的详细信息。以下是2019年旧金山市长伦敦·布里德（London Breed）对应的行。

[In ]:

sf2019.where('Job', 'Mayor')

Organization Group               | Department | Job Family                    | Job   | Salary | Overtime | Benefits | Total Compensation
General Administration & Finance | Mayor      | Administrative & Mgmt (Unrep) | Mayor | 342974 | 0        | 98012    | 440987

我们将研究最后一列 Total Compensation（总薪酬）。这是员工的工资加上市府为其退休和福利计划缴纳的部分。

一个日历年内的薪酬方案有时可能难以理解，因为它们取决于入职日期、员工是否在市内换工作等因素。例如，Total Compensation 列中的最低值看起来有点奇怪。

[In ]:

sf2019.sort('Total Compensation')

Organization Group                      | Department         | Job Family                    | Job                            | Salary | Overtime | Benefits | Total Compensation
Public Protection                       | Adult Probation    | Probation & Parole            | Deputy Probation Officer       | 0      | 0        | 0        | 0
Public Protection                       | Fire Department    | Clerical, Secretarial & Steno | Senior Clerk Typist            | 0      | 0        | 0        | 0
Public Protection                       | Juvenile Court     | Correction & Detention        | Counselor, Juvenile Hall PERS  | 0      | 0        | 0        | 0
Public Protection                       | Police             | Clerical, Secretarial & Steno | Clerk Typist                   | 0      | 0        | 0        | 0
Public Protection                       | Sheriff            | Correction & Detention        | Deputy Sheriff                 | 0      | 0        | 0        | 0
Public Works, Transportation & Commerce | Airport Commission | Sub-Professional Engineering  | StdntDsgn Train2/Arch/Eng/Plng | 0      | 0        | 0        | 0
Public Works, Transportation & Commerce | Airport Commission | Clerical, Secretarial & Steno | Executive Secretary 1          | 0      | 0        | 0        | 0
Public Works, Transportation & Commerce | Airport Commission | Payroll, Billing & Accounting | Senior Account Clerk           | 0      | 0        | 0        | 0
Public Works, Transportation & Commerce | Airport Commission | Housekeeping & Laundry        | Custodian                      | 0      | 0        | 0        | 0
Public Works, Transportation & Commerce | Airport Commission | Housekeeping & Laundry        | Custodian                      | 0      | 0        | 0        | 0
... (44515 rows omitted)

为了解释的清晰性，我们将重点关注那些全年大致相当于半职或以上工作的人员。按每小时约15美元的最低工资、每周20小时、52周计算，即超过15,000美元的工资。

[In ]:

sf2019 = sf2019.where('Salary', are.above(15000))

[In ]:

sf2019.num_rows

总体与参数

让这个刚超过37,000行的表格作为我们的总体。以下是该表中员工总薪酬的直方图。

[In ]:

sf_bins = np.arange(0, 726000, 25000)
sf2019.select('Total Compensation').hist(bins=sf_bins)

Histogram with 'Total Compensation' on the x-axis and 'Percent per unit' on the y-axis. The bars are roughly bell shaped centered above 100000 with a long tail out to 700000.

虽然大多数值低于300,000美元，但也有少数值高出不少。例如，首席投资官的总薪酬超过700,000美元。这就是为什么横轴在可见柱状图右侧延伸得很远。

[In ]:

sf2019.sort('Total Compensation', descending=True).show(2)

<IPython.core.display.HTML object>

假设我们关注的参数是总薪酬的中位数。

由于我们幸运地拥有来自总体的所有数据，我们可以直接计算该参数：

[In ]:

pop_median = percentile(50, sf2019.column('Total Compensation'))
pop_median

135747.0

所有员工的总薪酬中位数为135,747美元。

从实践角度来看，我们没有理由抽取样本来估计这个参数，因为我们直接知道它的值。但在本节中，我们将假装不知道这个值，并看看基于随机样本我们能多好地估计它。

在后面的章节中，我们将回归现实，在参数未知的情况下工作。现在，我们是全知的。

随机样本与估计值

让我们不放回地随机抽取500名员工样本，并将样本员工的总薪酬中位数作为我们对参数的估计值。

[In ]:

our_sample = sf2019.sample(500, with_replacement=False)
our_sample.select('Total Compensation').hist(bins=sf_bins)

Histogram with 'Total Compensation' on the x-axis and 'Percent per unit' on the y-axis. The shape of the histogram is about the same as in the previous histogram, being roughly bell shaped and centered above 100000 with a long right tail out to 700000, but the exact heights of the bars are different.

[In ]:

est_median = percentile(50, our_sample.column('Total Compensation'))
est_median

136835.0

样本量很大。根据平均律，样本的分布与总体的分布相似。因此，样本中位数与总体中位数相当接近，当然并非完全相同。

所以现在我们有了一个参数的估计值。但如果样本结果不同，估计值也会有不同的值。我们希望能量化估计值在不同样本之间可能变化的量。这种变异性的度量将帮助我们衡量我们估计参数的准确程度。

要了解如果样本结果不同，估计值会有多大差异，我们可以直接从总体中再抽取一个样本。但那将是作弊。我们试图模拟现实生活，在现实中我们手头不会有所有的总体数据。

我们必须以某种方式获得另一个随机样本，而无需再次从总体中抽样。

自助法：从样本中重抽样

我们拥有的是一个来自总体的大随机样本。正如我们所知，大随机样本很可能与它所来自的总体相似。这一观察使数据科学家能够“用自己的鞋带把自己提起来”（自助）：抽样过程可以通过“从样本中抽样”来复制。

以下是生成另一个与总体相似的随机样本的“自助法”步骤：

将原始样本视为总体。
从样本中有放回地随机抽取，次数与原始样本量相同。

以与原始样本量相同的次数进行重抽样很重要。原因是估计值的变异性取决于样本量。由于我们的原始样本由500名员工组成，我们的样本中位数基于500个值。要了解样本可能有多大不同，我们必须将其与其他容量为500的样本的中位数进行比较。

如果我们从容量为500的样本中无放回地随机抽取500次，我们只会得到相同的样本。通过有放回地抽取，我们创造了新样本与原始样本不同的可能性，因为某些员工可能被抽中多次，而其他员工则根本未被抽中。

为什么自助法有效

为什么这是一个好主意？根据平均律，原始样本的分布很可能与总体相似，而所有“重样本”的分布很可能与原始样本相似。因此，所有重样本的分布也很可能与总体相似。

自助法示意图，注释为“所有这些看起来都很相似，极有可能。左侧是显示总体的图形，箭头指向一个样本。样本看起来与总体相似，但不完全相同。从样本出发，展示了多个重样本，它们都有不同的分布，但也与样本和总体相似。

[In ]:

from IPython.display import Image
Image("../../../images/bootstrap_pic.png")

<IPython.core.display.Image object>

一个重抽样的中位数

回顾一下，sample 方法默认从表格中有放回地抽取行，当不指定样本量时，默认样本量等于表格的行数。这对自助法来说非常完美！以下是从原始样本中抽取的一个新样本，以及对应的样本中位数。

[In ]:

resample_1 = our_sample.sample()

[In ]:

resample_1.select('Total Compensation').hist(bins=sf_bins)

Histogram with 'Total Compensation' on the x-axis and 'Percent per unit' on the y-axis. The histogram is mostly bell shaped and centered above 100000 with a tail extending out to 700000. There are some notable differences from the previous two histograms. For example, the visible bars furtherest to the right are taller than in either previous histogram.

[In ]:

resampled_median_1 = percentile(50, resample_1.column('Total Compensation'))
resampled_median_1

141793.0

这个值是对总体中位数的一个估计。

通过反复重抽样，我们可以得到许多这样的估计值，从而得到估计值的经验分布。

[In ]:

resample_2 = our_sample.sample()
resampled_median_2 = percentile(50, resample_2.column('Total Compensation'))
resampled_median_2

135880.0

让我们整理这段代码，定义一个函数 one_bootstrap_median，该函数基于对我们称之为 our_sample 的原始随机样本进行自助法抽样，返回总薪酬的一个自助法中位数。

[In ]:

def one_bootstrap_median():
    resampled_table = our_sample.sample()
    bootstrapped_median = percentile(50, resampled_table.column('Total Compensation'))
    return bootstrapped_median

多次运行下面的单元格，观察自助法中位数如何变化。请记住，其中每一个都是对总体中位数的估计。

[In ]:

one_bootstrap_median()

132175.0

样本中位数的自助法经验分布

现在我们可以像往常一样通过运行 for 循环来多次重复自助法过程。在每次迭代中，我们将调用函数 one_bootstrap_median 来基于原始样本 our_sample 生成一个自助法中位数的值。然后将该自助法中位数追加到收集数组 bstrap_medians 中。

由于我们要求重复5000次，代码可能需要运行一段时间。它需要做大量的重抽样！

[In ]:

num_repetitions = 5000
bstrap_medians = make_array()
for i in np.arange(num_repetitions):
    bstrap_medians = np.append (bstrap_medians, one_bootstrap_median())

以下是5000个自助法中位数的经验直方图。绿点是总体参数：即整个总体的中位数，也就是我们试图估计的值。在这个例子中，我们碰巧知道它的值，但我们没有在自助法过程中使用它。

[In ]:

resampled_medians = Table().with_column('Bootstrap Sample Median', bstrap_medians)
median_bins=np.arange(120000, 160000, 2000)
resampled_medians.hist(bins = median_bins)

# Plotting parameters; you can ignore this code
parameter_green = '#32CD32'
plots.ylim(-0.000005, 0.00014)
plots.scatter(pop_median, 0, color=parameter_green, s=40, zorder=2);

Histogram with 'Bootstrap Sample Median' on the x-axis and 'Percent per unit' on the y-axis. Ther histogram has data between 120000 and 150000, roughly bell shaped although the highest bars aren't in the exact center of the distribution. There is a green point in the middle of the distribution.

重要的是要记住，绿点是固定的：即135,747美元，总体中位数。经验直方图是随机抽取的结果，并且相对于绿点会随机定位。

还要记住，所有这些计算的目的都是为了估计总体中位数（即绿点）。我们的估计值就是你上面直方图中看到的所有随机生成的样本中位数。我们希望这组估计值包含该参数。如果不包含，那么估计值就是偏离的。

估计值是否捕获了参数？

重抽样中位数的经验直方图多久能稳稳地覆盖绿点，而不是仅仅用尾部扫过绿点或根本不覆盖它？要回答这个问题，我们必须定义“稳稳覆盖”。我们将其理解为“重抽样中位数的中间95%包含绿点”。

以下是重抽样中位数“中间95%”区间的两端：

[In ]:

left = percentile(2.5, bstrap_medians)
left

129524.0

[In ]:

right = percentile(97.5, bstrap_medians)
right

143446.0

总体中位数135,747美元介于这两个数之间。该区间和总体中位数显示在下面的直方图中。

[In ]:

resampled_medians.hist(bins = median_bins)

# Plotting parameters; you can ignore this code
plots.ylim(-0.000005, 0.00014)
plots.plot([left, right], [0, 0], color='yellow', lw=3, zorder=1)
plots.scatter(pop_median, 0, color=parameter_green, s=40, zorder=2);

The previous histogram is shown again, including the green dot in the center of the distribution. There is now also a yellow bar that extends within the bulk of the data in the histogram, from about 130000 to approx 144000.

在我们的示例中，估计值的“中间95%”区间捕获了参数。但这是侥幸吗？

要了解区间包含参数的频率，我们必须反复运行整个过程。具体来说，我们将重复以下过程100次：

从总体中抽取一个容量为500的原始随机样本。
执行5000次自助法过程的重复，并生成重抽样中位数的“中间95%”区间。

我们将得到100个区间，并统计其中有多少个包含总体中位数。

剧透警告： 自助法的统计理论表明这个数字应该在95左右。它可能在90出头或90多，但我们预计不会偏离95太远。

我们首先编写一个函数 bootstrap_median，它接受两个参数：包含原始随机样本的表格名称，以及要抽取的自助法样本数量。它返回一个自助法中位数的数组，每个自助法样本对应一个。

[In ]:

def bootstrap_median(original_sample, num_repetitions):
    medians = make_array()
    for i in np.arange(num_repetitions):
        new_bstrap_sample = original_sample.sample()
        new_bstrap_median = percentile(50, new_bstrap_sample.column('Total Compensation'))
        medians = np.append(medians, new_bstrap_median)
    return medians

现在我们将编写一个 for 循环，调用此函数100次，并每次收集自助法中位数的“中间95%”。

下面的单元格将需要几分钟才能运行完成，因为它需要执行100次重复，每次从表格中随机抽样500次并生成5000个自助法样本。

[In ]:

# THE BIG SIMULATION: This one takes several minutes.

# Generate 100 intervals and put the endpoints in the table intervals

left_ends = make_array()
right_ends = make_array()

for i in np.arange(100):
    original_sample = sf2019.sample(500, with_replacement=False)
    medians = bootstrap_median(original_sample, 5000)
    left_ends = np.append(left_ends, percentile(2.5, medians))
    right_ends = np.append(right_ends, percentile(97.5, medians))

intervals = Table().with_columns(
    'Left', left_ends,
    'Right', right_ends
)

对于整个过程的100次重复中的每一次，我们得到一个中位数估计值的区间。

[In ]:

intervals

Left   | Right
125093 | 139379
129925 | 140757
133955 | 146369
129335 | 140847
132756 | 145429
130167 | 143200
125935 | 138491
131092 | 142472
128509 | 140462
131270 | 145998
... (90 rows omitted)

好的区间是那些包含我们试图估计的参数的区间。通常参数是未知的，但在本节中我们碰巧知道参数是什么。

[In ]:

pop_median

135747.0

100个区间中有多少个包含总体中位数？即左端低于总体中位数且右端高于总体中位数的区间数量。

[In ]:

intervals.where(
    'Left', are.below(pop_median)).where(
    'Right', are.above(pop_median)).num_rows

构建所有区间需要好几分钟，但如果你有耐心，可以再试一次。很可能，100个区间中约有95个会是好的区间：它们将包含参数。

很难在横轴上向你展示所有的区间，因为它们有大量重叠——毕竟，它们都在试图估计同一个参数。下图通过在纵向上堆叠来在同一坐标轴上显示每个区间。纵轴就是生成该区间的重复次数编号。

绿线是参数所在的位置。由于参数是固定的，它的位置也是固定的。

好的区间覆盖参数。通常大约有95个这样的区间。

如果一个区间没有覆盖参数，它就是无效的。无效区间是那些你可以在绿线周围看到“空隙”的区间。它们非常少——通常100个中大约有5个——但它们确实会发生。

任何基于抽样的方法都有可能出现偏差。基于随机抽样的方法的美妙之处在于，我们可以量化它们可能出错的频率。

[In ]:

replication_number = np.ndarray.astype(np.arange(1, 101), str)
intervals2 = Table(replication_number).with_rows(make_array(left_ends, right_ends))

plots.figure(figsize=(8,8))
for i in np.arange(100):
    ends = intervals2.column(i)
    plots.plot(ends, make_array(i+1, i+1), color='gold')
plots.scatter(pop_median, 0, color=parameter_green, s=40, zorder=2)
plots.plot(make_array(pop_median, pop_median), make_array(0, 100), color=parameter_green, lw=2)
plots.xlabel('Median (dollars)')
plots.ylabel('Replication')
plots.title('Population Median and Intervals of Estimates');

Graph titled 'Population Median and Intervals of Estimates' with 'Median (dollars)' on the x-axis and 'Replication' on the y-axis. 100 replications are shown, each showing a gold interval inside of the range from about 125000 to 155000. The intervals are all different within the range. A green line extends from a green dot just about 135000 up through the intervals. The green line intersects with most, but not all of the intervals.

总结一下模拟的结果：假设你通过以下过程估计总体中位数：

从总体中抽取一个大随机样本。
对你的随机样本进行自助法，并从新的随机样本中获取一个估计值。
将上述自助法步骤重复数千次，得到数千个估计值。
取出所有估计值的“中间95%”区间。

这就得到了一个估计区间。如果其他99个人重复整个过程，每次都从一个新的随机样本开始，那么你将得到100个这样的区间。这100个区间中约有95个会包含总体参数。

换句话说，这个估计过程大约在95%的情况下捕获了参数。

你可以用不同的值替换95%，只要不是100%即可。假设你将95%替换为80%，并将样本量保持在500。那么你的估计区间将比我们这里模拟的更短，因为“中间80%”的范围比“中间95%”更小。如果你不断重复这个过程，只有大约80%的区间会包含参数。