stats 设计

设计目标

基准测试计时是偏斜的、重尾的，并且与采集时刻相关。stats 给出能经受这些性质考验的决策：稳健的位置估计、基于成对测量的比较、显式的实际阈值，以及随机性由种子固定的区间。每个函数都是数组的纯变换，因此决策可以从事件流中保存的原始观测重新计算。

数学背景

顺序统计量与样本分位数

设 $x_1, \dots, x_n$ 为样本，$x_{(1)} \le \dots \le x_{(n)}$ 为其顺序统计量。本包在所有地方都使用同一个分位数定义，即 Hyndman 和 Fan 第 7 型的线性插值：¹1 R. J. Hyndman 和 Y. Fan，“Sample quantiles in statistical packages”，The American Statistician 50(4)，1996。第 7 型是 R 和 NumPy 的默认定义。

$$h = (n-1)\,p, \qquad Q(p) = x_{(\lfloor h\rfloor+1)} + (h-\lfloor h\rfloor)\bigl(x_{(\lfloor h\rfloor+2)} - x_{(\lfloor h\rfloor+1)}\bigr), \qquad 0 \le p \le 1 .$$

由该公式直接得到以下三个性质，下文会用到它们。

1. 范围。 $Q(p)$ 是两个相邻顺序统计量的凸组合，因此 $x_{(1)} \le Q(p) \le x_{(n)}$，且 $Q(0) = x_{(1)}$、$Q(1) = x_{(n)}$。
2. 单调性。 $Q$ 是非递减序列 $x_{(1)}, \dots, x_{(n)}$ 在节点 $p = k/(n-1)$ 处的分段线性插值，因此 $p \le p'$ 蕴含 $Q(p) \le Q(p')$。
3. 仿射同变性。 对于 $y_i = a + b\,x_i$（$b > 0$），排序与该映射可交换，凸组合与仿射映射可交换，因此 $Q_y(p) = a + b\,Q_x(p)$。把计时从 µs 换算为 ns，会使每个分位数、中位数、IQR 和 MAD 按同一因子缩放。

当 $p = 1/2$ 时，公式给出通常的中位数：$n = 2m+1$ 时为 $x_{(m+1)}$，$n = 2m$ 时为 $\tfrac12(x_{(m)} + x_{(m+1)})$。

位置与离散程度

summarize 并列报告两类估计量：

经典估计量	稳健估计量
$\bar x = \frac1n \sum_i x_i$	$\operatorname{med}(x) = Q(1/2)$
$s_n = \sqrt{\frac1n \sum_i (x_i - \bar x)^2}$	$\mathrm{IQR} = Q(3/4) - Q(1/4)$, $\mathrm{MAD} = \operatorname{med}_i \lvert x_i - \operatorname{med}(x)\rvert$

估计量的崩溃点是指：样本中可以被任意远地移动、而估计值不会因此被任意远地移动的最大比例。对 $\bar x$ 和 $s_n$ 它为 $0$（一个值就够了），对四分位数和 IQR 为 $1/4$，对中位数和 MAD 为 $1/2$。一次被调度出去的运行就能移动基准测试的均值；除非一半的运行都受到影响，否则它无法移动中位数。

$s_n$ 除以 $n$：它是把样本视为总体时的标准差，是一个描述性的量。作为 $\sigma^2$ 的估计量它是有偏的，因为

$$\begin{aligned} \sum_i (x_i - \bar x)^2 &= \sum_i (x_i - \mu)^2 - n(\bar x - \mu)^2, \\ \mathbb E\Bigl[\sum_i (x_i - \bar x)^2\Bigr] &= n\sigma^2 - n\cdot\frac{\sigma^2}{n} = (n-1)\,\sigma^2 , \end{aligned}$$

所以 $\mathbb E[s_n^2] = \frac{n-1}{n}\sigma^2$。需要无偏方差时，请把 $s_n^2$ 乘以 $n/(n-1)$。

MAD 以未缩放的形式报告。对于标准差为 $\sigma$ 的正态分布数据，$\operatorname{med}\lvert X - \mu\rvert = z_{3/4}\,\sigma$，其中 $z_{3/4} = \Phi^{-1}(3/4) \approx 0.6745$，因此 $1.4826\cdot\mathrm{MAD}$ 是 $\sigma$ 的估计，且 $\mathrm{IQR} \approx 1.349\,\sigma$。

成对测量

运行器在每个区组中按轮换的顺序对每个实现测量一次（参见 runner 设计）。把来自区组 $i$ 的基线计时 $b_i$ 和候选计时 $c_i$ 建模为

$$b_i = \mu_b + \beta_i + \varepsilon_i, \qquad c_i = \mu_c + \beta_i + \eta_i ,$$

其中 $\beta_i$ 是两者共享的区组效应（机器状态、频率、缓存内容），方差为 $\sigma_\beta^2$，$\varepsilon_i$、$\eta_i$ 是相互独立的噪声，方差分别为 $\sigma_\varepsilon^2$、$\sigma_\eta^2$。成对差值抵消了区组效应：

$$d_i = c_i - b_i = (\mu_c - \mu_b) + (\eta_i - \varepsilon_i), \qquad \operatorname{Var}(d_i) = \sigma_\varepsilon^2 + \sigma_\eta^2 ,$$

而跨区组 $i \ne j$ 取的差值方差为 $\sigma_\varepsilon^2 + \sigma_\eta^2 + 2\sigma_\beta^2$。当区组间的漂移比区组内的噪声更占主导时，配对能消除大部分方差。因此 compare_paired 总结的是差值 $d_i$，而从不分别总结两个样本。

百分位自助法

设 $\hat F_n$ 为差值的经验分布，$\hat\theta = \operatorname{med}(d)$。一个自助法重抽样样本 $d^*$ 从 $\hat F_n$ 中有放回地抽取 $n$ 个值；其中位数为 $\hat\theta^*$。把自助法分布记作 $\hat G(x) = P^*(\hat\theta^* \le x)$。置信度为 $1-\alpha$ 的百分位区间为

$$\bigl[\hat G^{-1}(\alpha/2),\ \hat G^{-1}(1-\alpha/2)\bigr].$$

在下述条件下它是精确的。假设存在一个递增映射 $\varphi$，使得 $W = \varphi(\hat\theta) - \varphi(\theta)$ 具有关于 $0$ 对称的分布 $H$，并且在重抽样下同一个 $H$ 也描述 $\varphi(\hat\theta^*) - \varphi(\hat\theta)$。那么 $\hat G^{-1}(q) = \varphi^{-1}\bigl(\varphi(\hat\theta) + H^{-1}(q)\bigr)$，且

$$\begin{aligned} \theta \le \hat G^{-1}(1-\tfrac\alpha2) &\iff -W \le H^{-1}(1-\tfrac\alpha2) \iff W \ge H^{-1}(\tfrac\alpha2), \\ \theta \ge \hat G^{-1}(\tfrac\alpha2) &\iff -W \ge H^{-1}(\tfrac\alpha2) \iff W \le H^{-1}(1-\tfrac\alpha2), \end{aligned}$$

这里用到了 $H^{-1}(q) = -H^{-1}(1-q)$。覆盖率为 $P\bigl(H^{-1}(\alpha/2) \le W \le H^{-1}(1-\alpha/2)\bigr) = 1-\alpha$。该区间从不需要知道 $\varphi$，这就是它具有变换保持性的原因。²2 B. Efron 和 R. J. Tibshirani，An Introduction to the Bootstrap，Chapman & Hall，1993，§13.3 与 §14。 一般情况下该条件只是近似成立，百分位区间的覆盖误差为 $n^{-1/2}$ 阶。

对于中位数，自助法分布可以精确写出。取 $n = 2m+1$ 个互不相同的值。$\hat\theta^* \le x_{(k)}$ 恰好在 $n$ 次抽取中至少有 $m+1$ 次不超过 $x_{(k)}$ 时成立，而每次抽取满足这一点的概率为 $k/n$，因此

$$\hat G\bigl(x_{(k)}\bigr) = \sum_{j=m+1}^{n} \binom{n}{j}\Bigl(\frac kn\Bigr)^{j}\Bigl(1-\frac kn\Bigr)^{n-j}.$$

$\hat G$ 是定义在数据点上的阶梯函数（当 $n$ 为偶数时，定义在相邻顺序统计量的中点上）。因此区间端点会落在观测到的差值上或其之间，并且在区组较少时以粗糙的步长移动。

两个结果的完整推导以及离散性界见附件：

成对差值中位数的百分位自助法

设计决策

以中位数作为默认中心

问题。 计时分布有一个硬性的下界（工作本身）和一条长长的右尾（中断、缺页、频率变化）。方案。 均值、截尾均值、中位数。选择。 每个决策都使用中位数，均值作为诊断量保留在 SummaryStats 中。理由。 中位数的崩溃点为 $1/2$，具有仿射同变性，而均值与中位数之间的距离提示存在值得在报告中查看的尾部。

成对差值，而非两个独立样本

问题。 机器状态在区组之间漂移。方案。 比较 $\operatorname{med}(c) - \operatorname{med}(b)$；比较成对差值。选择。 compare_paired 计算 $d_i = c_i - b_i$ 并总结 $\operatorname{med}(d)$。理由。 见上文的方差推导：区组效应在 $d_i$ 中抵消，而在非成对差值中保留。代价是一份契约：两个数组的索引 $i$ 必须来自同一个数据集、重复和区组。paired_deltas 截断到较短的数组，compare_paired 把长度不等标记为 Invalid，因此对齐错误会显现出来，而不是被悄悄修补。

相对差值与加速比

问题。 门禁以百分比表述，而人们读的是加速比。选择。 令 $m_b = \operatorname{med}(b)$，$m_d = \operatorname{med}(d)$，

$$r = 100\,\frac{m_d}{m_b}, \qquad s = \frac{m_b}{m_b + m_d} = \frac{1}{1 + m_d/m_b} = \frac{1}{1 + r/100}.$$

两者都是同两个中位数的函数，因此它们永远不会矛盾：在 $r > -100$ 上 $s$ 关于 $r$ 递减，对于阈值 $0 \le t < 100$

$$r \le -t \iff 1 + \frac{r}{100} \le 1 - \frac{t}{100} \iff s \ge \frac{1}{1 - t/100}.$$

阈值 $t = 2$ 恰好在 $s \ge 1.0204$ 时宣布候选为 Faster。当候选是基线的常数平移 $c_i = b_i + \delta$ 时，$m_b + m_d = \operatorname{med}(c)$，$s$ 就是两个中位数之比；一般而言，$m_b + m_d$ 是根据成对数据构建的典型候选时间的稳健估计。保护条件 $m_b = 0 \Rightarrow r = 0$ 和 $m_d = 0 \Rightarrow s = 1$ 使退化输入保持有限；$r = -100$（$m_b + m_d = 0$）仍会给出无穷大的加速比。

实际阈值，而非显著性检验

问题。 只要重复次数足够多，任何差异都会变得统计显著，包括没人会为之采取行动的 0.1 % 变化。方案。 $t$ 检验或 Wilcoxon 检验；等价性检验（TOST）；对点估计使用实际阈值。选择。 决策是三分规则

$$\text{Faster} \iff r \le -t, \qquad \text{Slower} \iff r \ge t, \qquad \text{Equivalent} \iff -t < r < t ,$$

在 Invalid（长度不等）和 Unknown（没有配对）之后检查。理由。 阈值以决策的单位陈述用户真正关心的问题（“它是否至少快 2 %？”），该规则可以由两个中位数复现，并且不会奖励运行更多的重复。对于 $t > 0$，三个区域构成 $\mathbb R$ 的一个划分。对于 $t = 0$，前两个区域在 $r = 0$ 处重叠，并且由于先检验 Faster，恰好持平会被报告为 Faster；请使用正的阈值。不确定性通过区间在决策旁边报告，而不是折叠进决策中。

对中位数差值的带种子百分位自助法

问题。 报告的读者需要看到中位数差值有多稳定，而且重新生成报告时这个数字必须保持不变。选择。 bootstrap_interval 用显式种子抽取 $B$ 个重抽样样本，并返回 $B$ 个重抽样中位数的第 7 型分位数 $Q^*(\alpha/2)$ 和 $Q^*(1-\alpha/2)$。理由。 百分位方法不需要中位数的方差公式（那需要密度估计），它具有变换保持性，而且代价低廉：$O(B\,n\log n)$。

随机流由一个 64 位 xorshift 步骤加一次乘法生成：

$$x \leftarrow x \oplus (x \gg 12), \quad x \leftarrow x \oplus (x \ll 25), \quad x \leftarrow x \oplus (x \gg 27), \quad x \leftarrow 2685821657736338717 \cdot x \bmod 2^{64}.$$

这些是 Marsaglia–Vigna xorshift64* 的常数，但相乘后的值会作为下一个状态反馈回去，因此 xorshift64* 经典的周期结论并不适用。仍然适用的是：每个 xorshift 步骤都是 $\mathrm{GF}(2)^{64}$ 上的可逆线性映射（一个幺幂三角矩阵），而乘以奇数常数在模 $2^{64}$ 下可逆，因此一步就是 64 位字上一个固定 $0$ 的置换。这就是种子 0 会被替换为常数 88172645463393265 的原因。生成器只使用回绕的 64 位整数算术，因此该流在每个目标上都完全相同。

状态 $x$ 通过 $j = x \bmod n$ 映射为索引。记 $2^{64} = qn + r$，其中 $0 \le r < n$，则余数 $j < r$ 被命中 $q+1$ 次，其他余数被命中 $q$ 次，因此对于均匀分布的状态

$$\frac{P(j)}{1/n} \in \Bigl[\frac{qn}{2^{64}},\ \frac{(q+1)n}{2^{64}}\Bigr], \qquad \Bigl\lvert\frac{P(j)}{1/n} - 1\Bigr\rvert \le \frac{n}{2^{64}} .$$

对于一百万个差值，偏差低于 $10^{-13}$，远低于区间的蒙特卡罗误差。

选择 $B$。 由 $B$ 个重抽样中位数估计的端点，以概率单位计的标准误为 $\sqrt{p(1-p)/B}$ 阶，其中 $p = \alpha/2$。对于 95 % 区间和 $B = 2000$，这约为 2.5 % 尾部质量的 $0.35$ 个百分点。对 95 % 区间请使用 $B \ge 1000$，并把种子和 $B$ 与实验配置一起保存。

离群值视图，而非删除离群值

问题。 报告希望图表中没有那个 50 ms 的缺页尖峰；而决策不能取决于某人选择隐藏了哪些点。选择。 filter_outliers 在显式的 OutlierPolicy 下返回过滤后的副本，运行器和事件流中都不会应用它。理由。 原始数据保持可审计，策略是一个可以被记录的具名值。两种非平凡策略在干净的正态数据上具有已知的误标率：

$$\begin{aligned} \text{Tukey:}&\quad Q(3/4) + 1.5\,\mathrm{IQR} \approx \mu + (0.6745 + 1.5\cdot 1.349)\,\sigma = \mu + 2.698\,\sigma, &\quad P(\lvert Z\rvert > 2.698) &\approx 0.70\,\%, \\ \text{MAD trim:}&\quad 3\,\mathrm{MAD} \approx 3\cdot 0.6745\,\sigma = 2.024\,\sigma, &\quad P(\lvert Z\rvert > 2.024) &\approx 4.3\,\% . \end{aligned}$$

由于 MAD 未缩放，在正态数据上 MAD 截尾大约比 Tukey 围栏激进四倍。当超过一半的值相同时它也会退化：此时 MAD 为 $0$，只有等于中位数的值会保留下来。

错误即值

无效的自助法输入返回 Err(BootstrapError)。否则，空的或含非有限值的样本会产生一个看似合理的全零或 NaN 区间。summarize 和 compare_paired 保持为全函数：它们返回调用者可以检测的中性值（count == 0、Unknown、Invalid）。

正确性与不变量

• 输入不会被修改。 summarize 对副本排序；filter_outliers 返回新数组；自助法对每个重抽样样本排序，从不对输入排序。
• 范围。 由 $Q$ 的范围性质，median、q1 和 q3 位于 $[x_{(1)}, x_{(n)}]$ 中。每个重抽样中位数都位于 $[\min d, \max d]$ 中，因为它是从 $d$ 中抽取的值的分位数；区间端点又是这些中位数的分位数，因此 $\min d \le$ low $\le$ high $\le \max d$。
• 顺序。 对于 $0 < \gamma < 100$，有 $\alpha/2 < 1 - \alpha/2$，由 $Q$ 的单调性得到 low $\le$ high。
• 确定性。 自助法是（值、种子、$B$、$\gamma$）的函数。对有限双精度数排序是确定性的，生成器使用回绕整数算术，因此相同的输入在每个目标上都给出逐位相同的边界。
• 浮点误差。 均值是从左到右的求和。令 $u = 2^{-53}$、$\gamma_k = ku/(1-ku)$，标准界 $\lvert \mathrm{fl}(\sum x_i) - \sum x_i\rvert \le \gamma_{n-1}\sum \lvert x_i\rvert$ 对正的计时值而言就是至多 $\gamma_{n-1} \approx n u$ 的相对误差。方差使用两遍公式 $\sum (x_i - \bar x)^2$，避免了单遍公式 $\sum x_i^2 - n\bar x^2$ 的相消问题。³3 N. J. Higham，Accuracy and Stability of Numerical Algorithms，第 2 版，SIAM，2002，§4.2 与 §1.9。
• 复杂度。 summarize 为 $O(n\log n)$；compare_paired 为 $O(n\log n)$；bootstrap_interval 为 $O(B\,n\log n + B\log B)$；filter_outliers 为 $O(n\log n)$。

被否决的方案

• BCa 与学生化自助法。 它们把覆盖误差降低到 $n^{-1}$ 阶，但需要为每个重抽样样本给出刀切法加速估计或方差估计；对于中位数，在基准测试常见的小 $n$ 下两者都不稳定。未实现。
• 分层自助法。 先重抽样数据集、再重抽样重复，更契合 HierarchicalDatasetsAndRepeats，但比较函数接收的是扁平数组。调用者仍然可以自行按数据集重抽样。
• Hodges–Lehmann 估计量。 两两 Walsh 平均的中位数在对称条件下更有效，但代价为 $O(n^2)$，且更难在报告中解释。
• 假设检验（$t$、Welch、Mann–Whitney）。 如上所述，不用于决策；它们的 $p$ 值回答的是另一个问题。
• 缩放的 MAD。 乘以 1.4826 隐含了正态性假设；因此报告未缩放的值，并在文档中说明该因子。

边界

• 决策只使用点估计和阈值；区间被报告，但不用于决策。没有等价性检验。
• 比较函数接收扁平数组。配对、阶段划分（探索性或验证性）和环境兼容性由调用者负责。
• 自助法区间的单位与差值相同；决策以百分比表示。要比较两者，请除以 $\operatorname{med}(b)$。
• 自助法假设差值是可交换的。轮换顺序使相邻区组相互依赖；区间并未对此建模。
• filter_outliers 从不自动应用，RunProtocol 中的 OutlierPolicy 是记录下来的意图，而不是一个动作。
• 只有自助法函数会拒绝非有限值。

Footnotes

1. R. J. Hyndman 和 Y. Fan，“Sample quantiles in statistical packages”，The American Statistician 50(4)，1996。第 7 型是 R 和 NumPy 的默认定义。 ↩

2. B. Efron 和 R. J. Tibshirani，An Introduction to the Bootstrap，Chapman & Hall，1993，§13.3 与 §14。 ↩

3. N. J. Higham，Accuracy and Stability of Numerical Algorithms，第 2 版，SIAM，2002，§4.2 与 §1.9。 ↩