百分位数（Percentile）

关于百分位数的定义，不同的课程、教材甚至不同软件都会存在微妙的区别（包括笔者上的《数理统计》课程中也给出了两种不同的定义），下面给出本课程对百分位数的定义：

对于一个集合，其 $p\%$ 分位数（ $0\leq p\leq 100$ ）指的是集合中满足至少有 $p\%$ 的样本不大于它的最小值。具体而言，若集合中有 $n$ 个值（可以重复），将其从小到大排序，最终取第 $k=\left\lceil n\times\dfrac{p}{100}\right\rceil$ 个数。（ $\lceil m \rceil$ 表示上取整）
python中可以使用numpy.percentile()计算百分位数（可参考官方文档），示例：
```
import numpy as np
a = np.array([25, 12, 15, 14, 19, 23, 25, 29, 33, 35])
print(np.percentile(a, 95,method = "inverted_cdf")) # 35
```
数理统计Cheat Sheet中的两种定义在numpy的method中分别对应average_inverted_cdf和linear。

Bootstrap（自助法）

通常情况下，一组样本只能得到一个估计量。为了衡量估计量的波动，需要在总体中多次抽样，然而现实情况中抽样往往需要耗费大量资源。为解决这一问题，统计学家提出使用Bootstrap法（即自助法），直接从已有样本中进行取样。
具体而言，我们可以将已有样本看作一个新的总体，对其进行有放回抽样，最终得到数量与已有样本相同的新样本。新样本与已有样本的差异源于每个个体都可能被重复抽样。
那么为什么Bootstrap法可以很好地衡量估计量的波动（方差）呢？原因很显然——Bootstrap得到的样本也可以看作总体得到的样本，因而其经验分布同样可以很好地反映总体的分布性质。

以对样本中位数的Bootstrap为例：

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
population = np.random.normal(loc=0, scale=1, size=100000) # 总体（标准正态分布）

# 抽取一个原始样本
n = 50
data = np.random.choice(population, size=n, replace=False)
sample_median = np.median(data)
print("样本中位数：", sample_median)

# Bootstrap
B = 5000
bootstrap_medians = []
for _ in range(B):
    sample = np.random.choice(data, size=n, replace=True)
    bootstrap_medians.append(np.median(sample))
bootstrap_medians = np.array(bootstrap_medians)

plt.figure(figsize=(8,5), dpi=300)
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.hist(
    bootstrap_medians,
    bins=40,
    density=True,
    alpha=0.5,
    label="Bootstrap分布"
)
plt.axvline(
    sample_median,
    color="blue",
    linestyle="--",
    linewidth=2,
    label="样本中位数"
)
plt.xlabel("样本中位数")
plt.ylabel("概率密度")
plt.title("Bootstrap分布")
plt.legend()
plt.show()

结果分布如图：

可以看到Boostrap样本得到的中位数的确分布在总体中位数

0

附近。

那么，Boostrap样本中位数分布一定能覆盖总体中位数吗？我们只考虑这个分布的中间95%数据，并重复Bootstrap 100次（得到100个Bootstrap分布），代码如下：

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
# 总体参数
true_median = 0
n = 50 
B = 5000 
R = 100 
lower_bounds = []
upper_bounds = []
covered = []

for _ in range(R):
    # 从总体重新抽一个样本
    sample = np.random.normal(loc=0, scale=1, size=n)
    # Bootstrap
    bootstrap_medians = []
    for _ in range(B):
        resample = np.random.choice(sample, size=n, replace=True)
        bootstrap_medians.append(np.median(resample))
    bootstrap_medians = np.array(bootstrap_medians)

    # Percentile Bootstrap 区间
    lower = np.percentile(bootstrap_medians,2.5)
    upper = np.percentile(bootstrap_medians,97.5)
    lower_bounds.append(lower)
    upper_bounds.append(upper)
    covered.append(lower <= true_median <= upper)

lower_bounds = np.array(lower_bounds)
upper_bounds = np.array(upper_bounds)
covered = np.array(covered)
coverage = covered.mean()
print(f"覆盖率 = {coverage:.3f}")

结果表明，其覆盖总体中位数的概率也在95%左右，如图所示： interval

换句话说，总体中位数的一个95%置信区间可以近似用Bootstrap分布的第2.5%和第97.5%分位数表示。【这一方法也称为Bootstrap百分位数法】

当然，除了中位数，我们也可以对其他统计量（如均值，比率等）进行Bootstrap估计，其对应 $1-\alpha$ 置信区间都可以用Bootstrap分布的 $\dfrac{\alpha}{2}$ 百分位数和 $1-\dfrac{\alpha}{2}$ 百分位数刻画。

总结一下Bootstrap的通用代码框架：

import numpy as np
def bootstrapper(sample, statistic, num_repetitions):
    """
    返回在 num_repetitions 个来自原始样本的 bootstrap 样本上计算的统计量。
    """
    bstrap_stats = []
    for i in np.arange(num_repetitions):
        # 步骤 1：对样本进行重抽样（有放回）
        indices = np.random.choice(len(sample), size=len(sample), replace=True)
        # 步骤 2：计算重抽样样本上的统计量
        bootstrap_stat = statistic(bootstrap_sample)
        # 累积统计量
        bstrap_stats = np.append(bstrap_stats, bootstrap_stat)
    return bstrap_stats

置信区间（Confidence Interval）

在国内的数理统计课程中，同样涉及到了置信区间的概念（可参考数理统计Cheat Sheet），但其主要使用的是枢轴量法进行构造。
而在这里，我们主要使用的是重采样模拟构造置信区间，这往往需要样本量足够大。但Bootstrap法构造置信区间的优势在于其不需要知道总体的分布，属于一种非参数估计。
Bootstrap采样次数一般和已有样本数相同，但采样次数越多，模拟产生的误差也会越低（不过边际效益会递减）。另外，总体分布本身在一定程度上也会影响Bootstrap区间估计的准确度。
对于一些易受异常值影响的统计量（如最小/最大次序统计量等），也不适合用Bootstrap进行估计。

或许之后非参数统计中也会涉及这些内容……

DATA8 CHAPTER 9

DATA8 CHAPTER 9

百分位数（Percentile）

Bootstrap（自助法）

置信区间（Confidence Interval）