置信区间
我们已经开发了一种使用随机抽样和自助法来估计参数的方法。我们的方法产生一个估计值的区间，以考虑随机样本中的随机变异。通过提供一个估计区间而不是仅仅一个估计值，我们为自己留出了一些余地。
在前面的例子中，我们看到我们的估计过程在大约95%的情况下产生了一个好的区间，“好的”区间是指包含参数的区间。我们说我们“有95%的信心”这个过程会产生一个好的区间。我们的估计区间被称为参数的“95%置信区间”（confidence interval），95%被称为区间的“置信水平”（confidence level）。
该方法被称为“自助法百分位数方法”（bootstrap percentile method），因为该区间是通过取出自助法估计值的两个百分位数形成的。
前一个例子中的情况有点不寻常。因为我们碰巧知道参数的值，所以我们能够检查一个区间是好的还是无效的，这反过来帮助我们看到了我们的估计过程在我们每次使用中大约有95次捕获了参数。
但通常，数据科学家不知道参数的值。这恰恰是他们首先想要估计参数的原因。在这种情况下，他们使用像我们开发的方法这样的技术为未知参数提供一个估计区间。由于统计理论和我们所见过的这类演示，数据科学家可以相信他们生成区间的过程在已知百分比的情况下会产生一个好的区间。
from datascience import *
%matplotlib inline
path_data = '../../../assets/data/'
import matplotlib.pyplot as plots
plots.style.use('fivethirtyeight')
import numpy as np估计总体中位数
现在我们将使用自助法来估计一个未知的总体中位数。你之前遇到过这个数据集。它来自一个大型医院系统的新生儿样本。尽管抽样是分多个阶段进行的，但我们将把它视为简单随机样本。Deborah Nolan和Terry Speed合著的Stat Labs提供了关于更大数据集（本数据集从中抽取）的详细信息。
表格 births 包含母婴对的以下变量：婴儿出生体重（盎司）、妊娠天数（母亲怀孕的天数）、母亲年龄（周岁）、母亲身高（英寸）、孕期体重（磅），以及母亲在怀孕期间是否吸烟。
births = Table.read_table(path_data + 'baby.csv')births.show(3)<IPython.core.display.HTML object>出生体重是新生儿健康的一个重要因素。较小的婴儿在出生头几天往往比较大的新生儿需要更多的医疗护理。因此，在婴儿出生前对出生体重有一个估计是很有帮助的。一种方法是考察出生体重与妊娠天数之间的关系。
这种关系的一个简单度量是出生体重与妊娠天数的比率。表格 ratios 包含了 baby 的前两列，以及一列比率值。该列的第一个条目计算如下：
$$ \frac{120~\mbox{盎司}}{284~\mbox{天}} ~\approx ~ 0.4225~ \mbox{盎司/天} $$
ratios = births.select('Birth Weight', 'Gestational Days').with_columns(
'Ratio BW:GD', births.column('Birth Weight')/births.column('Gestational Days')
)ratiosBirth Weight | Gestational Days | Ratio BW:GD
120 | 284 | 0.422535
113 | 282 | 0.400709
128 | 279 | 0.458781
108 | 282 | 0.382979
136 | 286 | 0.475524
138 | 244 | 0.565574
132 | 245 | 0.538776
120 | 289 | 0.415225
143 | 299 | 0.478261
140 | 351 | 0.39886
... (1164 rows omitted)以下是这些比率的直方图。
ratios.select('Ratio BW:GD').hist()Histogram with 'Ratio BW:GD' on the x-axis and 'Percent per unit' on the y-axis. The histogram is roughly bell shaped from about 0.3 to 0.6, with a right tail extending out to about 0.8.乍看之下，直方图看起来相当对称，密度在每天0.4盎司到每天0.45盎司区间达到最大值。但仔细观察会发现，有些比率相比之下相当大。比率的最大值刚超过每天0.78盎司，几乎是典型值的两倍。
ratios.sort('Ratio BW:GD', descending=True).take(0)Birth Weight | Gestational Days | Ratio BW:GD
116 | 148 | 0.783784中位数体现了典型比率的概念，因为它不受非常大或非常小的比率的影响。样本中的比率中位数约为每天0.429盎司。
percentile(50, ratios.column(2))0.42907801418439717但总体中的中位数是多少？我们不知道，所以我们将估计它。
我们的方法将与上一节完全相同。我们将对样本进行5,000次自助法，得到5,000个中位数估计值。我们的95%置信区间将是所有这些估计值的“中间95%”。
构建自助法置信区间
我们首先定义一个函数 one_bootstrap_median。它将对样本进行自助法，并返回自助法样本中的中位数比率。
def one_bootstrap_median():
resample = ratios.sample()
return percentile(50, resample.column('Ratio BW:GD'))运行下面的单元格，观察自助法比率的变化。请记住，其中每一个都是对总体中未知比率的估计。
one_bootstrap_median()0.43010752688172044现在我们可以使用 for 循环生成5000个自助法中位数。
# Generate medians from 5000 bootstrap samples
num_repetitions = 5000
bstrap_medians = make_array()
for i in np.arange(num_repetitions):
bstrap_medians = np.append(bstrap_medians, one_bootstrap_median())# Get the endpoints of the 95% confidence interval
left = percentile(2.5, bstrap_medians)
right = percentile(97.5, bstrap_medians)
make_array(left, right)array([0.42545455, 0.43272727])95%置信区间从每天约0.425盎司到每天约0.433盎司。我们估计总体中“出生体重与妊娠天数”比率的中位数在每天0.425盎司到每天0.433盎司之间的某个位置。
基于原始样本的估计值0.429恰好位于区间两端的一半位置，尽管一般情况不一定如此。
为了可视化我们的结果，让我们绘制自助法中位数的经验直方图，并将置信区间放在横轴上。
resampled_medians = Table().with_columns(
'Bootstrap Sample Median', bstrap_medians
)
resampled_medians.hist(bins=15)
plots.plot([left, right], [0, 0], color='yellow', lw=8);Histogram with 'Bootstrap Sample Median' on the x-axis and 'Percent per unit' on the y-axis. The histogram's three center bars are the tallest and extend from 0.4275 to 0.4325. Other bars are present from about 0.425 to 0.435这个直方图和区间与我们上一节绘制的相似，但有一个很大的不同——没有绿点显示参数的位置。我们不知道那个点应该在哪里，甚至不知道它是否在区间内。
我们只有一个估计区间。这是一个估计值的95%置信区间，因为生成它的过程大约在95%的情况下会产生一个好的区间。这当然比随机猜测比率要好！
请记住，这个区间是一个近似的95%置信区间。其计算涉及许多近似。近似效果不错，但并不精确。
估计总体平均值
我们对中位数所做的操作同样适用于平均值。假设我们想估计总体中母亲的平均年龄。一个自然的估计值是样本中母亲的平均年龄。以下是她们年龄的分布以及她们的平均年龄，约为27.2岁。
births.select('Maternal Age').hist()Histogram with 'Maternal Age' on the x-axis and 'Percent per unit' on the y-axis. The bars increase in height from 15 to about 25, then decreases through 45.np.average(births.column('Maternal Age'))27.228279386712096总体中母亲的平均年龄是多少？我们不知道这个参数的值。
让我们用自助法来估计这个未知参数。为此，我们将修改 bootstrap_median 的代码，转而定义函数 bootstrap_mean。代码是相同的，只是统计量是均值（即平均值）而非中位数，并收集在名为 bstrap_means 的数组中，而不是 bstrap_medians。
def one_bootstrap_mean():
resample = births.sample()
return np.average(resample.column('Maternal Age'))# Generate means from 5000 bootstrap samples
num_repetitions = 5000
bstrap_means = make_array()
for i in np.arange(num_repetitions):
bstrap_means = np.append(bstrap_means, one_bootstrap_mean())# Get the endpoints of the 95% confidence interval
left = percentile(2.5, bstrap_means)
right = percentile(97.5, bstrap_means)
make_array(left, right)array([26.90630324, 27.55962521])95%置信区间从约26.9岁到约27.6岁。也就是说，我们估计总体中母亲的平均年龄在26.9岁到27.6岁之间的某个位置。
注意两端与原始样本中约27.2岁的平均值有多么接近。样本量非常大——1,174名母亲——因此样本平均值变化不大。我们将在下一章进一步探讨这一观察。
下面是5,000个自助法平均年龄的经验直方图，以及总体平均年龄的95%置信区间。
resampled_means = Table().with_columns(
'Bootstrap Sample Mean', bstrap_means
)
resampled_means.hist(bins=15)
plots.plot([left, right], [0, 0], color='yellow', lw=8);Histogram with 'Bootstrap Sample Mean' on the x-axis and 'Percent per unit' on the y-axis. The histogram is roughly bell shaped from about 26.7 to 27.7 There is a yellow bar at the bottom of the graph from about 26.9 to 27.5, encompassing much of the data.再次，原始样本的平均值（27.23岁）接近区间的中心。这并不奇怪，因为每个自助法样本都是从同一个原始样本中抽取的。自助法样本的平均值大致对称地分布在它们所来自的样本平均值的两侧。
还要注意，重抽样均值的经验直方图大致呈对称的钟形，尽管样本年龄的直方图完全不对称：
births.select('Maternal Age').hist()A previously shown histogram is rendered again: Histogram with 'Maternal Age' on the x-axis and 'Percent per unit' on the y-axis. The bars increase in height from 15 to about 25, then decreases through 45.这是概率与统计学中中心极限定理的一个结果。在后面的章节中，我们将看到该定理的内容。
80%置信区间
你可以使用自助法样本均值构建任何置信水平的区间。例如，要为总体平均年龄构建一个80%置信区间，你需要取重抽样均值的“中间80%”。也就是说，你需要分布的两个尾部各占10%，因此端点将是重抽样均值的第10和第90百分位数。
left_80 = percentile(10, bstrap_means)
right_80 = percentile(90, bstrap_means)
make_array(left_80, right_80)array([27.01277683, 27.44293015])resampled_means.hist(bins=15)
plots.plot([left_80, right_80], [0, 0], color='yellow', lw=8);A previously shown histogram is rendered again with a different interval highlighted. The 'Bootstrap Sample Mean' is on the x-axis and is roughly bell shaped extending between 26.7 and 27.7. The highlighted interval is smaller, from about 27.0 to 27.4.这个80%置信区间比95%置信区间短得多。它只从约27.0岁到约27.4岁。虽然这是一组紧凑的估计，但你知道这个过程只有大约80%的时间会产生一个好的区间。
前面的过程产生了更宽的区间，但我们对其生成过程有更大的信心。
要获得高置信水平下的窄置信区间，你必须从更大的样本开始。我们将在下一章中了解原因。
估计总体比例
在样本中，39%的母亲在怀孕期间吸烟。
births.where('Maternal Smoker', are.equal_to(True)).num_rows / births.num_rows0.3909710391822828请记住，比例是0和1的平均值。因此，吸烟母亲的比例也可以使用数组操作如下计算。
smoking = births.column('Maternal Smoker')
np.count_nonzero(smoking) / len(smoking)0.3909710391822828总体中百分之多少的母亲在怀孕期间吸烟？这是一个未知参数，我们可以通过自助法置信区间来估计。步骤与我们估计总体中位数和平均值时的步骤类似。
在一个现在已经熟悉的过程中，我们首先定义一个函数 one_bootstrap_proportion，该函数对样本进行自助法并返回自助法样本中吸烟者的比例。然后，我们使用 for 循环多次调用该函数，并获得自助法比例的第2.5百分位数和第97.5百分位数。
def one_bootstrap_proportion():
resample = births.sample()
smoking = resample.column('Maternal Smoker')
return np.count_nonzero(smoking) / len(smoking)# Generate proportions from 5000 bootstrap samples
bstrap_proportions = make_array()
num_repetitions = 5000
for i in np.arange(num_repetitions):
bstrap_proportions = np.append(bstrap_proportions, one_bootstrap_proportion())# Get the endpoints of the 95% confidence interval
left = percentile(2.5, bstrap_proportions)
right = percentile(97.5, bstrap_proportions)
make_array(left, right)array([0.36286201, 0.41908007])置信区间从约36%到约42%。
resampled_proportions = Table().with_columns(
'Bootstrap Sample Proportion', bstrap_proportions
)
resampled_proportions.hist(bins=15)
plots.plot([left, right], [0, 0], color='yellow', lw=8);Histogram with 'Bootstrap Sample Proportion' on the x-axis and 'Percent per unit' on the y-axis. The data is centered around 0.39 and ranges from about 0.34 to 0.43. A yellow bar extends along the bottom of the histogram from about 0.36 to 0.42.使用自助法百分位数方法的注意事项
自助法是一种优雅而强大的方法。在使用之前，记住以下几点很重要。
-
从一个大的随机样本开始。如果不这样做，该方法可能不起作用。它的成功基于大随机样本（以及因此从样本中得到的重样本）与总体相似。平均律表明，只要随机样本足够大，这一点就很可能成立。
-
为了近似一个统计量的概率分布，最好尽可能多地重复重抽样过程。数千次重复将产生对样本中位数分布的良好近似，特别是当总体分布具有单峰且相当对称时。我们在示例中使用了5,000次重复，但通常建议使用10,000次。
-
自助法百分位数方法在基于大随机样本估计总体中位数或均值时效果良好。然而，与所有估计方法一样，它也有局限性。例如，在以下情况下，它预计不会表现良好：
- 目标是估计总体中的最小值或最大值，或非常低或非常高的百分位数，或受总体中稀有元素影响较大的参数。
- 统计量的概率分布不是大致钟形的。
- 原始样本非常小，比如少于10或15个。