GSPRT（一般化逐次確率比検定）

前提知識：SPRT、対数尤度比（LLR）

このページの要点

• GSPRT（Generalized SPRT）は、仮説中の未知パラメータを最尤推定量（MLE）で置き換えた SPRT の一般化
• 標準的な SPRT は仮説を完全に特定する必要があるが、GSPRT は部分的な仮説（スコアのみ指定）で動作する
• 誤り率が 0 に近づくとき、GSPRT は標準 SPRT と漸近的に同じ振る舞いをする
• Fishtest が採用している「SPRT」は実際には GSPRT であり、ShogiArena の SPRT 実装の理論的基盤でもある

標準 SPRT の限界

単純仮説と複合仮説

標準的な SPRT（Wald, 1945）は単純仮説を前提とします。 つまり、H0 と H1 のもとでの確率分布が完全に特定されている必要があります。

単純仮説（標準 SPRT）:
  H0: 分布は P0 である（完全に特定）
  H1: 分布は P1 である（完全に特定）

複合仮説（GSPRT）:
  H0: スコアの期待値は s0 である（分布の他のパラメータは不明）
  H1: スコアの期待値は s1 である（分布の他のパラメータは不明）

エンジンテストで立てる仮説は「スコアの期待値が \(s_0\)（または \(s_1\)）である」という形であり、引き分け率などのパラメータは未知のまま残ります。 分布が一つに定まらないため、単純仮説を設定できず、標準 SPRT を直接適用できません。

引き分け率を推定して単純仮説に戻す方法

素朴な解決策として、引き分け率をデータから推定し、それを固定して単純仮説を構成する方法があります（Fishtest のレガシー BayesElo モデル）。 ただし推定した引き分け率を真値のように扱うため、判定にバイアスが生じる可能性があります。

GSPRT は、未知パラメータを固定せず、仮説ごとに尤度を最大化する形で扱うことでこの問題を回避します。

GSPRT の定義

一般化対数尤度比

GSPRT では、仮説中の未知パラメータを制約付き最尤推定量で置き換えます。¹

\[ \Lambda_n = \frac{\sup_{\theta \in \Theta_1} L(\theta; X_1, \ldots, X_n)}{\sup_{\theta \in \Theta_0} L(\theta; X_1, \ldots, X_n)} \]

各記号の意味は次のとおりです。

• \(L(\theta; X_1, \ldots, X_n)\) はデータの尤度関数
• \(\Theta_0, \Theta_1\) はそれぞれ H0, H1 のもとでのパラメータ空間
• \(\sup\) は制約のもとでの最尤推定（MLE）に対応

対数を取ると、次の形になります。

\[ \text{LLR}_n = \ln \Lambda_n = \ell(\hat{\theta}_1) - \ell(\hat{\theta}_0) \]

ここで \(\hat{\theta}_0, \hat{\theta}_1\) はそれぞれ \(\Theta_0, \Theta_1\) 上の制約付き MLE、\(\ell\) は対数尤度関数です。

エンジンテストでの適用

エンジンテストでは、この枠組みが次のように対応します。

• パラメータ空間：多項分布のパラメータ（勝率、引き分け率、負け率）
• 制約：期待スコアが \(s_0\)（H0）または \(s_1\)（H1）であること
• MLE：観測された頻度分布のもとで、指定されたスコアを持つ多項分布のうち最も尤もらしいもの

多項分布の制約付き MLE

定式化

観測された経験分布を \(\hat{p} = (\hat{p}_0, \hat{p}1, \ldots, \hat{p}{N-1})\) とします（\(N\) はカテゴリ数）。 各カテゴリの値を \(a_i = i / (N-1)\)（三項：\(i \in {0, 1, 2}\)、五項：\(i \in {0, 1, 2, 3, 4}\)）とします。

「期待値が \(s\) である多項分布のうち、経験分布に対する尤度を最大化する分布」を求めます。

世俗方程式（Secular Equation）

制約付き MLE は、次の世俗方程式を解くことで得られます。²

\[ \sum_{i} \frac{\hat{p}_i \cdot (a_i - s)}{1 + x \cdot (a_i - s)} = 0 \]

この方程式を \(x\) について解き、MLE 分布を構成します。

\[ p_i^{\text{MLE}} = \frac{\hat{p}_i}{1 + x \cdot (a_i - s)} \]

def MLE_expected(pdfhat, s):
    """期待値が s である制約のもとで、
    経験分布 pdfhat に対する MLE 分布を計算する。"""
    pdf1 = [(ai - s, pi) for ai, pi in pdfhat]
    x = secular(pdf1)  # 世俗方程式を解く
    pdf_MLE = [(ai, pi / (1 + x * (ai - s))) for ai, pi in pdfhat]
    return pdf_MLE


def secular(pdf):
    """世俗方程式 sum(pi * ai / (1 + x * ai)) = 0 を解く。"""
    v = min(ai for ai, pi in pdf)
    w = max(ai for ai, pi in pdf)
    # x の探索範囲: (-1/w, -1/v)
    lower_bound = -1 / w
    upper_bound = -1 / v
    x = brentq(lambda x: sum(pi * ai / (1 + x * ai) for ai, pi in pdf),
               lower_bound + epsilon, upper_bound - epsilon)
    return x

世俗方程式の意味

世俗方程式は、「観測データに最もフィットしつつ、期待値が \(s\) であるという制約を満たす分布」を求めるラグランジュ乗数法の結果です。

パラメータ \(x\) はラグランジュ乗数に対応し、\(x = 0\) のとき MLE は経験分布そのもの（制約なしの場合）になります。

GSPRT-LLR の計算

正確な公式

H0（スコア = \(s_0\)）と H1（スコア = \(s_1\)）に対する GSPRT-LLR は次のように書けます。

\[ \text{LLR} = N \sum_{i} \hat{p}i \cdot \ln\frac{p{1,i}^{\text{MLE}}}{p_{0,i}^{\text{MLE}}} \]

ここで \(p_{0,i}^{\text{MLE}}\), \(p_{1,i}^{\text{MLE}}\) はそれぞれ H0, H1 のもとでの制約付き MLE です。

def LLR(pdf, s0, s1):
    """GSPRT の一般化対数尤度比（N で割ったもの）を計算する。"""
    pdf0 = MLE_expected(pdf, s0)
    pdf1 = MLE_expected(pdf, s1)
    return sum(pi * (math.log(p1i) - math.log(p0i))
               for (_, pi), (_, p0i), (_, p1i) in zip(pdf, pdf0, pdf1))

2 次近似公式

ブラウン運動近似を用いると、GSPRT-LLR は次の式で近似できます。³

\[ \text{LLR} \approx \frac{(s_1 - s_0)(2\hat{s} - s_0 - s_1)}{2 \hat{\sigma}^2 / N} \]

ここで \(\hat{s}\) は観測スコア、\(\hat{\sigma}^2\) はスコアの分散、\(N\) はサンプル数です。

分母の \(N\) を分子に移すと、等価な次の表現になります。

\[ \text{LLR} \approx \frac{N(s_1 - s_0)(2\hat{s} - s_0 - s_1)}{2 \hat{\sigma}^2} \]

def LLR_alt2(pdf, s0, s1):
    """2 次近似による GSPRT-LLR（N で割ったもの）。"""
    s, var = stats(pdf)
    return (s1 - s0) * (2 * s - s0 - s1) / var / 2.0

この近似は、\(\hat{s}\) が \(s_0\) と \(s_1\) の間にあるとき特に精度が高いです。

分散比による近似

\[ \text{LLR} \approx \frac{N}{2} \ln\frac{r_0}{r_1} \]

ここで \(r_j = \sum_i \hat{p}_i (a_i - s_j)^2\) は仮説 \(s_j\) のもとでの残差 2 乗和です。

def LLR_alt(pdf, s0, s1):
    """分散比近似による GSPRT-LLR（N で割ったもの）。"""
    r0, r1 = [sum(prob * (value - s) ** 2 for value, prob in pdf)
              for s in (s0, s1)]
    return 0.5 * math.log(r0 / r1)

ブラウン運動近似

GSPRT をブラウン運動として扱う

GSPRT の LLR 過程は、十分なサンプル数のもとではドリフト付きブラウン運動で近似できます。⁴

\[ \text{LLR}n \approx \mu \cdot n + \sigma{\text{LLR}} \cdot W_n \]

ここで \(W_n\) は標準ウィーナー過程、\(\mu\) と \(\sigma_{\text{LLR}}^2\) はそれぞれ LLR ジャンプのドリフトと分散です。

ドリフトと分散

真のスコアが \(s\) のとき、ブラウン運動近似でのドリフトと分散は次のようになります。

\[ \mu = \frac{(s - (s_0 + s_1)/2)(s_1 - s_0)}{v}, \quad \sigma^2 = \frac{(s_1 - s_0)^2}{v} \]

ここで \(v\) はスコアの分散です。

期待対局数の公式

ブラウン運動近似を用いると、SPRT テストの期待対局数を計算できます。 これは SPRT calculator の基礎となっています。

オーバーシュート補正

離散更新によるオーバーシュート

SPRT の理論的な境界値は、LLR が連続的に変化することを前提としています。 しかし実際の LLR は 1 局または 1 ペアごとに離散的に更新されるため、境界値をちょうど踏むのではなく、飛び越える（オーバーシュートする）ことがあります。

LLR
 +2.94 ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ 理論的境界
       │                 ╱╱
       │              ╱╱╱
       │           ╱╱╱     ← オーバーシュート！LLR が +3.1 で判定
       │        ╱╱╱
 +2.00 ─ ─ ╱╱╱ ─ ─ ─ ─ ─ ─
       └──────────────────→ ゲーム数

オーバーシュートが大きいと、実際の誤り率が設定値（α, β）より大きくなります。

Siegmund のオーバーシュート補正

Fishtest は、Siegmund (1985) の理論に基づく動的オーバーシュート補正を実装しています。⁵ LLR が境界付近にあるとき、オーバーシュートの期待値を推定し、実効的な境界を内側へ調整します。

\[ \text{lower}{\text{eff}} = \text{lower} + o_0, \quad \text{upper}{\text{eff}} = \text{upper} - o_1 \]

オーバーシュートは、これまでのジャンプ幅から次のように推定します。

\[ o_0 = -\frac{\sum \delta_0^2}{2 \sum \delta_0}, \quad o_1 = \frac{\sum \delta_1^2}{2 \sum \delta_1} \]

ここで \(\delta_0\) は LLR が下降して新しい最小値を記録したときのジャンプ、\(\delta_1\) は上昇して新しい最大値を記録したときのジャンプです。

# Fishtest のオーバーシュート補正
if LLR_ < o["ref0"]:
    delta = LLR_ - o["ref0"]
    o["m0"] += delta       # δ の累積和
    o["sq0"] += delta ** 2  # δ² の累積和
    o["ref0"] = LLR_

if LLR_ > o["ref1"]:
    delta = LLR_ - o["ref1"]
    o["m1"] += delta
    o["sq1"] += delta ** 2
    o["ref1"] = LLR_

# 実効的な境界の計算
o0 = -o["sq0"] / o["m0"] / 2 if o["m0"] != 0 else 0
o1 = o["sq1"] / o["m1"] / 2 if o["m1"] != 0 else 0

# 判定
if llr < lower_bound + o0:
    state = "rejected"   # H0 採択
elif llr > upper_bound - o1:
    state = "accepted"   # H1 採択

漸近的性質

GSPRT と SPRT の関係

Li, Liu, and Ying (2014) により、GSPRT が次の漸近的性質を持つことが証明されています。¹

1. 漸近的に SPRT と同等：誤り率 \(\alpha, \beta \to 0\) のとき、GSPRT の判定領域は標準 SPRT の判定領域に収束する
2. 一様な漸近的最適性：任意の真のパラメータ値に対して、期待サンプル数が漸近的に最小となる
3. Elo 差の大きさに非依存：検定される Elo 差が小さい必要はない

実用的な意味

これらの性質により、GSPRT を標準 SPRT と同様に運用できます。

• 判定境界 \(\ln(\beta/(1-\alpha))\) と \(\ln((1-\beta)/\alpha)\) はそのまま適用可能
• \(\alpha = \beta = 0.05\) という標準的な設定で、理論保証が十分に機能する

正規化 Elo での GSPRT

nElo を使う場合、検定統計量は「期待値」ではなく「t 値」になります。⁶

t 値に基づく MLE

\[ t = \frac{\mu - 0.5}{\sigma} \]

この t 値に対する制約付き MLE は、固定小数点法による反復計算で求めます。

def MLE_t_value(pdfhat, ref, s):
    """t 値 = (mu - ref) / sigma が s であるという制約のもとで MLE を計算する。"""
    N = len(pdfhat)
    pdf_MLE = uniform(pdfhat)  # 一様分布で初期化
    for _ in range(10):  # 固定小数点反復
        mu, var = stats(pdf_MLE)
        sigma = var ** 0.5
        pdf1 = [
            (ai - ref - s * sigma * (1 + ((mu - ai) / sigma) ** 2) / 2, pi)
            for ai, pi in pdfhat
        ]
        x = secular(pdf1)
        pdf_MLE = [(pdfhat[i][0], pdfhat[i][1] / (1 + x * pdf1[i][0]))
                   for i in range(N)]
    return pdf_MLE

nElo での LLR 近似

\[ \text{LLR} \approx \frac{N}{2} \ln\frac{1 + (\text{nt} - \text{nt}_0)^2}{1 + (\text{nt} - \text{nt}_1)^2} \]

ここで \(\text{nt}\) は観測された正規化 t 値、\(\text{nt}_0, \text{nt}_1\) は閾値の正規化 t 値です。

三項分布と五項分布の自動選択

Fishtest の実装では、入力データの形式によってモデルが自動的に選択されます。

五項分布ではスコアが 0〜2.0 の範囲を取り、互角に対応する期待値は 0.5 ではなく 1.0 になります。 MLE の計算式は同一で、カテゴリ数とスコアの範囲だけが異なります。

実装リファレンス

参考文献

• Michel Van den Bergh, "The SPRT for a Brownian Motion" — ブラウン運動としての SPRT の理論
• Statistical Methods and Algorithms in Fishtest — Fishtest の統計手法の包括的ドキュメント

次に読む

• 五項分布モデル：ペアゲームを単位とする分析手法