stats の設計

設計目標

ベンチマークの計時値は歪んでおり、裾が重く、計測された時点と相関しています。stats はこれらの性質に耐える判断を提供します。頑健な位置の推定、対応のある計測に基づく比較、明示的な実用的な閾値、そしてランダム性がシードで固定された区間です。すべての関数は配列の純粋な変換なので、判断はイベントストリームに保存された生の観測から再計算できます。

数学的背景

順序統計量と標本分位点

$x_1, \dots, x_n$ を標本、$x_{(1)} \le \dots \le x_{(n)}$ をその順序統計量とします。このパッケージはあらゆる箇所で 1 つの分位点の定義、すなわち Hyndman と Fan のタイプ 7 の線形補間を使います:¹1 R. J. Hyndman and 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 .$$

この式から直ちに 3 つの性質が導かれ、以下で使われます。

1. 値域。 $Q(p)$ は隣り合う 2 つの順序統計量の凸結合なので、$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. アフィン同変性。 $b > 0$ として $y_i = a + b\,x_i$ とすると、ソートはこの写像と可換で、凸結合はアフィン写像と可換なので、$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 は 2 系統の推定量を並べて報告します:

古典的	頑健
$\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$（1 つの値で十分）、四分位数と IQR では $1/4$、中央値と MAD では $1/2$ です。スケジュールから外された 1 回の実行はベンチマークの平均を動かしますが、実行の半数が影響を受けない限り中央値は動かせません。

$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$ の正規分布に従うデータでは、$z_{3/4} = \Phi^{-1}(3/4) \approx 0.6745$ として $\operatorname{med}\lvert X - \mu\rvert = z_{3/4}\,\sigma$ なので、$1.4826\cdot\mathrm{MAD}$ が $\sigma$ を推定し、$\mathrm{IQR} \approx 1.349\,\sigma$ となります。

対応のある計測

ランナーは各実装をブロックごとに 1 回、ローテーションした順序で計測します（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$ を要約し、2 つの標本を別々に要約することはありません。

パーセンタイル・ブートストラップ

$\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].$$

これは次の条件の下で厳密です。$W = \varphi(\hat\theta) - \varphi(\theta)$ が $0$ について対称な分布 $H$ を持ち、再標本化の下での $\varphi(\hat\theta^*) - \varphi(\hat\theta)$ も同じ $H$ で記述されるような増加写像 $\varphi$ が存在すると仮定します。このとき $\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 and R. J. Tibshirani, An Introduction to the Bootstrap, Chapman & Hall, 1993, §13.3 and §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$ で、アフィン同変であり、平均と中央値の距離はレポートで確認すべき裾の存在を示すからです。

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}.$$

どちらも同じ 2 つの中央値の関数なので、食い違うことはありません。$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）、点推定値に対する実用的な閾値。選択。 判断は次の 3 分岐の規則です

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

これは Invalid（長さの不一致）と Unknown（ペアなし）の後に検査されます。理由。 閾値は、ユーザーが問うている問い（「少なくとも 2 % 速いか？」）を判断の単位で述べるものであり、規則は 2 つの中央値から再現でき、反復を増やすことに報いることもありません。$t > 0$ では 3 つの領域が $\mathbb R$ を分割します。$t = 0$ では最初の 2 つが $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}$ を法として可逆なので、1 ステップは $0$ を固定する 64 ビットワードの置換です。シード 0 が定数 88172645463393265 に置き換えられるのはそのためです。生成器はラップアラウンドする 64 ビット整数算術だけを使うため、ストリームはどのターゲットでも同一です。

状態 $x$ は $j = x \bmod n$ によってインデックスに写されます。$0 \le r < n$ として $2^{64} = qn + r$ と書くと、剰余 $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}} .$$

100 万個の差でも偏りは $10^{-13}$ 未満で、区間のモンテカルロ誤差よりはるかに小さいです。

$B$ の選び方。 $B$ 個の再標本化した中央値から推定した端点は、$p = \alpha/2$ として、確率の単位で $\sqrt{p(1-p)/B}$ のオーダーの標準誤差を持ちます。95 % 区間と $B = 2000$ では、これは 2.5 % の裾の確率質量に対しておよそ $0.35$ パーセントポイントです。95 % 区間には $B \ge 1000$ を使い、シードと $B$ を実験の設定とともに保存してください。

外れ値の削除ではなく外れ値のビュー

問題。 レポートでは 50 ms のページフォールトのスパイク 1 つに邪魔されないプロットが欲しい一方で、判断が、誰かが隠すことにした点に依存してはなりません。選択。 filter_outliers は明示的な OutlierPolicy の下でフィルタリングしたコピーを返し、ランナーやイベントストリームでそれを適用するものは何もありません。理由。 生のデータは監査可能なまま残り、ポリシーは記録できる名前付きの値になります。2 つの自明でないポリシーは、きれいな正規分布のデータに対して既知の誤検出率を持ちます:

$$\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 のフェンスのおよそ 4 倍積極的です。また、値の半数を超えるものが一致すると退化します。そのとき MAD は $0$ となり、中央値に等しい値だけが残ります。

値としてのエラー

不正なブートストラップの入力は Err(BootstrapError) を返します。そうしなければ、空の標本や有限でない標本から、ゼロや NaN からなるもっともらしく見える区間が得られてしまいます。summarize と compare_paired は全域関数のままです。呼び出し側が検査できる中立的な値（count == 0、Unknown、Invalid）を返します。

正しさと不変条件

• 入力は変更されません。 summarize はコピーをソートし、filter_outliers は新しい配列を返し、ブートストラップは各再標本をソートして入力はソートしません。
• 値域。 $Q$ の値域の性質により、median、q1、q3 は $[x_{(1)}, x_{(n)}]$ に収まります。再標本化した中央値はどれも $d$ から引いた値の分位点なので $[\min d, \max 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$ の相対誤差になります。分散は 2 パスの式 $\sum (x_i - \bar x)^2$ を使い、1 パスの $\sum x_i^2 - n\bar x^2$ で生じる桁落ちを避けます。³3 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, §4.2 and §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 and Y. Fan, “Sample quantiles in statistical packages”, The American Statistician 50(4), 1996. タイプ 7 は R と NumPy の既定です。 ↩

2. B. Efron and R. J. Tibshirani, An Introduction to the Bootstrap, Chapman & Hall, 1993, §13.3 and §14. ↩

3. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, §4.2 and §1.9. ↩