uechi.io/source/_posts/2019/welch-t-test.md

---
title: プログラムの速度改善が誤差かどうかを統計的に調べる
date: 2019-10-03 17:21:00 +09:00
---

**Welch の t 検定**を用いて 2 つのベンチマークの分布の平均が等しい（速度差は誤差の範疇）か、あるいは異なる（=有意な速度改善が成されている）かどうかを判定したい。

ベンチマーク用に TypeScript プログラムを用意した。

```ts a.ts
function a() {
  const noise = Math.random() - 0.5;
  const offset = 1.0;
  const t = noise * 2 + offset;
  setTimeout(() => console.log(t), t * 1000);
}
a();
```

```ts b.ts
function b() {
  const noise = Math.random() - 0.5;
  const offset = 2.0;
  const t = noise * 2 + offset;
  setTimeout(() => console.log(t), t * 1000);
}
b();
```

まず[hyperfine](https://github.com/sharkdp/hyperfine)で 2 つの プログラムのベンチマークを取り、`result.json`に保存する。

```shell
hyperfine 'ts-node a.ts' 'ts-node b.ts' -r 50 --warmup 3 --export-json ab.json
```

`result.json`の中身は以下のようになる。

```json result.json
{
  "results": [
    {
      "command": "ts-node a.ts",
      "mean": 1.9369869248950002,
      "stddev": 0.6074252496423262,
      "median": 2.005230080295,
      "user": 1.549546345,
      "system": 0.08031985000000001,
      "min": 0.8807363742950001,
      "max": 2.830435366295,
      "times": [
        1.4010462692949999,
        2.830435366295,
        1.010024359295,
        1.159667609295,
        1.8311979602950001,
        ...
      ]
    },
    {
      "command": "ts-node b.ts",
      "mean": 2.833931665055,
      "stddev": 0.6505564501747996,
      "median": 2.7373719187950005,
      "user": 1.5474132649999999,
      "system": 0.07978893000000001,
      "min": 1.938184970295,
      "max": 3.946562622295,
      "times": [
        2.2806011012950003,
        2.0140897212950004,
        2.1835023382950003,
        2.304886362295,
        3.8122057912950003,
        ...
      ]
    }
  ]
}
```

> t 検定はサンプルが正規分布に従っているという仮定を置いているため、大数の法則から本当はもっと試行回数を増やした方が良い。

この`result.json`の`times`配列を受け取り、2 つの分布間に有意差があるかどうかを判定する。

```ts
import fs from "fs";
import { jStat } from "jstat";

const log = console.log;

const sum = (x: number[]) => x.reduce((a: number, b: number) => a + b); // 総和
const sqsum = (x: number[], mu: number) =>
  x.reduce((a: number, b: number) => a + (b - mu) ** 2); // 自乗誤差の総和

function ttest(X: number[], Y: number[]) {
  const Xn = X.length; // サンプル数
  const Yn = Y.length;
  log(`Xn = ${Xn}`);
  log(`Yn = ${Yn}`);

  const X_mu = sum(X) / Xn; // 平均
  const Y_mu = sum(Y) / Yn;
  log(`X_mu = ${X_mu}`);
  log(`Y_mu = ${Y_mu}`);

  const X_sigma = sqsum(X, X_mu) / (Xn - 1); // 不偏分散
  const Y_sigma = sqsum(Y, Y_mu) / (Yn - 1);
  log(`X_sigma = ${X_sigma}`);
  log(`Y_sigma = ${Y_sigma}`);
  const t = (X_mu - Y_mu) / Math.sqrt(X_sigma / Xn + Y_sigma / Yn); // t値
  log(`t = ${t}`);
  const df =
    (X_sigma + Y_sigma) ** 2 /
    (X_sigma ** 2 / (Xn - 1) + Y_sigma ** 2 / (Yn - 1)); // 自由度
  log(`df = ${df}`);
  return jStat.studentt.cdf(-Math.abs(t), df) * 2.0; // p値
}

const filename = process.argv.slice(2)[0];
const result = JSON.parse(fs.readFileSync(filename).toString());
const X = result.results[0].times;
const Y = result.results[1].times;
const p = ttest(X, Y);
log(`p = ${p}`);
log(`p < 0.05 = ${p < 0.05}`);
log(p < 0.05 ? "Possibly some difference there" : "No difference");
```

ここで`X_mu`は分布 X の平均、`X_sigma`は分布 X の不偏分散だ。

$$
\begin{eqnarray}
\mu_X &=& \frac{1}{n_X} \sum^{n_X}_i X_i\\
\sigma_X &=& \frac{1}{n_X-1}\sum^{n_X}_i (X_i - \mu_X)^2
\end{eqnarray}
$$

これを X と Y 両方に対して求めます。さらに以下のようにして t を求める。

$$
t = \frac{\mu_X - \mu_Y}{\sqrt{\frac{\sigma_X}{n_X} + \frac{\sigma_Y}{n_Y}}}
$$

t 分布の累積密度関数 (Cumlative Distribution Function; CDF) を定義する。面倒すぎたので[jstat](https://github.com/jstat/jstat)の`studentt.cdf`を使った。コードを見ると、分子の積分は[シンプソンの公式](https://ja.wikipedia.org/wiki/シンプソンの公式)を使って近似していた。

$$
\text{CDF} =\frac{\int_0^{\frac{v}{t^2+v}}\frac{r^{\frac{v}{2}-1}}{\sqrt{1-r}}dr{}}{\text{exp}(\ln(\Gamma(\frac{v}{2}))+\ln(\Gamma(0.5))+\ln(\Gamma(\frac{v}{2}+0.5)))}
$$

CDF を用いて p 値を求める。両側検定をするので 2 を掛ける。t 分布の自由度 (degree of freedom; df) は$n-1$なので、両分布の自由度を$n_X+n_Y-2$で与える。本当は

$$
\text{df} = \frac{(\sigma_X + \sigma_Y)^2}{
    \frac{\sigma_X^2}{n_X - 1} + \frac{\sigma_Y^2}{n_Y - 1}}
$$

で求める必要があるが、さぼって近似した。

$$
p = \text{CDF}(-|t|, n_X+n_Y-2) \times2
$$

## 結果

異なる実行時間を示すプログラム`a`,`b`を比較すると、2 つの分布の平均が異なることが示唆された。

```
❯ ts-node test.ts ab.json
Xn = 10
Yn = 10
X_mu = 1.8022945422950003
Y_mu = 2.9619571628950006
X_sigma = 0.6067285795623545
Y_sigma = 0.6593856215802901
t = -3.2590814831310353
df = 17.968919419652778
-0.0001571394779906754
p = 0.004364964634417297
p < 0.05 = true
Possibly some difference there
```

p 値が 0.05 未満となり、帰無仮説「２つの分布は等しい」が棄却されたので「２つの分布は等しくない」ことがわかった。では、同じプログラム同士でベンチマークを取るとどうなるか？

```
❯ ts-node test.ts aa.json
Xn = 10
Yn = 10
X_mu = 1.7561671737900002
Y_mu = 1.9892996860899999
X_sigma = 0.5127362786380443
Y_sigma = 0.442053230382934
t = -0.754482245774979
df = 17.901889803947558
-27359.526584574112
p = 0.4603702896905685
p < 0.05 = false
No difference
```

p 値が 0.05 未満ではないため、帰無仮説は棄却されず、つまり「２つの分布は等しい」ことがわかった。

ウェルチの t 検定はスチューデントの t 検定と違って等分散性（２つの分布の分散が等しいこと）を仮定しないため、とても取り扱いやすい検定だ。もちろん等分散性のある分布でも使用できるので、基本的にはウェルチの方法を使う方針で良さそうだ。

## 参考文献

- [Welch's t-test - Rosetta Code](https://rosettacode.org/wiki/Welch%27s_t-test)