Mann-Whitneyの $U$ 検定
ここでいうMann-Whitneyの $U$ 検定は、統計のところで書いたWilcoxon-Mann-Whitney検定（順位和検定）と同じものです。PythonではSciPyの stats.mannwhitneyu で行うことが多いので、こういう名前の見出しにしました。
Wilcoxon-Mann-Whitney検定と同じ問題を解いてみましょう。Rでは、正規分布近似で $p = 0.007741$、厳密な検定で $p = 0.00669$ でした。
まずは正規分布近似です：
import numpy as np from scipy import stats x = [1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1] y = [3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4] stats.mannwhitneyu(x, y)
MannwhitneyuResult(statistic=np.float64(32.5), pvalue=np.float64(0.007740786705107166))
$U = 32.5$、$p = 0.007741$ で、Rと同じ結果になりました。
厳密なテストをしてみましょう。でも、これは間違いです：
stats.mannwhitneyu(x, y, method="exact")
MannwhitneyuResult(statistic=np.float64(32.5), pvalue=np.float64(0.015207520078969802))
紛らわしいのですが、method="exact" はタイがないときにしか使えません。タイがあるときに使うと警告もなく間違った答えを返してしまいます💢。ちなみにデフォルトは、片方の大きさが8以下かつタイのないときに厳密な方法を使い、それ以外は正規分布近似になります。
この問題はタイがあるので method="exact" は使えず、次のように permutation test を使います。n_resamples=np.inf で全部の置換を試すので厳密になります。
stats.mannwhitneyu(x, y, method=stats.PermutationMethod(n_resamples=np.inf))
MannwhitneyuResult(statistic=np.float64(32.5), pvalue=np.float64(0.00872571454210975))
あれ？ Rの結果 $p = 0.00669$ と少し違いますね。
片側検定にしてみましょう：
stats.mannwhitneyu(x, y, alternative="less",
method=stats.PermutationMethod(n_resamples=np.inf))
MannwhitneyuResult(statistic=np.float64(32.5), pvalue=np.float64(0.004362857271054875))
ちょうど先ほどの $p=0.00872571454210975$ の半分の結果となりました。つまり、stats.mannwhitneyu() の両側検定は単に片側検定を2倍にしているだけです。一方で、Rは $|U - E(U)| \geq |U_{\text{obs}} - E(U)|$ となる確率を求めます（詳細はWilcoxon-Mann-Whitney検定参照）。
両側検定の定義の問題を避けたいなら、片側検定をして、$p$ 値を 0.05 でなく 0.025 と比較するのもよさそうです。
なお、両群のサイズが大きいと n_resamples=np.inf は使えません。適当な値にしてください。
rng = np.random.default_rng(12345) stats.mannwhitneyu(x, y, method=stats.PermutationMethod(n_resamples=100000, rng=rng))
PythonでRと同じ両側検定をする例（GPT-5.4作）：
import numpy as np
from scipy.stats import permutation_test, rankdata
x = np.array([1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1], dtype=float)
y = np.array([3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4], dtype=float)
def mann_whitney_u(a, b):
z = np.r_[a, b]
r = rankdata(z, method="average") # average ranks for ties
n1 = len(a)
return r[:n1].sum() - n1 * (n1 + 1) / 2
# exact permutation distribution under H0
res = permutation_test(
(x, y),
mann_whitney_u,
permutation_type="independent",
n_resamples=np.inf, # exact
alternative="less" # ignored for our manual two-sided calculation
)
u_obs = res.statistic
u_null = res.null_distribution
u_mean = len(x) * len(y) / 2
# R-style / exactRankTests-style two-sided p-value
p_r_style = np.mean(np.abs(u_null - u_mean) >= np.abs(u_obs - u_mean))
print("U =", u_obs)
print("R-style exact two-sided p =", p_r_style)
SciPyに依存しないR互換の実装例（GPT-5.4作）：
from math import comb
from collections import defaultdict, namedtuple
MannWhitneyURStyleResult = namedtuple(
"MannWhitneyURStyleResult",
["statistic", "pvalue"]
)
def _average_scaled_ranks(values):
"""
values の average rank を 2 倍した整数で返す。
例: rank 3.5 -> 7
"""
n = len(values)
order = sorted(range(n), key=lambda i: values[i])
ranks2 = [None] * n
i = 0
while i < n:
j = i + 1
v = values[order[i]]
while j < n and values[order[j]] == v:
j += 1
# 1-based ranks: i+1, ..., j
# その平均 rank の 2 倍
avg_rank_times_2 = (i + 1) + j
for k in range(i, j):
ranks2[order[k]] = avg_rank_times_2
i = j
return ranks2
def _exact_rank_sum_distribution(x, y):
"""
帰無仮説の下での exact な rank-sum 分布を返す。
ties は average rank を使って条件付き exact 分布として扱う。
Returns
-------
dist_w2 : dict
key = 2 * W (W は x 側の rank sum)
value = その値をとる組合せ数
ranks2 : list[int]
pooled data の各要素の 2 * average rank
"""
x = list(x)
y = list(y)
pooled = x + y
n1 = len(x)
ranks2 = _average_scaled_ranks(pooled)
# 同じ rank 値ごとにまとめる
multiplicities = defaultdict(int)
for r2 in ranks2:
multiplicities[r2] += 1
# DP:
# dp[k][s] = 「ここまで見た rank 群から k 個選んで rank-sum(×2)=s にする方法数」
dp = [defaultdict(int) for _ in range(n1 + 1)]
dp[0][0] = 1
seen = 0
for r2, count in sorted(multiplicities.items()):
new_dp = [defaultdict(int) for _ in range(n1 + 1)]
for k_prev in range(min(n1, seen) + 1):
for s_prev, ways_prev in dp[k_prev].items():
max_take = min(count, n1 - k_prev)
for q in range(max_take + 1):
new_dp[k_prev + q][s_prev + q * r2] += ways_prev * comb(count, q)
dp = new_dp
seen += count
return dp[n1], ranks2
def mannwhitneyu_rstyle(x, y, alternative="two-sided"):
"""
R-style exact Mann-Whitney / Wilcoxon rank-sum test.
Parameters
----------
x, y : sequence of numbers
2標本
alternative : {"two-sided", "less", "greater"}
- "less" : x が y より小さい方向
- "greater" : x が y より大きい方向
- "two-sided": R-style の central two-sided exact p-value
Returns
-------
MannWhitneyURStyleResult(statistic, pvalue)
statistic は x 側の U
"""
x = list(x)
y = list(y)
if len(x) == 0 or len(y) == 0:
raise ValueError("x and y must both be non-empty")
if alternative not in {"two-sided", "less", "greater"}:
raise ValueError("alternative must be 'two-sided', 'less', or 'greater'")
n1 = len(x)
n2 = len(y)
dist_w2, ranks2 = _exact_rank_sum_distribution(x, y)
# 観測された x 側の rank sum W と U
w2_obs = sum(ranks2[:n1]) # 2 * W
u2_obs = w2_obs - n1 * (n1 + 1) # 2 * U
total = comb(n1 + n2, n1)
if alternative == "less":
p = sum(count for w2, count in dist_w2.items() if w2 <= w2_obs) / total
elif alternative == "greater":
p = sum(count for w2, count in dist_w2.items() if w2 >= w2_obs) / total
else:
# R-style two-sided:
# P(|U - E[U]| >= |u_obs - E[U]|)
mean_u2 = n1 * n2 # 2 * E[U]
dev_obs = abs(u2_obs - mean_u2)
p = sum(
count
for w2, count in dist_w2.items()
if abs((w2 - n1 * (n1 + 1)) - mean_u2) >= dev_obs
) / total
return MannWhitneyURStyleResult(statistic=u2_obs / 2, pvalue=p)
実行例：
x = [1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1] y = [3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4] print(mannwhitneyu_rstyle(x, y)) print(mannwhitneyu_rstyle(x, y, alternative="less")) print(mannwhitneyu_rstyle(x, y, alternative="greater"))
MannWhitneyURStyleResult(statistic=32.5, pvalue=0.006690222999955131) MannWhitneyURStyleResult(statistic=32.5, pvalue=0.004362857271054875) MannWhitneyURStyleResult(statistic=32.5, pvalue=0.9978420155247454)
Wilcoxon-Mann-Whitney検定で紹介した漸化式（タイのない場合に厳密）：
def p(n, m, u):
if n * m == 0:
return int(u == 0)
else:
return (n * p(n-1, m, u-m) + m * p(n, m-1, u)) / (n+m)
正規分布との比較：
import matplotlib.pyplot as plt from scipy import stats import numpy as np n = 9 m = 9 x = np.arange(n * m + 1) a = [p(n, m, u) for u in x] b = stats.norm.pdf(x, loc=n*m/2, scale=np.sqrt(n*m*(n+m+1)/12)) plt.plot(x, a) # 漸化式 plt.plot(x, b) # 正規分布近似