2025年にSchüürhuis, Konietschke, Brunnerらによって提案されたノンパラメトリック検定（シューアハウス・コニーチュケ・ブルンナー検定）を紹介します。原著論文の中では $C^2$ 検定と表記されていましたが、カイ二乗検定との混同を避けるため、本記事ではSKB検定と書くことにします。
簡単な説明
SKB検定は、2000年に発表されたブルンナー・ムンチェル検定（BM検定）の改良版です。BM検定はマン・ホイットニーのU検定より優れたノンパラメトリック検定ですが、有意水準を $\alpha=0.005$ にすると精度が悪くなるという問題がありました。SKB検定は $\alpha=0.005$ でも正しく動作するだけでなく、$\alpha=0.05$ のときもBM検定より少し良くなることが数値シミュレーションによって示されています。ただし、SKB検定は比較したい2群の少なくとも一方のデータ数が15未満の場合に非推奨とされているので、注意して下さい。
SKB検定の式はかなり複雑ですが、実装するために最低限必要なことを以下にまとめます。
2群のデータ $x_1,\ldots,x_n$ と $y_1,\ldots,y_m$ が与えられたとします。$n$ は $x$ の要素数、$m$ は $y$ の要素数です。このとき、以下のような $n\times m$ の行列 $A$ を用意します。
$$A_{ij}=\left\{\begin{array}{lc}1&\text{if}\quad x_i<y_j,\\1/2&\text{if}\quad x_i=y_j,\\0&\text{if}\quad x_i>y_j.\end{array}\right.$$
$A$ の平均値を計算します。これは確率的優位性 (stochastic superiority) と呼ばれています。この値が1に近ければ $y$ が優位、0に近ければ $x$ が優位と見なされます。U検定やBM検定でも使われている統計量です。
$$\hat\theta=\frac{1}{nm}\sum_{i=1}^n\sum_{j=1}^mA_{ij}.$$
この $\hat\theta$ の値は確率的に揺らぎます。その分散を見積もったものが以下の式です。
$$s^2=\frac{(n-1)s_x^2+(m-1)s_y^2-nm\hat\theta(1-\hat\theta)+t_{xy}/4}{n(n-1)m(m-1)}.$$
ここで $s_x^2$ は行列 $A$ の行和の分散、$s_y^2$ は列和の分散、$t_{xy}$ は $x$ と $y$ の間のタイ（同じ値が出ること）の数です。詳しくは以下のソースコードを見て下さい。
検定統計量は以下で計算できます。
$$C=\frac{\hat\theta-1/2}{s}\cdot 2\sqrt{\hat\theta(1-\hat\theta)}.$$
これが標準正規分布 $\mathcal{N}(0,1)$ に従うので、その生存関数 $S$ を使って $p$ 値が出てきます。
$$p=2S(|C|).$$
もし $s=0$ だったときは $C=\sqrt{\min(n,m)}$ に置き換えてから $p$ 値を計算します。
ソースコード
import numpy as np
from scipy import stats
def SKB_test(x, y):
x_arr = np.asarray(x)
y_arr = np.asarray(y)
n = len(x_arr) # xの長さ
m = len(y_arr) # yの長さ
A = np.heaviside(y_arr - x_arr.reshape(-1,1), 0.5) # n x m 行列
theta = A.mean() # Aの平均
sx2 = A.sum(axis=1).var(ddof=1) # 行和の分散
sy2 = A.sum(axis=0).var(ddof=1) # 列和の分散
t_xy = (A==0.5).sum() # タイの数
s2 = ((n-1)*sx2 + (m-1)*sy2 - n*m*theta*(1-theta) + t_xy/4) /\
(n*(n-1)*m*(m-1)) # 不偏分散
if s2 > 0:
C = ((theta-0.5)/s2**0.5) * 2*(theta*(1-theta))**0.5 # 検定統計量
else:
C = np.min([n,m])**0.5
p = 2 * stats.norm().sf(np.abs(C)) # p値
return C, p
適用例
元論文に載っていたサンプルデータに適用してみます。
def check_SKB_test():
# 元論文に載っていたサンプルデータ
data1 = [1]*16 + [2]*5 + [4]*1
data2 = [1]*4 + [2]*1 + [3]*5 + [4]*7 + [5]*2
# SKB検定を実行
C, p = SKB_test(data1, data2)
# 結果を表示
print(f'C^2 = {C**2:.1f}')
print(f'p = {p:.6f}')
以下が出力結果です。元論文と同じ値になりました。
C^2 = 14.9
p = 0.000116
補足
元論文では行列 $A$ を使わず、$\hat\theta$, $s_x^2$, $s_y^2$, $t_{xy}$, $s^2$ を順位を用いて定義していました。等価です。理由を説明すると長くなるので、気になる方は 私のarXivの解説 を見て下さい。
元論文では $C$ でなく $C^2$ を検定統計量としていました。$C\sim\mathcal{N}(0,1)$ なので $C^2\sim\chi^2_1$ です。1は自由度です。ここからC2検定という名前になったようです。上記の説明では簡単のため $C$ を使いました。
[2026/09/28追記] 少なくとも一方のグループのデータ数が15未満だった場合、並び替えBM検定 (Neubert & Brunner, 2007) の方がマシなようです。また、両方とも30個以上あり、データ分布が極端に変な形をしていなければ、中心極限定理によりウェルチのt検定が使えます。図で表すと以下のような感じです。
[2026/09/28追記] 95%信頼区間に関しては注意が必要です。元論文では、SKB検定に対応する95%信頼区間の方がBM検定に対応する95%信頼区間より優れていると主張していますが、logit変換したものとの比較を載せていません。私が調べたところ以下のような結果になりました。$n=m=15$, $x\sim\mathcal{N}(0,1)$, $y\sim\mathcal{N}(\mu,1)$, $\mu=\{0,0.5,1,1.5,2,2.5\}$ としました。試行回数は $10^4$ です。これを見る限り、logit変換の方が良さそうです。
[2026/09/28追記] SKB検定はBM検定と比べて主に2か所変わっています。一つは $\hat\theta$ の分散の推定量 $s^2$ の式です。もう一つはスケール因子 $2\sqrt{\hat\theta(1-\hat\theta)}$ を掛けたことです。改善に寄与しているのは主に後者のようです。下図は、$n=m=15$, $x,y\sim\mathcal{N}(0,1)$ で、$\hat\theta$ とその分散の推定量の関係を表したものです。上側の図を見ると、$\hat\theta$ が $0.5$ から離れるほど、分散が小さく推定されてしまっています。真の値は一つなので、これは良くない挙動です。この依存性は概ね二次関数の形をしています。また、この時点でBM検定とSKB検定の差は大きくないことも分かります。下側の図はスケール因子で補正したものです。$\hat\theta$ の位置によらず、その分散の推定値が同程度の値になっていることが分かります。つまり、もともと分散が過小評価されていたケースがあり、それを補正したため、全体としても改善したと考えて良さそうです。
