Brunner-Munzel検定のPython実装
Brunner-Munzel検定をPython（NumPy/SciPy使用）の実装です。
下記bmtest関数に2群のデータを入力すると、検定統計量W、p値（両側）、自由度、および帰無仮説で示される確率の推定値を出力します。
片側検定にする場合（どちらにずれているかが重要な場合）、出力されたp値を0.5倍して、Wの値が正負で判断します。
注意: 少しもかぶりもしない2群を入力すると、分散が0となって、検定統計量Wや自由度の計算中に Runtimewarning; divide by zero が出ます。
# coding: utf-8
#
# bmtest(x1, x2)
# Calculate Brunner-Munzel-test scores.
#
# Parameters:
# x1, x2: array_like
# Numeric data values from sample 1, 2.
#
# Returns:
# w:
# Calculated test statistic.
# p_value:
# Two-tailed p-value of test.
# dof:
# Degree of freedom.
# p:
# "P(x1 < x2) + 0.5 P(x1 = x2)" estimates.
#
# References:
# * https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
import numpy
import scipy.stats
def bmtest(x1, x2):
if isinstance(x1, numpy.ndarray): x1 = x1.tolist()
if isinstance(x2, numpy.ndarray): x2 = x2.tolist()
n1, n2 = len(x1), len(x2)
R = scipy.stats.rankdata(x1 + x2, method = 'average')
R1, R2 = R[:n1], R[n1:]
r1_mean, r2_mean = numpy.mean(R1), numpy.mean(R2)
p = (r2_mean - r1_mean) / (n1 + n2) + 0.5
Ri1 = scipy.stats.rankdata(x1, method = 'average')
Ri2 = scipy.stats.rankdata(x2, method = 'average')
var1 = numpy.var([r - ri for r, ri in zip(R1, Ri1)], ddof = 1)
var2 = numpy.var([r - ri for r, ri in zip(R2, Ri2)], ddof = 1)
w = ((n1 * n2) * (r2_mean - r1_mean)) / ((n1 + n2) * numpy.sqrt(n1 * var1 + n2 * var2))
dof = (n1 * var1 + n2 * var2) ** 2 / ((n1 * var1) ** 2 / (n1 - 1) + (n2 * var2) ** 2 / (n2 - 1))
c = scipy.stats.t.cdf(abs(w), dof) if not numpy.isinf(w) else 0.0
p_value = min(c, 1.0 - c) * 2.0
return (w, p_value, dof, p)
JUGEMテーマ：コンピュータ
