MMD（最大平均差異）とは？カーネルで分布の違いを測る

学習時に集めたデータと、いま本番環境に流れ込んでいるデータは、本当に「同じ分布」から来ているのでしょうか。センサーを1台交換した、撮影する季節が変わった、ユーザー層が入れ替わった——こうした変化はモデルの精度を静かに蝕みますが、平均や分散を眺めているだけでは気づけません。平均も分散もぴったり一致したまま、分布の形だけが二峰性に変わっている、ということが現実には起こるからです。

MMD（Maximum Mean Discrepancy、最大平均差異） は、この「2つのデータ集合は同じ分布から来たのか」という問いに、密度推定を一切せずに答えるための道具です。カーネルを1つ選ぶだけで、平均も分散も高次モーメントも、さらには分布の形そのものの違いまでを、たった1つの実数にまとめてくれます。

MMDのアイデアは、この図の3段階に尽きます。まず①で分かるとおり、私たちが持っているのは2つのサンプル集合だけで、その背後にある密度関数は見えません。そこで②のように各データ点を「山（カーネル関数）」に変えて平均し、分布そのものを1本の関数にまとめます。この関数は無限次元空間の1点とみなせるので、③のように2つの分布の違いが単なる「2点間の距離」になります。密度推定もヒストグラムの箱分けも一切登場しない、というのがMMDの出発点です。

MMDが活きる場面は大きく2つあります。1つは 二標本検定（two-sample test） です。「新旧2つの製造ラインの計測値は同じ分布か」「生成モデルが吐いたサンプルは本物のデータと見分けがつくか」を統計的に判定できます。もう1つは ドメイン適応（domain adaptation） です。学習済みモデルの中間特徴について「訓練ドメインの特徴分布」と「本番ドメインの特徴分布」のMMDを損失に足し込むと、ネットワークは両者を揃えるように特徴を学び、ラベルのない本番データでも精度が保たれるようになります。ほかにも、分布ドリフト検知、因果探索で使われる独立性検定（HSIC）、生成モデルの評価指標など、MMDは現代の機械学習のあちこちに顔を出します。

本記事の内容

• 「分布の違いを測る」という問題を、なぜ密度推定なしで解きたいのか
• 平均埋め込み $\mu_P$ の定義と、MMDのカーネル期待値による3項表現の導出（省略なし）
• 特性カーネルなら $\mathrm{MMD}=0 \iff P=Q$ になる理由（Bochnerの定理から）
• 不偏なU統計量とバイアスのあるV統計量の違い、帯域幅の選び方
• 置換検定によるp値の計算と、t検定との検出力比較
• ドメイン適応の損失としてMMDを使う実験

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• カーネル法の理論（カーネルトリック）を導出して理解する
• 再生核ヒルベルト空間（RKHS）とは？わかりやすく解説
• KLダイバージェンスとJSダイバージェンス — 確率分布間の距離と機械学習への応用
• 並べ替え検定（置換検定）の理論と実装をわかりやすく解説

「同じ分布か」を測るのはなぜ難しいのか

$m$ 個のサンプル $\{x_1,\dots,x_m\}$ と $n$ 個のサンプル $\{y_1,\dots,y_n\}$ が手元にあって、前者は分布 $P$ から、後者は分布 $Q$ から独立に得られたとします。知りたいのは $P=Q$ かどうか。素朴に思いつく方法をいくつか並べて、どこで詰まるかを見ておきましょう。

方法1：平均を比べる。 これが2標本t検定です。実装も解釈も簡単ですが、平均しか見ていません。$P=\mathcal{N}(0,1)$ と $Q$ が「$\pm\sqrt{3}/2$ の2つの山からなる、平均0・分散1の二峰分布」だったとき、t検定はまったく無力です。実際、本記事の後半の実験では、この二峰ケースでt検定の検出力は 0.030 と、有意水準 $\alpha=0.05$ とほぼ変わりませんでした。「差を見つけられない」のではなく「見に行っていない」のです。

方法2：ヒストグラムや密度推定で分布そのものを比べる。 カーネル密度推定(KDE) で $\hat p, \hat q$ を作り、$\int |\hat p – \hat q|$ などを計算する手です。1次元ならうまくいきますが、次元 $d$ が上がると密度推定に必要なサンプル数が指数的に増える「次元の呪い」に真正面からぶつかります。100次元の特徴ベクトルの密度を推定してから比べる、というのは現実的ではありません。

方法3：KLダイバージェンスを計算する。 $\mathrm{KL}(P\|Q)=\int p\log(p/q)$ は分布間の隔たりの標準的な尺度ですが、密度比 $p/q$ を必要とします。サンプルしか持っていない状況では、結局どこかで密度（か密度比）を推定しなければなりません。しかもサポートがずれていると $\mathrm{KL}$ は $+\infty$ に発散してしまい、「どのくらい違うか」の情報が失われます。

左の図が、モーメントを追いかける方針の限界を端的に示しています。標準正規分布と二峰分布は見た目がまるで違うのに、右の図では1次（平均）も2次（標準偏差）も3次（歪度）も小数第2位まで一致してしまい、4次（尖度）でようやく 3.01 対 1.88 という差が現れます。「では4次まで見ればいいのか」と思うかもしれませんが、4次まで一致して5次で違う分布もいくらでも作れるので、この追いかけっこには終わりがありません。一方でガウスカーネルのMMD$^2$ は、各群2000点から 0.0139 という明確な正の値を1回の計算で返します。カーネルを1つ選ぶだけで、無限個のモーメントをまとめて見に行ける のがMMDの強みです。

こうして眺めると、欲しい性質がはっきりしてきます。(i) 平均だけでなく分布の形全体を見る、(ii) 密度推定を経由しない、(iii) サンプルの内積計算だけで済むので高次元でも動く。 この3つを同時に満たすのがMMDです。次の節では、そこへ至る発想——「分布の違いを目撃してくれる関数を探す」——から始めます。

積分確率距離：分布の違いを「目撃者」に語らせる

2つの分布が違うことを、どうやって示せばよいでしょうか。裁判に例えるなら、$P$ と $Q$ の違いを目撃した 証人 を1人連れてくればよいわけです。ここでの証人とは関数 $f$ のことで、$f$ の期待値が $P$ と $Q$ で食い違えば、それが「2つは別物だ」という証拠になります。

たとえば $f(x)=x$ という証人は、平均の違いだけを証言します。$f(x)=x^2$ なら2次モーメントの違いを証言します。1人の証人では見落としが出るので、たくさんの証人を集めて、いちばん強く証言してくれた人の証言の大きさ を「違いの大きさ」と定義しましょう。これを数式にしたのが 積分確率距離（Integral Probability Metric, IPM） です。

$$ \begin{equation} \mathrm{IPM}_{\mathcal{F}}(P,Q) = \sup_{f \in \mathcal{F}} \left( \mathbb{E}_{X \sim P}[f(X)] – \mathbb{E}_{Y \sim Q}[f(Y)] \right) \end{equation} $$

$\mathcal{F}$ は「証人候補の集まり」で、これをどう選ぶかで距離の性格が決まります。

ここで $\mathcal{F}$ を広く取りすぎると、$\sup$ を計算するのが絶望的に難しくなります。全変動距離が「値としては自然だが推定できない」のはそのためです。逆に狭くしすぎると、$P\ne Q$ なのに全員が沈黙する（$\mathrm{IPM}=0$ になる）恐れが出てきます。

左の図に描かれているのが「証人候補」です。RKHSの単位球から選んだ関数はどれも滑らかで、暴れ回るような関数は入っていません。右の棒グラフを見ると、ランダムに選んだ5人の証言は $+0.082$ から $-0.020$ までばらばらで、なかには「違いを見つけられない（負の証言をする）」候補もいます。そのなかで最も強く証言する $f^\ast$ の値 0.144 が、この2分布のMMDです。証人候補を無限に集めたときの最大値 という定義が、こうして具体的な数値に落ちるわけです。

MMDの巧妙さは、$\mathcal{F}$ を 再生核ヒルベルト空間の単位球 に取った点にあります。RKHSは「滑らかな関数の集まり」であると同時に、内積が $k(x,y)$ という具体的な関数で計算できる空間です。この構造のおかげで、後で見るように $\sup$ が 手で解けてしまい、しかもカーネルを適切に選べば「誰も沈黙しない」ことも保証できます。ではその主役、平均埋め込みを導入しましょう。

平均埋め込み $\mu_P$：分布をヒルベルト空間の1点に潰す

カーネル法では、1つの点 $x$ を特徴写像 $\phi(x)=k(\cdot,x)$ によって関数（＝RKHSの1点）に写します。$x$ という数値が、$x$ を中心に盛り上がった「山の形」に化けるイメージです。ガウスカーネルなら、$\phi(x)$ は $x$ を中心とする釣鐘型の関数そのものです。

では、点ではなく 分布 を写すにはどうすればよいでしょう。素直に、分布に従って山を無限個置いて足し合わせる（＝期待値を取る）のです。KDEで「各データ点に山を置いて足す」のとまったく同じ発想を、母集団レベルでやることになります。

$$ \begin{equation} \mu_P := \mathbb{E}_{X \sim P}\left[ k(\cdot, X) \right] \in \mathcal{H} \end{equation} $$

これを 平均埋め込み（mean embedding） と呼びます。$\mu_P$ は関数であり、点 $x$ での値は $\mu_P(x) = \mathbb{E}_{X\sim P}[k(x,X)]$ です。ガウスカーネルなら、これはまさに「$P$ から無限個サンプルを取ってKDEを作った極限の形」になっています。

薄い線が「1つのデータ点に置いた山」、太い線がそれらを平均した $\hat\mu_P$、破線が理論的に計算した真の $\mu_P$ です。左の $n=8$ では太線が破線からはっきりずれており、右肩にデータが偏った分だけ山の位置もずれています。ところが $n=400$ になると太線と破線はほとんど重なります。分布 $P$ の情報が、標本数を増やすほど正確に1本の関数へ集約されていく ことが目で確認できます。この収束の速さこそ、後でMMDの推定が $O(1/\sqrt{n})$ で効く根拠になります。

$\mu_P$ がちゃんと $\mathcal{H}$ の元として存在するための条件も押さえておきましょう。Bochner積分の一般論から、

$$ \mathbb{E}_{X\sim P}\left[\sqrt{k(X,X)}\right] < \infty $$

であれば $\mu_P \in \mathcal{H}$ が保証されます。ガウスカーネルのように $k(x,x)=1$ と有界なカーネルでは、どんな $P$ に対しても自動的に成り立つので心配は要りません。

平均埋め込みが本当に強力なのは、次の等式が成り立つからです。RKHSの任意の関数 $f$ について、その期待値が $\mu_P$ との内積で書けます。

$$ \begin{equation} \langle f, \mu_P \rangle_{\mathcal{H}} = \mathbb{E}_{X\sim P}[f(X)] \end{equation} $$

導出は再生性 $\langle f, k(\cdot,x)\rangle_{\mathcal{H}} = f(x)$ を使えば1行です。内積と期待値（積分）の順序を交換すると、

$$ \langle f, \mu_P \rangle_{\mathcal{H}} = \left\langle f, \mathbb{E}_{X}[k(\cdot,X)] \right\rangle_{\mathcal{H}} = \mathbb{E}_{X}\left[ \langle f, k(\cdot,X) \rangle_{\mathcal{H}} \right] = \mathbb{E}_{X}[f(X)] $$

となります（順序交換はBochner積分の性質から正当化されます）。

この式の意味を噛みしめてください。分布 $P$ の「あらゆる証人 $f$ に対する証言 $\mathbb{E}_P[f]$」という無限個の情報が、$\mu_P$ というたった1つのベクトルにすべて畳み込まれている のです。証人 $f$ を1人指名するというのは、$\mu_P$ を $f$ の方向に射影して読み取る操作にすぎません。ならば「最も強く証言する証人」を探す作業は、単に $\mu_P – \mu_Q$ というベクトルの向きを見ればよいはずです。次節でそれを実行します。

MMDの定義と3項表現の導出

準備が整いました。証人クラスをRKHSの単位球に取ると、MMDは次のように定義されます。

$$ \begin{equation} \mathrm{MMD}(P,Q) := \sup_{\|f\|_{\mathcal{H}} \le 1} \left( \mathbb{E}_{X\sim P}[f(X)] – \mathbb{E}_{Y\sim Q}[f(Y)] \right) \end{equation} $$

まず、括弧の中身を平均埋め込みで書き換えます。前節の $\langle f,\mu_P\rangle = \mathbb{E}_P[f]$ を両方の項に適用すると、

$$ \mathbb{E}_{P}[f(X)] – \mathbb{E}_{Q}[f(Y)] = \langle f, \mu_P \rangle_{\mathcal{H}} – \langle f, \mu_Q \rangle_{\mathcal{H}} = \langle f, \mu_P – \mu_Q \rangle_{\mathcal{H}} $$

内積の線形性を使っただけです。すると求めたいのは「単位ベクトル $f$ と固定ベクトル $\mu_P-\mu_Q$ の内積の最大値」になりました。これはCauchy–Schwarzの不等式そのもので、

$$ \langle f, \mu_P – \mu_Q \rangle_{\mathcal{H}} \le \|f\|_{\mathcal{H}} \, \|\mu_P – \mu_Q\|_{\mathcal{H}} \le \|\mu_P – \mu_Q\|_{\mathcal{H}} $$

等号は $f$ が $\mu_P-\mu_Q$ と同じ向きの単位ベクトルのとき、すなわち

$$ \begin{equation} f^\ast = \frac{\mu_P – \mu_Q}{\|\mu_P – \mu_Q\|_{\mathcal{H}}} \end{equation} $$

のときに達成されます。したがって、

$$ \begin{equation} \boxed{\ \mathrm{MMD}(P,Q) = \|\mu_P – \mu_Q\|_{\mathcal{H}}\ } \end{equation} $$

「証人を全員試して最大の証言を探す」という無限次元の最適化が、たった1つのノルムの計算に化けました。 これがMMDの第一の魔法です。

最大値を与える $f^\ast$ は 証人関数（witness function） と呼ばれます。分子を書き下すと

$$ f^\ast(x) \propto \mu_P(x) – \mu_Q(x) = \mathbb{E}_{X\sim P}[k(x,X)] – \mathbb{E}_{Y\sim Q}[k(x,Y)] $$

で、ガウスカーネルなら 「$P$ のKDE」から「$Q$ のKDE」を引いた形 です。$f^\ast(x)$ が大きく正になる $x$ の付近は $P$ が過剰、大きく負になる付近は $Q$ が過剰、という具合に、どこで分布が食い違っているかを可視化してくれます。単なるスカラー値だけでなく「違いの現場」まで教えてくれるのは、実務でかなりありがたい性質です。

上段では2つの平均埋め込みが微妙にしかずれていないように見えますが、下段でその差を取ると構造がはっきり浮かび上がります。$x\approx 0$ 付近で $f^\ast$ が正に振れているのは、二峰分布の谷にあたる中央部分に $P$ だけがデータを持っているからです。逆に $x\approx \pm 1.4$ で負に沈んでいるのは、そこに $Q$ の2つの山があるからです。MMDはスカラー値を返すだけでなく、証人関数を描けば「どこが食い違っているか」まで指し示してくれます。 分布ドリフトの原因調査では、この可視化が実務上いちばん役に立つ場面も少なくありません。

二乗して3項表現に展開する

$\mathrm{MMD}$ そのものより、二乗した $\mathrm{MMD}^2$ のほうが扱いやすくなります。ノルムの二乗を内積で展開しましょう。

$$ \mathrm{MMD}^2(P,Q) = \|\mu_P – \mu_Q\|_{\mathcal{H}}^2 = \langle \mu_P – \mu_Q, \mu_P – \mu_Q \rangle_{\mathcal{H}} $$

内積の双線形性で3項に分けると、

$$ \mathrm{MMD}^2 = \langle \mu_P, \mu_P \rangle – 2\langle \mu_P, \mu_Q \rangle + \langle \mu_Q, \mu_Q \rangle $$

あとは各項を計算します。第1項について、$\mu_P$ の定義を両方に代入し、内積と期待値の順序を交換すると（$X, X’$ は $P$ から独立に取った2つのコピー）、

$$ \langle \mu_P, \mu_P \rangle_{\mathcal{H}} = \left\langle \mathbb{E}_X[k(\cdot,X)], \mathbb{E}_{X’}[k(\cdot,X’)] \right\rangle = \mathbb{E}_{X,X’}\left[ \langle k(\cdot,X), k(\cdot,X’) \rangle_{\mathcal{H}} \right] $$

ここで 再生性をもう一度使います。$\langle k(\cdot,X), k(\cdot,X’)\rangle_{\mathcal{H}} = k(X,X’)$ でしたから、

$$ \langle \mu_P, \mu_P \rangle_{\mathcal{H}} = \mathbb{E}_{X,X’\sim P}\left[ k(X,X’) \right] $$

第2項・第3項もまったく同じ手順で計算できます。まとめると、

$$ \begin{equation} \boxed{\ \mathrm{MMD}^2(P,Q) = \mathbb{E}_{X,X’}[k(X,X’)] – 2\,\mathbb{E}_{X,Y}[k(X,Y)] + \mathbb{E}_{Y,Y’}[k(Y,Y’)]\ } \end{equation} $$

ここで $X,X’ \overset{\text{iid}}{\sim} P$、$Y,Y’ \overset{\text{iid}}{\sim} Q$ です。

この式が MMDの実用上の顔 です。注目すべきは、右辺に登場するのがカーネル値の期待値だけで、密度も、RKHSも、無限次元の $\mu_P$ も一切出てこない ことです。サンプルさえあれば、カーネル行列を作って平均を取るだけで推定できます。積分確率距離の $\sup$ という抽象的な定義が、行列演算にまで降りてきました。

直感的な読み方もしておきましょう。ガウスカーネルは「2点がどれくらい近いか」を測る類似度です。第1項は「$P$ 内どうしの平均的な近さ」、第3項は「$Q$ 内どうしの平均的な近さ」、第2項は「$P$ と $Q$ の間の平均的な近さ」。$P$ と $Q$ が重なっていれば3つとも同じ値になって相殺し0になる、逆に2つの雲が離れていれば「内側の近さ」より「相互の近さ」が小さくなり、正の値が残ります。「クラスタ内類似度からクラスタ間類似度を引く」という、非常に素直な量なのです。

3項表現を実際のカーネル行列で見ると、左の図のように行列が2×2のブロックに分かれます。左上と右下が「同じ分布の中どうしの類似度」、右上と左下が「PとQをまたいだ類似度」です。右の棒グラフが示すとおり、この例では $K_{XX}$ の平均 0.4593 と $K_{YY}$ の平均 0.4130 に対して、$K_{XY}$ の平均は 0.4234 とわずかに小さい。この「内側どうしの近さ」と「またいだ近さ」のわずかな差を2倍して引くだけで、MMD² = 0.0254 が得られます。 難解な $\sup$ の定義が、行列を3つのブロックに分けて平均するという操作にまで降りてきたことになります。

では、この量が実際にどんな値になるのか。次はガウス分布どうしで閉じた形を計算してみます。

具体例：ガウス分布どうしのMMDを手計算する

ガウスカーネルと1次元ガウス分布の組み合わせなら、$\mathrm{MMD}^2$ は完全に閉じた形で書けます。理論値と実測値を突き合わせられるので、理解の確認にちょうどよい題材です。

必要な部品は次の1つの積分だけです。$Z \sim \mathcal{N}(m, v)$ のとき、

$$ \mathbb{E}\left[ \exp\left(-\frac{Z^2}{2\sigma^2}\right) \right] = \frac{\sigma}{\sqrt{\sigma^2+v}} \exp\left( -\frac{m^2}{2(\sigma^2+v)} \right) $$

を導いておきます。定義通りに積分を書くと

$$ I = \int \frac{1}{\sqrt{2\pi v}} \exp\left(-\frac{(z-m)^2}{2v}\right)\exp\left(-\frac{z^2}{2\sigma^2}\right) dz $$

指数部をまとめ、$z$ について平方完成します。$A := \frac{1}{v}+\frac{1}{\sigma^2} = \frac{\sigma^2+v}{v\sigma^2}$ と置くと、指数部は

$$ -\frac{1}{2}\left[ A z^2 – \frac{2mz}{v} + \frac{m^2}{v} \right] = -\frac{A}{2}\left( z – \frac{m}{vA} \right)^2 + \frac{m^2}{2v^2 A} – \frac{m^2}{2v} $$

定数項をまとめます。 $\frac{m^2}{2v^2A} = \frac{m^2 \sigma^2}{2v(\sigma^2+v)}$ を使うと、

$$ \frac{m^2}{2v^2A} – \frac{m^2}{2v} = \frac{m^2}{2v}\left( \frac{\sigma^2}{\sigma^2+v} – 1 \right) = -\frac{m^2}{2(\sigma^2+v)} $$

残ったガウス積分は $\int \exp(-\frac{A}{2}(z-c)^2)dz = \sqrt{2\pi/A}$ なので、

$$ I = \frac{1}{\sqrt{2\pi v}}\sqrt{\frac{2\pi}{A}}\exp\left(-\frac{m^2}{2(\sigma^2+v)}\right) = \frac{1}{\sqrt{vA}}\exp(\cdots) = \frac{\sigma}{\sqrt{\sigma^2+v}}\exp\left(-\frac{m^2}{2(\sigma^2+v)}\right) $$

これで部品が揃いました。$P=\mathcal{N}(\mu_1,s_1^2)$、$Q=\mathcal{N}(\mu_2,s_2^2)$、カーネルの帯域幅を $\sigma$ とします。$X,X’\sim P$ が独立なら $X-X’\sim\mathcal{N}(0,2s_1^2)$ なので、上の公式に $m=0, v=2s_1^2$ を代入して

$$ \mathbb{E}[k(X,X’)] = \frac{\sigma}{\sqrt{\sigma^2+2s_1^2}} $$

$X\sim P, Y\sim Q$ なら $X-Y\sim\mathcal{N}(\mu_1-\mu_2, s_1^2+s_2^2)$ なので、$m=\mu_1-\mu_2, v=s_1^2+s_2^2$ を代入して

$$ \mathbb{E}[k(X,Y)] = \frac{\sigma}{\sqrt{\sigma^2+s_1^2+s_2^2}}\exp\left(-\frac{(\mu_1-\mu_2)^2}{2(\sigma^2+s_1^2+s_2^2)}\right) $$

3項表現に入れれば完成です。

$$ \begin{equation} \mathrm{MMD}^2 = \frac{\sigma}{\sqrt{\sigma^2+2s_1^2}} + \frac{\sigma}{\sqrt{\sigma^2+2s_2^2}} – \frac{2\sigma}{\sqrt{\sigma^2+s_1^2+s_2^2}}\exp\left(-\frac{(\mu_1-\mu_2)^2}{2(\sigma^2+s_1^2+s_2^2)}\right) \end{equation} $$

数値を入れてみましょう。すべて $\sigma=1$ とします。

ケース1（平均差） $P=\mathcal{N}(0,1)$、$Q=\mathcal{N}(0.5,1)$。 第1項・第2項はともに $1/\sqrt{3}=0.57735$。第3項の中身は $\frac{1}{\sqrt{3}}\exp(-0.25/6) = 0.57735 \times 0.95919 = 0.55379$。よって

$$ \mathrm{MMD}^2 = 0.57735 + 0.57735 – 2\times 0.55379 = 0.04712 $$

ケース2（分散差） $P=\mathcal{N}(0,1)$、$Q=\mathcal{N}(0,2^2)$。 第1項 $=1/\sqrt{3}=0.57735$、第2項 $=1/\sqrt{1+8}=0.33333$、第3項 $=2/\sqrt{1+1+4}=2/\sqrt{6}=0.81650$。よって

$$ \mathrm{MMD}^2 = 0.57735 + 0.33333 – 0.81650 = 0.09419 $$

ケース3（同一分布） $\mu_1=\mu_2$、$s_1=s_2$ なら3項が打ち消し合って厳密に $0$ になります。式の上でも $\frac{2\sigma}{\sqrt{\sigma^2+2s_1^2}}$ が2回引かれる形なので、確かに $0$ です。

実際に標本から推定してみると、各群4000点×100回の平均で ケース1が 0.04665、ケース2が 0.09283 と測れました。理論値 0.04712 / 0.09419 に対して、いずれも標準誤差（0.0005 程度）の数倍以内に収まっており、閉じた式が正しいことを確認できます。

左の曲線は、同じ2分布に対して帯域幅 $\sigma$ を振ったときのMMD$^2$ の理論値です。$\sigma$ が小さすぎても大きすぎても値が0に潰れ、$\sigma\approx 1$ 付近で最大になる山なりの形をしています。MMDの値は分布のペアだけでなくカーネルにも依存するので、$\sigma$ の違う値どうしを比較しても意味がありません。 右の棒グラフでは $\sigma=1$ に固定した理論値と実測値が誤差棒の範囲で一致しており、同時に「平均差 0.5 のケース1（0.047）より、平均が完全に一致しているケース2（0.094）のほうが値が大きい」という逆転も確認できます。

この具体例には、もう1つ重要な教訓があります。ケース2では平均がまったく同じ（どちらも0）なのに $\mathrm{MMD}^2$ は 0.094 と、平均が 0.5 ずれているケース1の 0.047 より大きい のです。つまりガウスカーネルのMMDは、平均以外の違いをきちんと拾っています。一方、線形カーネル $k(x,y)=xy$ を使うと $\mu_P = \mathbb{E}[X]$ そのものになるので $\mathrm{MMD}^2 = (\mu_1-\mu_2)^2$、ケース2では $0$ になってしまいます。カーネル選びが本質的だということです。

そこで次の疑問が生まれます。どんなカーネルを使えば「$\mathrm{MMD}=0$ なら本当に $P=Q$」と言い切れるのか？ ここがMMD理論の核心です。

特性カーネル：MMD $=0 \iff P=Q$ はいつ成り立つか

$\mathrm{MMD}(P,Q) = \|\mu_P – \mu_Q\|_{\mathcal{H}}$ なので、$\mathrm{MMD}=0$ は $\mu_P = \mu_Q$ と同じことです。したがって知りたいのは、「分布 $P$ を $\mu_P$ に写す写像が単射か」 という一点に尽きます。単射なら $\mu_P=\mu_Q \Rightarrow P=Q$、単射でなければ「違う分布なのに同じ点に潰れる」ペアが存在してしまいます。

$$ \begin{equation} P \mapsto \mu_P \ \text{が単射} \quad \Longleftrightarrow \quad k \ \text{は\textbf{特性カーネル}（characteristic kernel）} \end{equation} $$

単射でないカーネルの例は、すでに見ました。線形カーネルでは $\mu_P$ が平均だけなので、平均が同じ分布はすべて同じ点に潰れます。2次多項式カーネル $k(x,y)=(x^\top y + 1)^2$ なら1次と2次のモーメントまでは区別できますが、それ以上は潰れます。多項式カーネルは「有限個のモーメントしか見ない証人団」なのです。

では、なぜガウスカーネルは大丈夫なのでしょうか。並進不変カーネル $k(x,y)=\psi(x-y)$ の場合、Bochnerの定理 を使うと美しく証明できます。定理によれば、連続な正定値関数 $\psi$ は有限非負測度 $\Lambda$（スペクトル測度）のFourier変換として書けます。

$$ \psi(\delta) = \int_{\mathbb{R}^d} e^{-i\omega^\top \delta}\, d\Lambda(\omega) $$

これを3項表現に代入していきます。$P$ の特性関数を $\varphi_P(\omega)=\mathbb{E}_{X\sim P}[e^{-i\omega^\top X}]$ と書くと、第1項は

$$ \mathbb{E}_{X,X’}[k(X,X’)] = \int \mathbb{E}_{X,X’}\left[e^{-i\omega^\top (X-X’)}\right] d\Lambda(\omega) $$

$X$ と $X’$ が独立なので期待値が積に分かれ、$\mathbb{E}[e^{i\omega^\top X’}] = \overline{\varphi_P(\omega)}$ ですから、

$$ \mathbb{E}_{X,X’}[k(X,X’)] = \int \varphi_P(\omega)\overline{\varphi_P(\omega)}\, d\Lambda(\omega) = \int |\varphi_P(\omega)|^2 d\Lambda(\omega) $$

まったく同様に、交差項は $\int \varphi_P \overline{\varphi_Q}\, d\Lambda$、第3項は $\int |\varphi_Q|^2 d\Lambda$ になります。3つをまとめると、複素数の $|a-b|^2 = |a|^2 – 2\mathrm{Re}(a\bar b) + |b|^2$ という展開がそのまま現れます。

$$ \begin{equation} \mathrm{MMD}^2(P,Q) = \int_{\mathbb{R}^d} \left| \varphi_P(\omega) – \varphi_Q(\omega) \right|^2 d\Lambda(\omega) \end{equation} $$

これは非常に示唆的な式です。MMDとは、2つの分布の特性関数の差を、カーネルのスペクトル測度で重み付けして測った $L^2$ 距離だった のです。

ここから特性カーネルの条件がただちに読み取れます。$\mathrm{MMD}=0$ なら被積分関数は $\Lambda$-a.e. で $0$、つまり $\Lambda$ の台の上で $\varphi_P=\varphi_Q$。もし $\Lambda$ の台が $\mathbb{R}^d$ 全体 なら、特性関数が連続であることと合わせて $\varphi_P \equiv \varphi_Q$ が全域で成り立ち、特性関数の一意性定理から $P=Q$ が従います。

ガウスカーネル $k(x,y)=\exp(-\|x-y\|^2/(2\sigma^2))$ のスペクトル測度は $\mathcal{N}(0,\sigma^{-2}I)$ で、密度はどこでも正、つまり台は $\mathbb{R}^d$ 全体です。だから ガウスカーネルは特性カーネル です。ラプラスカーネル $\exp(-\|x-y\|_1/\sigma)$ のスペクトル測度もCauchy型で台は全域なので、これも特性カーネルです。

逆に、$\Lambda$ の台に「穴」があると、その周波数帯だけで違う2つの分布を区別できなくなります。sinc型のカーネル（帯域制限）が特性でないのはこのためです。

この3枚が、Bochner表現の意味をそのまま絵にしたものです。左を見ると、2つの分布の特性関数は $\omega\approx 0$ 付近ではどちらも1に張り付いていて区別できず、$\omega\approx \pm 3$ で最も大きく食い違います。中央の図では、その「差の2乗」（黒）と、帯域幅ごとのスペクトル測度（色付き）を重ねました。$\sigma=3$ のスペクトル測度は $\omega=0$ の狭い範囲に集中しているため、差が大きい $\omega\approx\pm 3$ の情報をまったく拾えません。 右の図はその帰結で、Bochner表現の積分値（黒線）と標本から計算したMMD$^2$（緑点）は $\sigma=0.2$ で 0.0212 対 0.0205、$\sigma=1.5$ で 0.00136 対 0.00145 と一致し、$\sigma$ を大きくすると値が0へ潰れていきます。特性カーネルであっても「重みをどこに置くか」で見える違いが変わる、という後の帯域幅の議論の伏線がここにあります。

最後の行について補足すると、統計学で古くから使われるエネルギー距離 $2\mathbb{E}\|X-Y\| – \mathbb{E}\|X-X’\| – \mathbb{E}\|Y-Y’\|$ は、距離から作られる特定のカーネルに対する $\mathrm{MMD}^2$ と一致することが知られています（Sejdinovic et al., 2013）。エネルギー距離はMMDの特別な場合だった わけで、この統一的な視点はMMDの理論的な射程の広さを物語っています。

ここまでは母集団レベル（真の期待値）の話でした。実際に手元にあるのは有限個のサンプルです。次は推定の話に移ります。

有限標本からの推定：U統計量とV統計量

3項表現の各期待値を、標本平均で置き換えればよさそうです。ところが「対角成分をどう扱うか」で2つの流儀が生まれ、性質が大きく変わります。

V統計量（バイアスあり）

いちばん素直なのは、経験分布の平均埋め込み $\hat\mu_P = \frac{1}{m}\sum_i k(\cdot,x_i)$ をそのまま使う方法です。

$$ \begin{equation} \widehat{\mathrm{MMD}}^2_b = \|\hat\mu_P – \hat\mu_Q\|^2_{\mathcal{H}} = \frac{1}{m^2}\sum_{i,j}^{m} k(x_i,x_j) + \frac{1}{n^2}\sum_{i,j}^{n} k(y_i,y_j) – \frac{2}{mn}\sum_{i,j} k(x_i,y_j) \end{equation} $$

これを V統計量 と呼びます。$i=j$ の対角項（$k(x_i,x_i)$）も含めて平均を取っているのがポイントです。ノルムの二乗そのものなので、必ず $\ge 0$ になります。 この「絶対に負にならない」という性質は、MMDを損失関数として最小化する場面（ドメイン適応など）で重要になります。

バイアスを計算してみる

V統計量にはバイアスがあります。$m=n$ として第1項の期待値を計算しましょう。$m^2$ 個の項のうち、$i=j$ の対角が $m$ 個、$i\ne j$ が $m(m-1)$ 個です。対角では同じ点どうしのカーネル値 $k(X,X)$、非対角では独立な2点のカーネル値 $k(X,X’)$ が現れる ことに注意して、

$$ \mathbb{E}\left[\frac{1}{m^2}\sum_{i,j} k(x_i,x_j)\right] = \frac{m}{m^2}\mathbb{E}[k(X,X)] + \frac{m(m-1)}{m^2}\mathbb{E}[k(X,X’)] $$

$= \frac{1}{m}\mathbb{E}[k(X,X)] + \left(1-\frac{1}{m}\right)\mathbb{E}[k(X,X’)]$ と整理できます。真に欲しかったのは $\mathbb{E}[k(X,X’)]$ ですから、差分がバイアスです。第3項（交差項）は $i,j$ が別々の標本なので独立で不偏。よって全体では

$$ \begin{equation} \mathbb{E}\left[\widehat{\mathrm{MMD}}^2_b\right] = \mathrm{MMD}^2 + \frac{1}{m}\Big(\mathbb{E}[k(X,X)] – \mathbb{E}[k(X,X’)]\Big) + \frac{1}{n}\Big(\mathbb{E}[k(Y,Y)] – \mathbb{E}[k(Y,Y’)]\Big) \end{equation} $$

ガウスカーネルなら $k(x,x)=1$、また $\mathbb{E}[k(X,X’)] < 1$ なので、バイアスは必ず正、つまりV統計量は分布の違いを常に過大評価します。しかも $O(1/m)$ で減衰します。

具体的に見ておきましょう。$P=Q=\mathcal{N}(0,1)$、$\sigma=1$ のとき、前節の公式から $\mathbb{E}[k(X,X’)] = 1/\sqrt{1+2}=1/\sqrt{3}=0.5774$。$m=n$ なら

$$ \text{バイアス} = \frac{2}{m}\left(1 – \frac{1}{\sqrt{3}}\right) = \frac{0.8453}{m} $$

$m=100$ なら 0.00845。真値は 0 なのに、V統計量は平均 0.008 程度の値を返します。後の実験でも $m=100$ で 0.00824 と、この理論値をほぼ言い当てました。検定に使うには致命的です。 真値0でも正の値が出るので、素朴に「0より大きいから違う分布だ」と判定できません。

U統計量（不偏）

そこで、対角を除く という修正を入れます。

$$ \begin{equation} \widehat{\mathrm{MMD}}^2_u = \frac{1}{m(m-1)}\sum_{i \ne j}^{m} k(x_i,x_j) + \frac{1}{n(n-1)}\sum_{i \ne j}^{n} k(y_i,y_j) – \frac{2}{mn}\sum_{i,j} k(x_i,y_j) \end{equation} $$

対角を除いたので、残るのはすべて「独立な2点のペア」です。したがって各項が不偏となり、

$$ \mathbb{E}\left[\widehat{\mathrm{MMD}}^2_u\right] = \mathrm{MMD}^2(P,Q) $$

が厳密に成り立ちます。これが U統計量 です。ただし引き換えに、$\widehat{\mathrm{MMD}}^2_u$ は負の値を取り得ます。$P=Q$ のときは真値0のまわりで正負にばらつくので、これは不偏性の当然の帰結です。「MMDの二乗が負になった」と慌てないでください。

使い分けは明快です。検定にはU統計量、最適化の目的関数にはV統計量。 前者は「0との比較」が問題になるので不偏性が要り、後者は「勾配を流して小さくする」だけなので $\ge 0$ で滑らかなほうが都合がよいのです。

計算量と線形時間推定量

どちらも $m+n$ 点のカーネル行列を作るので $O((m+n)^2)$ の時間とメモリが必要です。$m=10^5$ になると $10^{10}$ 要素で現実的ではありません。そこで、サンプルを2個ずつのペアに区切って平均を取る 線形時間推定量 が用意されています。$z_i=(x_i,y_i)$ とし、

$$ h(z, z’) := k(x,x’) + k(y,y’) – k(x,y’) – k(x’,y) $$

と置くと、$\mathbb{E}[h] = \mathrm{MMD}^2$ が成り立ちます（3項表現をそのまま2点で近似した形です）。これを重ならないペアだけで平均すれば、

$$ \widehat{\mathrm{MMD}}^2_\ell = \frac{2}{m}\sum_{i=1}^{\lfloor m/2 \rfloor} h(z_{2i-1}, z_{2i}) $$

が得られます。$O(m)$ 時間・$O(1)$ メモリで不偏。ただしペアを大量に捨てているので分散は大きく、同じ検出力を得るにはサンプル数が余分に要ります。「時間を取るか、サンプル効率を取るか」のトレードオフ です。

推定量が決まったところで、次はカーネルの帯域幅 $\sigma$ をどう決めるかという実務上の最大の悩みに進みます。

帯域幅の選び方：メディアンヒューリスティック

ガウスカーネルの帯域幅 $\sigma$ は、MMDの「解像度」を決めるツマミです。極端な2つの場合を考えると、その役割がはっきりします。

$\sigma \to \infty$（幅が広すぎる）：$\exp(-\|x-y\|^2/(2\sigma^2)) \approx 1 – \frac{\|x-y\|^2}{2\sigma^2}$ とTaylor展開でき、定数項は3項表現で打ち消し合います。残るのは $\|x-y\|^2$ の期待値の組み合わせで、これを整理すると $\mathrm{MMD}^2 \approx \frac{1}{\sigma^2}\|\mathbb{E}[X]-\mathbb{E}[Y]\|^2$ に比例する形になります。つまり広すぎるガウスカーネルは、線形カーネルと同じく平均しか見なくなります。 特性カーネルであっても、有限標本では平均以外の違いがノイズに埋もれてしまうのです。

$\sigma \to 0$（幅が狭すぎる）：$k(x_i,x_j)\to 0$（$i\ne j$）となり、どの2点も「無関係」になります。U統計量はほぼ全項が0となり、統計量が真値を拾えません。V統計量は対角しか残らず、$\mathrm{MMD}^2 \to 1+1-0=2$ という、分布によらない定数に張り付きます。どちらも「見分けがつかない」方向です。

つまり $\sigma$ は「データ点どうしの典型的な距離」と同じくらいのオーダーに取るのが自然、ということになります。これを実装したのが メディアンヒューリスティック です。

$$ \begin{equation} \sigma^2 = \frac{1}{2}\,\mathrm{median}\left\{ \|z_i – z_j\|^2 : i < j \right\} \end{equation} $$

ここで $\{z_i\}$ は $X$ と $Y$ を混ぜた全サンプルです。$1/2$ を掛けているのは、この $\sigma$ を使うとカーネルの指数部が典型的なペアでちょうど $-1/2$ 程度になり、$k$ が $e^{-1/2}\approx 0.6$ という「飽和も消滅もしない」値になるからです。実装は3行で書け、チューニング不要なので、実務ではまずこれを使います。

ただし メディアンヒューリスティックは最適ではありません。本記事の実験では、二峰分布を見分けるタスクで、メディアンヒューリスティック（$\sigma\approx 0.71$）の検出力が 0.355 だったのに対し、$\sigma=0.2$ では 0.570 と 1.6 倍になりました。二峰の山の間隔は $\sqrt{3}\approx 1.73$ ですが、山そのものの幅は 0.5。「見分けたい構造のスケール」に $\sigma$ を合わせたほうが検出力が上がる わけです。逆に $\sigma=1.5$ 以上では検出力が 0.055 まで落ち、有意水準 0.05 とほぼ変わらなくなります。前述の「広すぎると平均しか見えない」が現実に効いています。

実務的な対処は2つあります。1つは マルチカーネルMMD：複数の $\sigma$（たとえばメディアンの $\{1/4,1/2,1,2,4\}$ 倍）のカーネルを平均して使う方法で、どのスケールの違いも取りこぼしにくくなります。ドメイン適応の代表的手法DAN（Deep Adaptation Networks）はこれを採用しています。もう1つは カーネル自体を学習する 方法で、データを2分割して片方で検出力を最大化する $\sigma$（あるいはニューラルネットで表現した深層カーネル）を選び、もう片方で検定します。

帯域幅が決まれば、あとは「この $\widehat{\mathrm{MMD}}^2_u$ の値は偶然にしては大きいのか」を判定するだけです。次節でその手続きを固めます。

二標本検定：帰無分布と置換検定

検定の形は標準的です。帰無仮説 $H_0: P=Q$、対立仮説 $H_1: P\ne Q$ とし、$\widehat{\mathrm{MMD}}^2_u$ が大きければ $H_0$ を棄却します。問題は 「どのくらい大きければ棄却してよいか」の閾値 です。

理論的には、Gretton et al. (2012) が漸近分布を与えています。$H_1$ のもとでは中心極限定理が効いて

$$ \sqrt{m}\left( \widehat{\mathrm{MMD}}^2_u – \mathrm{MMD}^2 \right) \xrightarrow{d} \mathcal{N}(0, 4\sigma_u^2) $$

と正規分布に収束します。ところが $H_0$（$P=Q$）のもとでは、$\mathrm{MMD}^2$ が最小値0を取るため1次の変動が消えてしまい、退化したU統計量 になります。このとき $m$ 倍したスケールで、正規分布ではなく無限個のカイ二乗の重み付き和に収束します。

$$ \begin{equation} m\, \widehat{\mathrm{MMD}}^2_u \xrightarrow{d} \sum_{l=1}^{\infty} \lambda_l \left( z_l^2 – 1 \right), \qquad z_l \overset{\text{iid}}{\sim} \mathcal{N}(0,1) \end{equation} $$

ここで $\lambda_l$ は、中心化カーネル

$$ \tilde k(x,x’) = k(x,x’) – \mathbb{E}_{x}[k(x,x’)] – \mathbb{E}_{x’}[k(x,x’)] + \mathbb{E}_{x,x’}[k(x,x’)] $$

が定める積分作用素の固有値です。この $\lambda_l$ は真の分布 $P$ に依存するので、表を引いて閾値を得る、という使い方ができません。 ここがt検定との決定的な違いです。

そこで実務では 置換検定（permutation test） を使います。発想は単純で、「$H_0$ が正しければ、$X$ と $Y$ のラベルはただの名札にすぎず、貼り替えても統計量の分布は変わらないはずだ」という交換可能性を利用します。手順は次の通りです。

1. 観測値から $\widehat{\mathrm{MMD}}^2_u$ を計算し、$T_{\text{obs}}$ とする。
2. $X$ と $Y$ をすべて混ぜて $N=m+n$ 個のプールを作る。
3. プールをランダムに $m$ 個と $n$ 個に分け直し、$\widehat{\mathrm{MMD}}^2_u$ を計算する。これを $B$ 回繰り返して $T_{(1)},\dots,T_{(B)}$ を得る。
4. p値を次で計算する。

$$ \begin{equation} \hat p = \frac{1 + \#\{ b : T_{(b)} \ge T_{\text{obs}} \}}{B+1} \end{equation} $$

分子・分母の $+1$ は「観測値そのものも1つの並べ替え結果として数える」ための補正で、これがあると $H_0$ のもとで $\hat p$ が正しく一様分布に近づき、第一種の過誤が名目水準を超えません。

実装上のコツを1つ。カーネル行列 $K$ は最初に1回だけ作れば十分です。 並べ替えは行と列を同じ順序で入れ替えるだけなので、K[np.ix_(idx, idx)] とインデックス操作すれば、カーネル関数の再評価が要りません。$B=200$ 回の並べ替えでも、カーネル計算のコストは1回分で済みます。

理論も手続きも揃いました。ここからは実際に動かして、これまでの主張を1つずつ数値で確かめていきます。

Pythonでの実装

準備：カーネルとMMD推定量

まず、ガウスカーネル・線形カーネル・U統計量・V統計量・メディアンヒューリスティックを実装します。

import numpy as np

def gaussian_kernel(A, B, sigma):
    """ガウスカーネル行列 K[i,j] = exp(-||a_i - b_j||^2 / (2 sigma^2))"""
    d2 = ((A[:, None, :] - B[None, :, :]) ** 2).sum(axis=-1)
    return np.exp(-d2 / (2.0 * sigma ** 2))

def linear_kernel(A, B):
    return A @ B.T

def mmd2_unbiased(X, Y, kfun):
    """U統計量（不偏）: 対角を除いて平均を取る"""
    m, n = len(X), len(Y)
    Kxx, Kyy, Kxy = kfun(X, X), kfun(Y, Y), kfun(X, Y)
    sxx = (Kxx.sum() - np.trace(Kxx)) / (m * (m - 1))
    syy = (Kyy.sum() - np.trace(Kyy)) / (n * (n - 1))
    return sxx + syy - 2.0 * Kxy.mean()

def mmd2_biased(X, Y, kfun):
    """V統計量（バイアスあり）: 対角も含めた素直な平均"""
    return kfun(X, X).mean() + kfun(Y, Y).mean() - 2.0 * kfun(X, Y).mean()

def median_sigma(Z):
    """メディアンヒューリスティック: sigma^2 = median(||z_i - z_j||^2) / 2"""
    d2 = ((Z[:, None, :] - Z[None, :, :]) ** 2).sum(axis=-1)
    iu = np.triu_indices(len(Z), k=1)
    return np.sqrt(np.median(d2[iu]) / 2.0)

mmd2_unbiased で np.trace を引いているのが、理論で見た「対角を除く」操作そのものです。mmd2_biased はそれをせず素直に mean() を取るだけで、この1行の違いがバイアスの有無を生みます。median_sigma では上三角のみ（k=1）を使い、自分自身との距離0を中央値の計算から除いています。ここを忘れると中央値が小さくなりすぎるので注意してください。

実験1：平均差・分散差・形状差の3ケース

理論で予告した「線形カーネルは平均しか見ない」を実験で確認します。$P=\mathcal{N}(0,1)$ に対し、平均だけずらしたもの、分散だけ広げたもの、平均も分散も一致させたまま二峰にしたものの3種類を用意します。二峰分布は $\pm\sqrt{3}/2$ に中心を置き、標準偏差 0.5 のノイズを乗せることで、平均 0・分散 $3/4+1/4=1$ とぴったり揃えました。

import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

rng = np.random.default_rng(0)
n, a = 300, np.sqrt(3) / 2
P = rng.normal(0.0, 1.0, (n, 1))
cases = {
    "同一分布": rng.normal(0.0, 1.0, (n, 1)),
    "平均差": rng.normal(0.5, 1.0, (n, 1)),
    "分散差": rng.normal(0.0, 2.0, (n, 1)),
    "形状差(二峰)": (np.where(rng.random(n) < 0.5, -a, a) + rng.normal(0, 0.5, n))[:, None],
}

names, vg, vl = [], [], []
for name, Q in cases.items():
    s = median_sigma(np.vstack([P, Q]))
    mg = mmd2_unbiased(P, Q, lambda A, B: gaussian_kernel(A, B, s))
    ml = mmd2_unbiased(P, Q, linear_kernel)
    names.append(name); vg.append(mg); vl.append(ml)
    print(f"{name:14s} 平均差={float(P.mean()-Q.mean()):+.4f} "
          f"標準偏差差={float(P.std()-Q.std()):+.4f} "
          f"MMD2(ガウス)={mg:.5f} MMD2(線形)={ml:.5f}")

出力は次のようになります。

同一分布           平均差=-0.0262 標準偏差差=+0.0426 MMD2(ガウス)=-0.00175 MMD2(線形)=-0.00596
平均差             平均差=-0.5088 標準偏差差=+0.0914 MMD2(ガウス)=0.05191 MMD2(線形)=0.25257
分散差             平均差=+0.0755 標準偏差差=-0.9603 MMD2(ガウス)=0.08336 MMD2(線形)=-0.01085
形状差(二峰)       平均差=-0.0104 標準偏差差=-0.0658 MMD2(ガウス)=0.02336 MMD2(線形)=-0.00728

この表がMMDの本質を凝縮しています。まず 線形カーネル（右端）は、平均差ケースでのみ 0.253 という大きな値を返し、分散差・形状差では $-0.011$、$-0.007$ とほぼ0 です。理論通り「平均しか見ていない」ことが数値で確認できました。一方 ガウスカーネル（左隣）は3ケースすべてで正の値（0.052 / 0.083 / 0.023）を返し、同一分布ケースの $-0.0018$ とは明確に桁が違います。とくに形状差ケースは平均も分散も一致しているので、t検定でも分散のF検定でも検出できませんが、ガウスカーネルのMMDは同一分布の約13倍の値を出して違いを捉えています。

左の棒グラフを見ると、灰色（線形カーネル）が突出しているのは平均差のケースだけで、分散差では $-0.011$、形状差では $-0.007$ と、同一分布のケース（$-0.006$）と区別がつきません。線形カーネルは「平均という1人の証人しか呼んでいない」ため、それ以外の違いに対しては構造的に無力です。 対する緑（ガウスカーネル）は3つのケースすべてで正の値を返し、右の図に示した「平均も分散も揃った二峰分布」でさえ 0.023 と、同一分布の $-0.002$ から明確に離れています。カーネルの選択が検出できる違いの種類を決める、という理論上の主張が数値で裏づけられました。

また、ガウスカーネルの値は理論値ともよく合っています。平均差ケースの理論値は $\sigma=1$ のとき 0.047 でしたが、ここではメディアンヒューリスティックで $\sigma\approx 0.66$ が選ばれているため 0.052 とやや大きめ。帯域幅が違えば値も変わるので、MMDの絶対値どうしを異なる $\sigma$ で比較してはいけない ことも読み取れます。

実験2：U統計量とV統計量のバイアス

次に、$P=Q$（真値は厳密に0）という条件で2つの推定量の平均を比べ、理論バイアス $\frac{2}{m}(1-1/\sqrt{3})$ と照合します。

import numpy as np

ms = [20, 50, 100, 200, 400]
kf = lambda A, B: gaussian_kernel(A, B, 1.0)
for m in ms:
    r = np.random.default_rng(1)
    us, vs = [], []
    for _ in range(400):  # 400回のモンテカルロ
        X = r.normal(0, 1, (m, 1)); Y = r.normal(0, 1, (m, 1))
        us.append(mmd2_unbiased(X, Y, kf))
        vs.append(mmd2_biased(X, Y, kf))
    theory = 2 / m * (1 - 1 / np.sqrt(3))
    print(f"m={m:4d}  E[U]={np.mean(us):+.5f}  E[V]={np.mean(vs):+.5f}  理論バイアス={theory:.5f}")

m=  20  E[U]=-0.00001  E[V]=+0.04211  理論バイアス=0.04226
m=  50  E[U]=-0.00151  E[V]=+0.01531  理論バイアス=0.01691
m= 100  E[U]=-0.00019  E[V]=+0.00824  理論バイアス=0.00845
m= 200  E[U]=-0.00026  E[V]=+0.00396  理論バイアス=0.00423
m= 400  E[U]=+0.00002  E[V]=+0.00213  理論バイアス=0.00211

左の図では、赤い折れ線（V統計量）が黒い破線（理論バイアス $0.845/m$）にぴたりと重なり、青い折れ線（U統計量）は $m$ によらず0の線上に張りついています。真値は厳密に0なので、赤と黒が重なることは「V統計量のずれが偶然ではなく理論どおりのバイアスである」ことの証明 になっています。右の棒グラフは同じ内容を数値で並べたもので、$m=400$ では 0.0021 対 0.0021 とほぼ完全に一致します。$m$ を4倍にするとバイアスが約4分の1になる、$O(1/m)$ の減衰もはっきり読み取れます。

3列を見比べると、理論が数値でぴたりと再現されているのがわかります。U統計量の平均はどの $m$ でも $\pm 0.0015$ 以内で0を挟んでおり、不偏性が確認できます（$m=50$ の $-0.00151$ はモンテカルロ誤差）。対して V統計量は常に正で、しかもその値は理論式 $\frac{2}{m}(1-1/\sqrt{3})=0.845/m$ とよく一致 しています（$m=100$ で 0.00824 対 0.00845、相対誤差 2.5%）。

実務的な含意は重大です。$m=20$ のときV統計量は 0.042 を返しますが、これは実験1の「形状差(二峰)」で得られた 0.023 より大きな値です。サンプルが少ないと、同じ分布でも「二峰との違い」より大きな値が出てしまう。 V統計量をそのまま検定に使ってはいけない理由がここにあります。

実験3：置換検定でp値を出す

理論的な閾値が使えない以上、並べ替えで帰無分布を作ります。カーネル行列を1回だけ計算して使い回すのがポイントです。

import numpy as np
from scipy import stats

def mmd2u_from_K(K, m):
    """既に計算済みのカーネル行列からU統計量を計算"""
    n = K.shape[0] - m
    Kxx, Kyy, Kxy = K[:m, :m], K[m:, m:], K[:m, m:]
    return ((Kxx.sum() - np.trace(Kxx)) / (m * (m - 1))
            + (Kyy.sum() - np.trace(Kyy)) / (n * (n - 1)) - 2.0 * Kxy.mean())

def mmd_permutation_test(X, Y, B=200, rng=None, sigma=None):
    Z = np.vstack([X, Y]); m, N = len(X), len(X) + len(Y)
    s = median_sigma(Z) if sigma is None else sigma
    K = gaussian_kernel(Z, Z, s)          # カーネル行列は1回だけ計算
    obs = mmd2u_from_K(K, m)
    null = np.empty(B)
    for b in range(B):
        idx = rng.permutation(N)          # ラベルを貼り替える
        null[b] = mmd2u_from_K(K[np.ix_(idx, idx)], m)
    p = (1 + np.sum(null >= obs)) / (B + 1)   # +1補正付きのp値
    return obs, p, null, s

この関数を、同一分布・平均差 0.4・二峰（平均分散一致）の3シナリオに適用し、t検定と並べます。

import numpy as np
from scipy import stats

rng = np.random.default_rng(11)
n, a = 100, np.sqrt(3) / 2
scen = {
    "同一分布":          (rng.normal(0, 1, (n, 1)), rng.normal(0, 1, (n, 1))),
    "平均差0.4":         (rng.normal(0, 1, (n, 1)), rng.normal(0.4, 1, (n, 1))),
    "二峰(平均分散一致)": (rng.normal(0, 1, (n, 1)),
                          (np.where(rng.random(n) < .5, -a, a) + rng.normal(0, .5, n))[:, None]),
}
for name, (X, Y) in scen.items():
    obs, p, null, s = mmd_permutation_test(X, Y, B=200, rng=rng)
    tp = stats.ttest_ind(X.ravel(), Y.ravel()).pvalue
    print(f"{name:20s} sigma={s:.3f}  MMD^2_u={obs:+.5f}  p(MMD)={p:.4f}  p(t検定)={tp:.4f}")

同一分布                 sigma=0.635  MMD^2_u=-0.00159  p(MMD)=0.4677  p(t検定)=0.8993
平均差0.4                sigma=0.659  MMD^2_u=+0.03422  p(MMD)=0.0050  p(t検定)=0.0008
二峰(平均分散一致)        sigma=0.745  MMD^2_u=+0.01394  p(MMD)=0.0647  p(t検定)=0.7617

置換検定の考え方が、この3枚でそのまま見えます。灰色のヒストグラムが「$H_0$ が正しいと仮定したときに統計量が取りうる値の分布」で、いずれも0を少しまたぐ位置に山があり、U統計量が負の値も取ることが確認できます。左の同一分布では赤い観測線が山のど真ん中にあり、中央の平均差では山の右端をはるかに越えています。右の二峰ケースは観測線が裾にかかっており、p値 0.065 という「惜しい」状態が視覚的に理解できます。 t検定のp値 0.76 がまったく反応していないのと対照的です。

3行から3つのことが読み取れます。同一分布ではMMDの推定値が $-0.0016$ と負になり（U統計量なので当然です）、p値 0.47 で棄却されません。平均差ケースでは両者とも有意ですが、t検定のほうがp値が小さい（0.0008 対 0.0050）。これは平均シフトという「t検定が想定した通りの違い」なので、専用の検定に軍配が上がるのは自然です。そして 二峰ケースでは、t検定のp値が 0.76 とまったく反応していないのに対し、MMDは 0.065 まで下がっています。この1回の試行では惜しくも 0.05 を切りませんでしたが、方向性は明確です。次にこれを多数回繰り返して検出力として定量化します。

実験4：検出力の比較

同じ手続きを300回繰り返し、$\alpha=0.05$ で棄却された割合（検出力）を測ります。4つのシナリオを用意しました。

import numpy as np
from scipy import stats

a = np.sqrt(3) / 2
def gen_null(r, k):  return r.normal(0, 1, (k, 1)), r.normal(0, 1, (k, 1))
def gen_mean(r, k):  return r.normal(0, 1, (k, 1)), r.normal(0.4, 1, (k, 1))
def gen_bimod(r, k):
    X = r.normal(0, 1, (k, 1))
    Y = (np.where(r.random(k) < .5, -a, a) + r.normal(0, .5, k))[:, None]
    return X, Y
def gen_laplace(r, k):
    return r.normal(0, 1, (k, 1)), r.laplace(0, 1 / np.sqrt(2), (k, 1))

def power(gen, trials=300, n=100, B=200, alpha=0.05, seed=7, sigma=None):
    r = np.random.default_rng(seed); hm = ht = 0
    for _ in range(trials):
        X, Y = gen(r, n)
        _, p, _, _ = mmd_permutation_test(X, Y, B, r, sigma)
        hm += (p <= alpha)
        ht += (stats.ttest_ind(X.ravel(), Y.ravel()).pvalue <= alpha)
    return hm / trials, ht / trials

for nm, g in [("H0: 同一分布", gen_null), ("平均差0.4", gen_mean),
              ("二峰(平均分散一致)", gen_bimod), ("ラプラス(平均分散一致)", gen_laplace)]:
    pm, pt = power(g)
    print(f"{nm:22s} MMD検定={pm:.3f}  t検定={pt:.3f}")

H0: 同一分布            MMD検定=0.040  t検定=0.030
平均差0.4              MMD検定=0.643  t検定=0.837
二峰(平均分散一致)      MMD検定=0.390  t検定=0.030
ラプラス(平均分散一致)  MMD検定=0.337  t検定=0.050

この棒グラフが記事全体の結論を1枚に凝縮しています。左端の $H_0$ では両者とも破線（$\alpha=0.05$）を超えず、検定として妥当であることが分かります。平均差では灰色（t検定）0.837 が緑（MMD）0.643 を上回りますが、右の2つでは灰色が破線の高さに張りついたまま で、二峰では 0.030、ラプラスでは 0.050 と有意水準そのものです。そこで緑だけが 0.39 / 0.34 と立ち上がっている——この非対称性こそが、MMDを使うべき場面と使うべきでない場面の境界線です。

この4行が、MMDを使う理由と、使わない理由の両方を語っています。

1行目（$H_0$）：どちらも第一種の過誤が 0.040 / 0.030 と名目水準 0.05 を超えていません。置換検定の $+1$ 補正が正しく効いており、MMD検定は妥当な検定になっています。

2行目（平均差）：t検定 0.837 に対しMMD 0.643 で、t検定が勝ちます。効果量 $d=0.4$、片群 $n=100$ のときt検定の理論検出力は約 0.805 で、実測 0.837 はこれと整合します。MMDは「あらゆる違いを探す」ぶん、特定の違いに特化した検定には検出力で劣る——これは避けられない代償です。平均シフトしか疑っていないなら、素直にt検定を使うべき なのです。

3行目・4行目（形状だけが違う）：ここで完全に立場が逆転します。二峰でもラプラスでも、平均と分散は $P$ と一致させてあるので t検定の検出力は 0.030 / 0.050 と、有意水準そのもの＝まったく無力 です。対してMMDは 0.390 / 0.337 と、明確に違いを見つけています。t検定が「気づけない」のではなく「見に行っていない」違いを、MMDは拾い上げているわけです。

まとめると、MMDは「何が違うか分からないが、とにかく違いを見つけたい」場面のための道具 です。分布ドリフト検知や生成モデルの評価は、まさにその状況にあたります。

実験5：帯域幅が検出力を左右する

理論節で述べた「$\sigma$ が広すぎると平均しか見えない」を確かめます。二峰ケース（平均・分散一致）で $\sigma$ を振って検出力を測ります。

import numpy as np

r = np.random.default_rng(3)
meds = [median_sigma(np.vstack(gen_bimod(r, 100))) for _ in range(50)]
print("メディアンヒューリスティックのsigma平均 =", round(float(np.mean(meds)), 3))

for s in [0.05, 0.1, 0.2, 0.4, 0.75, 1.5, 3.0, 10.0, None]:
    pm, _ = power(gen_bimod, trials=200, B=150, seed=3, sigma=s)
    label = "メディアン" if s is None else f"{s:6.2f}"
    print(f"  sigma={label}  検出力={pm:.3f}")

メディアンヒューリスティックのsigma平均 = 0.705
  sigma=  0.05  検出力=0.325
  sigma=  0.10  検出力=0.435
  sigma=  0.20  検出力=0.570
  sigma=  0.40  検出力=0.560
  sigma=  0.75  検出力=0.320
  sigma=  1.50  検出力=0.055
  sigma=  3.00  検出力=0.050
  sigma= 10.00  検出力=0.045
  sigma=メディアン  検出力=0.355

この1枚は、実務でいちばん効くパラメータが $\sigma$ であることを示しています。曲線のピークは点線で示した「二峰の山そのものの幅（標準偏差 0.5）」のすぐ左側にあり、見分けたい構造のスケールに帯域幅を合わせたときが最強 という原則がそのまま現れています。右側では急落して灰色の破線（$\alpha=0.05$）に張りつき、$\sigma=1.5$ 以上では検定がまったく働きません。赤い星のメディアンヒューリスティック（$\sigma\approx 0.71$）は 0.355 で、最適値 0.570 には届かないものの崖の手前で踏みとどまっており、「無調整の既定値」としては妥当な位置にいることが分かります。

検出力は $\sigma$ に対して 山なり になり、$\sigma=0.2 \sim 0.4$ で最大 0.57 に達します。二峰の各山の標準偏差が 0.5 であることを思い出すと、「見分けたい構造のスケール」に $\sigma$ を合わせたときが最強 だと納得できます。$\sigma$ を 1.5 以上に広げると検出力は 0.055、3.0 で 0.050、10.0 で 0.045 と、有意水準 0.05 と区別がつかなくなります。これは理論節で導いた「$\sigma\to\infty$ では平均差しか見なくなる」の直接的な証拠です（このケースでは平均差が0なので、検出力も $\alpha$ に落ちる）。

一方、$\sigma$ を 0.05 まで狭めても検出力は 0.325 と下がります。カーネルがデルタ関数に近づき、どの2点も「無関係」と判定されて情報が消えるためです。

そしてメディアンヒューリスティック（$\sigma\approx 0.705$）は 0.355。最適値の 0.57 には及ばないものの、無情報な選択（0.05や大きな値）よりはずっとまし という位置づけです。「まずメディアンで動かし、検出力が欲しければマルチカーネルやカーネル学習に進む」という実務の定石が、この一枚から見えてきます。

実験6：ドメイン適応の損失としてMMDを使う

最後に、MMDを検定ではなく 損失関数 として使います。設定はこうです。ソースドメイン（ラベルあり）は二日月形の2クラスデータ、ターゲットドメイン（ラベルなし）は同じデータを 35 度回転して平行移動したもの。ソースだけで学習した分類器をターゲットに当てると、当然うまくいきません。

そこで、中間特徴 $\phi(x)$ のソース側分布とターゲット側分布のMMDを損失に足します。

$$ \mathcal{L} = \underbrace{\frac{1}{m}\sum_i \ell\big(g(\phi(x_i^s)), y_i^s\big)}_{\text{ソースの分類損失}} + \lambda\, \underbrace{\widehat{\mathrm{MMD}}^2_b\big( \{\phi(x^s)\}, \{\phi(x^t)\} \big)}_{\text{特徴分布を揃える正則化}} $$

ここで V統計量（バイアスあり）を使う のがポイントでした。最小化したい量なので $\ge 0$ が保証されているほうが安全ですし、バイアス項は $\phi$ に依存しない定数ではないものの、$\ge 0$ かつ微分可能で勾配が素直に流れます。

import numpy as np, torch, torch.nn as nn
from sklearn.datasets import make_moons

def make_domains(seed=0, n=400, deg=35, shift=(0.6, 0.4)):
    Xs, ys = make_moons(n_samples=n, noise=0.12, random_state=seed)
    Xt, yt = make_moons(n_samples=n, noise=0.12, random_state=seed + 100)
    th = np.deg2rad(deg)
    R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
    Xt = Xt @ R.T + np.array(shift)          # 回転＋平行移動でドメインシフトを作る
    return Xs.astype(np.float32), ys.astype(np.int64), Xt.astype(np.float32), yt.astype(np.int64)

class Net(nn.Module):
    def __init__(self):
        super().__init__()
        self.feat = nn.Sequential(nn.Linear(2, 32), nn.ReLU(), nn.Linear(32, 2))
        self.clf = nn.Linear(2, 2)
    def forward(self, x):
        f = self.feat(x)
        return f, self.clf(f)

def mmd2_biased_torch(X, Y, sigmas):
    """マルチカーネルMMD（V統計量）。torchで微分可能"""
    Z = torch.cat([X, Y], 0)
    d2 = torch.cdist(Z, Z) ** 2
    m = len(X)
    K = sum(torch.exp(-d2 / (2 * s ** 2)) for s in sigmas) / len(sigmas)
    return K[:m, :m].mean() + K[m:, m:].mean() - 2 * K[:m, m:].mean()

def median_sigma_t(Z):
    d = torch.cdist(Z, Z)
    return (d[d > 0].median() / np.sqrt(2)).clamp(min=1e-3)

学習ループでは、毎ステップで特徴空間のメディアン帯域幅を測り直し、その 0.5 / 1 / 2 倍の3つのカーネルを平均して使います。特徴が学習中に伸縮するので、固定 $\sigma$ だと途中でスケールが合わなくなるためです。

import numpy as np, torch, torch.nn as nn

def train(lam, seed=0, epochs=600):
    torch.manual_seed(seed); np.random.seed(seed)
    Xs, ys, Xt, yt = make_domains(seed)
    Xs_t, ys_t, Xt_t = torch.tensor(Xs), torch.tensor(ys), torch.tensor(Xt)
    net = Net(); opt = torch.optim.Adam(net.parameters(), lr=2e-3)
    ce = nn.CrossEntropyLoss()
    for _ in range(epochs):
        fs, ls = net(Xs_t); ft, _ = net(Xt_t)
        loss = ce(ls, ys_t)                       # ソースの分類損失（ラベルを使う）
        if lam > 0:
            base = median_sigma_t(torch.cat([fs, ft]).detach()).item()
            loss = loss + lam * mmd2_biased_torch(fs, ft, [base * c for c in (0.5, 1, 2)])
        opt.zero_grad(); loss.backward(); opt.step()
    net.eval()
    with torch.no_grad():
        fs, ls = net(Xs_t); ft, lt = net(Xt_t)
        acc_s = (ls.argmax(1).numpy() == ys).mean()
        acc_t = (lt.argmax(1).numpy() == yt).mean()   # ターゲットのラベルは評価にのみ使用
        sg = median_sigma_t(torch.cat([fs, ft])).item()
        mmd_feat = mmd2_biased_torch(fs, ft, [sg]).item()
    return acc_s, acc_t, mmd_feat, fs.numpy(), ft.numpy()

# 比較のため、生の入力空間でのMMDも測っておく
Xs, ys, Xt, yt = make_domains(0)
sg0 = median_sigma_t(torch.tensor(np.vstack([Xs, Xt]))).item()
print("入力空間のMMD^2 =",
      round(mmd2_biased_torch(torch.tensor(Xs), torch.tensor(Xt), [sg0]).item(), 4))

for lam in [0.0, 0.5, 1.0, 2.0, 5.0]:
    r = [train(lam, seed=s) for s in range(5)]      # 5シードの平均
    print(f"lambda={lam:4.1f}  ソース精度={np.mean([a[0] for a in r]):.3f}"
          f"  ターゲット精度={np.mean([a[1] for a in r]):.3f}"
          f"  特徴のMMD^2={np.mean([a[2] for a in r]):.4f}")

入力空間のMMD^2 = 0.1771
lambda= 0.0  ソース精度=0.996  ターゲット精度=0.527  特徴のMMD^2=0.1675
lambda= 0.5  ソース精度=0.994  ターゲット精度=0.714  特徴のMMD^2=0.0134
lambda= 1.0  ソース精度=0.992  ターゲット精度=0.715  特徴のMMD^2=0.0099
lambda= 2.0  ソース精度=0.980  ターゲット精度=0.725  特徴のMMD^2=0.0067
lambda= 5.0  ソース精度=0.934  ターゲット精度=0.741  特徴のMMD^2=0.0030

①では回転と平行移動によって青（ソース）と赤（ターゲット）がずれている様子が見えます。②の左（$\lambda=0$）では中間特徴も同じようにずれたままで、差し込みの主成分方向のヒストグラムを見ると青だけが右側にもう1つ山を持っています。これが「ソースで完璧な決定境界がターゲットでは通用しない」状態の正体です。②の右（$\lambda=2$）では2色のヒストグラムがほぼ重なり、MMD$^2$ も 0.18 から 0.008 へ落ちました。③のグラフでは、$\lambda$ を上げるほど赤い破線（特徴のMMD$^2$）が下がって緑（ターゲット精度）が上がる一方、青（ソース精度）が $\lambda=5$ で 0.93 まで削られていく トレードオフが読み取れます。

この表には、MMD正則化の効果と限界が両方はっきり出ています。$\lambda=0$（ソースのみで学習）ではターゲット精度が 0.527 と、2クラス問題のチャンスレベル 0.5 とほぼ変わりません。 ソース精度は 0.996 なので「学習が失敗した」のではなく、「ソースで完璧に働く境界がターゲットではまったく通用しない」という典型的なドメインシフトです。そして特徴空間のMMD$^2$ は 0.1675 と、入力空間の 0.177 とほぼ同じ。普通に学習すると、ネットワークはドメインの違いを特徴空間に持ち込んだままにする のです。

$\lambda$ を入れると状況が一変します。$\lambda=0.5$ で特徴のMMD$^2$ は 0.1675 → 0.0134 と 12分の1以下 に落ち、ターゲット精度は 0.527 → 0.714 へ大きく改善しました。ターゲットのラベルは学習に一切使っていない ことに注意してください。「特徴分布を揃えろ」という制約だけで、これだけの転移が起きています。

一方、$\lambda$ を上げ続けても効果は頭打ちです。$\lambda=5$ では特徴のMMD$^2$ が 0.0030 まで下がり、ターゲット精度も 0.741 とわずかに上がりますが、ソース精度が 0.996 → 0.934 と犠牲になっています。分布を揃えることを優先しすぎて、クラスを分ける情報まで潰し始めたのです。$\lambda$ は「揃える強さ」と「分類できる情報を残す強さ」のトレードオフを調整するツマミであり、実務では検証セットで選ぶ必要があります。

さらに重要な限界として、改善後のターゲット精度も 0.71 止まり で、ソースの 0.99 には遠く及びません。理由は明快です。MMDが揃えるのは特徴の 周辺分布 $p(\phi(x))$ だけであり、クラス条件付き分布 $p(\phi(x)\mid y)$ が揃う保証はないからです。回転が大きいと「ソースのクラス0の雲」と「ターゲットのクラス1の雲」を重ねても周辺分布としては揃ってしまい、ラベルが入れ替わったまま整合してしまいます。「分布を揃えれば転移する」は必要条件であって十分条件ではない——この事実は、条件付き分布まで揃えるCDAN、疑似ラベルを併用する自己学習、DANN（敵対的ドメイン適応） といった後続研究の出発点になりました。

MMDの限界と発展

実験を通じて見えてきた注意点を、理論的な背景と合わせて整理しておきます。

(1) 高次元で検出力が落ちる。 MMDの値そのものは次元によらず定義できますが、同じ $\mathrm{MMD}$ の大きさを検出するのに必要なサンプル数は次元とともに増えます。直感的には、高次元では「典型的なペアの距離」がどんどん大きくなり、しかもその分散が相対的に小さくなる（距離の集中）ため、メディアンヒューリスティックで選んだ $\sigma$ のもとではどのペアも似たようなカーネル値になってしまうからです。生の1000次元ベクトルに直接MMDを当てるより、まず学習した特徴空間に落としてからMMDを測る のが実務の定石です。

(2) カーネル選択が検出力を支配する。 実験5で見たように、$\sigma$ の選び方だけで検出力が 0.05 から 0.57 まで動きました。マルチカーネルMMD（複数の $\sigma$ の凸結合）や、深層カーネル（ニューラルネットで特徴を作ってからガウスカーネル）を学習する手法が提案されています。ただし 同じデータでカーネルを選んで検定すると第一種の過誤が壊れる ので、データを分割するか、選択を含めた補正が必要です。

(3) 分布が遠すぎると勾配が消える。 ガウスカーネルは $\|x-y\|$ が $\sigma$ より十分大きいと $k\approx 0$ に飽和します。すると $\mathrm{MMD}^2 \approx \mathbb{E}[k(X,X’)] + \mathbb{E}[k(Y,Y’)]$ という「相手の位置に依存しない」値になり、$Q$ を $P$ に近づける方向の勾配がほぼ0になります。生成モデルの学習でMMDをそのまま損失にすると学習が進まないことがあるのはこのためで、Wasserstein距離 が好まれる場面がある理由の1つでもあります。適応的な帯域幅やマルチカーネルはこの問題の緩和策です。

(4) 計算量 $O(n^2)$。 カーネル行列がボトルネックです。線形時間推定量、ブロック分割（$B$ 個のブロックでU統計量を計算して平均）、Random Fourier Featuresによる近似（$\varphi$ を有限次元でサンプリングして平均埋め込みをベクトルとして持つ）などが使われます。

(5) 発展：HSICと独立性検定。 同時分布 $P_{XY}$ と周辺分布の積 $P_X \otimes P_Y$ のあいだのMMDを取ると、HSIC（Hilbert–Schmidt Independence Criterion） になります。$\mathrm{HSIC}=0 \iff X \perp Y$ が特性カーネルのもとで成り立つので、これは「非線形な相関」まで検出できる独立性の指標です。因果探索・特徴選択・公平性制約などに広く使われています。「二標本検定の道具が、そのまま独立性検定になる」——平均埋め込みという枠組みの汎用性を示す好例です。

まとめ

本記事では、MMD（最大平均差異）について、定義から導出、推定、検定、応用まで解説しました。

• 定義：MMDは「RKHSの単位球を証人クラスに取った積分確率距離」であり、$\sup$ はCauchy–Schwarzで解けて $\mathrm{MMD}(P,Q)=\|\mu_P-\mu_Q\|_{\mathcal{H}}$ という平均埋め込みの距離になる
• 3項表現：再生性を2回使うと $\mathrm{MMD}^2 = \mathbb{E}[k(X,X’)] – 2\mathbb{E}[k(X,Y)] + \mathbb{E}[k(Y,Y’)]$ と、カーネル期待値だけで書ける。密度推定は一切不要
• 特性カーネル：$\mathrm{MMD}=0 \iff P=Q$ が成り立つのは $P\mapsto\mu_P$ が単射なとき。並進不変カーネルではBochner表現 $\mathrm{MMD}^2=\int|\varphi_P-\varphi_Q|^2 d\Lambda$ から、スペクトル測度の台が全域なら特性 と分かる。ガウス・ラプラスはOK、線形・多項式は有限次モーメントしか見ない
• U統計量とV統計量：対角を除けば不偏（負にもなり得る、検定向き）、含めればバイアス $+\frac{2}{m}(1-\mathbb{E}[k(X,X’)])$ が乗るが常に $\ge 0$（損失関数向き）。$m=100$ でバイアス 0.0085 という理論値を実測で確認した
• 帯域幅：メディアンヒューリスティックは無難だが最適ではない。二峰の識別で最適 $\sigma$ の検出力 0.57 に対しメディアンは 0.355、$\sigma$ が大きすぎると 0.05（＝無力）まで落ちる
• 検定：帰無分布が分布依存の重み付きカイ二乗和になるため、実務では置換検定を使う。カーネル行列を1回作って添字を並べ替えるだけで済む
• 検出力：平均シフトではt検定に負ける（0.643 対 0.837）が、平均・分散が一致した二峰やラプラスではt検定が無力（0.03〜0.05）なのに対しMMDは 0.34〜0.39 を出す
• ドメイン適応：特徴分布のMMDを損失に足すと、ターゲットのラベルなしで精度 0.527 → 0.715 へ改善。ただし周辺分布しか揃えないので万能ではない

MMDは「カーネルを1つ選ぶだけで、分布の全情報を1つのスカラーに落とす」という、シンプルさと理論的な強さを兼ね備えた道具です。二標本検定、ドメイン適応、独立性検定（HSIC）、生成モデルの評価と、応用範囲は今も広がり続けています。

次のステップとして、以下の記事も参考にしてください。

• 再生核ヒルベルト空間（RKHS）とは？わかりやすく解説 — 平均埋め込みの土台となる空間の性質
• カーネル関数の種類と性質 — ガウス過程で使うカーネルを徹底解説 — 特性カーネルかどうかを見分けるための基礎
• 最適輸送とWasserstein距離をわかりやすく：Kantorovich双対とSinkhornまで — MMDと並ぶもう1つの積分確率距離
• DANN（敵対的ドメイン適応） — 勾配反転でドメイン不変な特徴を学ぶ — MMD正則化の敵対的な代替手法
• 並べ替え検定（置換検定）の理論と実装をわかりやすく解説 — p値の計算に使った枠組みの一般論
• KLダイバージェンスとJSダイバージェンス — 確率分布間の距離と機械学習への応用 — 密度比ベースの距離との比較

参考文献

• A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, A. Smola, “A Kernel Two-Sample Test,” Journal of Machine Learning Research, 13(25):723–773, 2012.
• B. K. Sriperumbudur, A. Gretton, K. Fukumizu, B. Schölkopf, G. R. G. Lanckriet, “Hilbert Space Embeddings and Metrics on Probability Measures,” JMLR, 11:1517–1561, 2010.
• D. Sejdinovic, B. Sriperumbudur, A. Gretton, K. Fukumizu, “Equivalence of distance-based and RKHS-based statistics in hypothesis testing,” Annals of Statistics, 41(5):2263–2291, 2013.
• M. Long, Y. Cao, J. Wang, M. I. Jordan, “Learning Transferable Features with Deep Adaptation Networks,” ICML, 2015.
• A. Ramdas, S. J. Reddi, B. Póczos, A. Singh, L. Wasserman, “On the Decreasing Power of Kernel and Distance Based Nonparametric Hypothesis Tests in High Dimensions,” AAAI, 2015.