JavaScript での Smirnov-Grubbs 検定 の実装

# Reactive stat 内部の Smirnov-Grubbs 検定の実装 です

外れ値は、データセットの他の値から大きく外れた値のことを指します。
外れ値は、データ分析に大きな影響を与える可能性があるため、適切に検出し、処理することが重要です。

この記事では、Smirnov-Grubbs 検定を用いて外れ値を検出する方法を JavaScript で実装する方法を説明します。

## はじめに

弊社では **[ブラウザだけで使える無料統計ソフト Reactive stat](https://www.emuyn.net/stats)** を提供しています。
信頼性の高い R で統計解析し、その結果を AI が解説します！
その背景には、統計に苦労している医療者の助けになりたい、という気持ちがあります。
最終的な統計解析は R で行うのですが、レスポンスとサーバー負荷の改善のため、一部の統計計算はブラウザ内で行っています。
そのうち、汎用性の高い部分を共有させていただきたいと思います。
できるだけ R の出力と整合性を持つように javascript をインプリメントしてあるので、参考になれば幸いです。

## 使用したライブラリ

このコードでは、jStatライブラリを使用しています。jStatは、JavaScriptで統計計算を行うための便利なライブラリです。
jStatの詳細については、[公式ドキュメント](https://github.com/jstat/jstat) を参照してください。
全関数の一覧は [jStat v1.9.3 Documentation](https://jstat.github.io/all.html) を参照してください。

## Smirnov-Grubbs検定の概要

この検定は、正規分布を仮定して外れ値の検定をより精密に行います。

1. データセット中で最も大きな値と最も小さい値を確認し、平均からの距離が最も大きい値 ($(X_{\text{max}}$ または $(X_{\text{min}}$)) を選択します
2. 次に統計量 $G$ を以下の式で計算します：
    $$ G = \frac{ |X - \bar{X}|}{sd} $$
    ここで、
    - $\bar{X}$ はサンプル平均
    - $sd$ はサンプル標準偏差（標本標準偏差）
    - $X$ は$(X_{\text{max}}$ または $(X_{\text{min}}$)です。
3. $G$ の値が臨界値よりも大きいかどうかを確認します。臨界値は、希望する有意水準とサンプルサイズに依存し、以下の式で計算されます：
    $$ G_{\text{crit}} = \frac{(n-1) \sqrt{t^2}}{\sqrt{n(n-2 + t^2)} } $$
    ここで、
    - $n$ はサンプルサイズ
    - $t$ は自由度 $(n-2)$ のStudent's t分布の $(1-\alpha / (2n))$パーセンタイル（両側検定）です。
4. もし $ G > G_{\text{crit}} $ であれば、そのデータポイントは外れ値と判定されます。そのデータポイントをデータセットから除外し、次のステップへ進みます。

除外後のデータセットに対して、ステップ1から4を繰り返します。これをp値が設定した有意水準より大きくなる（通常は0.05以上）、または外れ値がこれ以上検出されなくなるまで続けます。



## コード

ライブラリ込みの html にしてありますので、
[jsfiddle](https://jsfiddle.net/) などにコピペして簡単に試せます。

```html:html
<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  <title>Smirnov-Grubbs Outlier Detection</title>
  <script src="https://cdn.jsdelivr.net/npm/jstat@latest/dist/jstat.min.js"></script>
</head>
<body>
  <h1>Smirnov-Grubbs Outlier Detection Example</h1>
  <p>Check the console for the result.</p>
  
  <script>
    // Smirnov-Grubbs 検定
    function getSmirnovGrubbsOutliers(data, confLevel = 0.95) {
        let x = data.map(x => Number(x)).filter(val => !isNaN(val));
        let _x = [...x];
        let delVal = [];

        while (true) {
            let n = x.length;
            if (n < 3) {
                break;
            }

            let r = [jStat.min(x), jStat.max(x)];
            let meanX = jStat.mean(x);
            let sdX = jStat.stdev(x, true);
            let t = r.map(val => Math.abs(val - meanX) / sdX);
            let q2 = t.map(val => (n - 2) / (Math.pow((n - 1), 2) / Math.pow(val, 2) / n - 1));
            q2 = q2.map(val => (val < 0 ? 0 : val));
            let q = q2.map(val => Math.sqrt(val));
            let p = q.map(val => n * (1 - jStat.studentt.cdf(val, n - 2)));

            if (t[0] < t[1]) {
                if (p[1] < 1 - confLevel) {
                    delVal.push(r[1]);
                    x = x.filter(val => val !== r[1]);
                    continue;
                }
            } else {
                if (p[0] < 1 - confLevel) {
                    delVal.push(r[0]);
                    x = x.filter(val => val !== r[0]);
                    continue;
                }
            }
            break;
        }
        return {
            outliers: delVal.sort(), // 重複無しのため実際の数以下となる
            removed: _x.filter(x => delVal.includes(x)),
            cleaned: x.sort(),
        };
    }

    // サンプルデータ
    const data = [23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 100, 0];

    // 外れ値の検出
    const sg = getSmirnovGrubbsOutliers(data);
    console.log('Outliers:', sg.outliers);
    console.log('Removed values:', sg.removed);
    console.log('Cleaned data:', sg.cleaned);
  </script>
</body>
</html>
```

## 検証用 R コード

[rdrr](https://rdrr.io/snippets/) にコピペして簡単に試せます。

```R:R
# Smirnov-Grubbs Outlier Detection

# サンプルデータ
data <- c(23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 23, 22, 21, 25, 30, 31, 23, 22, 100, 0)

# 外れ値の検出
outliers <- outliers::grubbs.test(data)

# 結果の表示
print(outliers)
```

## 注意点

- このコードは、jStatライブラリに依存しています。jStatを読み込まないと動作しません
- 浮動小数点数の計算では、わずかな誤差が生じる可能性があります
- Smirnov-Grubbs 検定は、データが正規分布に従っていることを前提としています。データが正規分布に従っていない場合、結果が信頼できないものになる可能性があります
- R の grubbs.test() 関数は、一度に1つの外れ値しか検出しません。複数の外れ値を検出するには、外れ値を除去した後、再度検定を実行する必要があります
- JavaScript のコードでは、信頼水準 confLevel を設定できます。デフォルトは0.95（95%）です。R の grubbs.test() 関数では、信頼水準を設定するオプションはなく、外れ値として選定されたデータの p値を確認することで、外れ値として扱うかどうかを決定します。

## 最後に

コードの内容はよく吟味し、R との食い違いが極力ないように気を付けていますが、もし間違いに気づかれた場合には指摘していただけますとありがたいです。

ただ、Chat GPT に読ませてその出力を検証もせず「問題があるので修正します」として貼り付けられると、他の方を混乱させてしまうのと、弊社の製品に対する信頼を損ねることにもつながりかねませんので、どうかそのようなことのないようにお願いいたします。

検証のための R のコードもつけてありますので、ご活用いただけますと幸いです。