「サンプルサイズを十分大きく取れば、標本平均は真の期待値に収束する」——大数の法則が保証してくれるこの安心感は、データ分析の出発点としてあまりに自然です。ところが世の中には、$n$ を $10^6$ にしても $10^9$ にしても標本平均が暴れ続け、決して収束しない確率分布が存在します。その代表格が コーシー分布（Cauchy distribution） です。
コーシー分布は、見た目こそ正規分布によく似た左右対称の釣鐘型をしていますが、裾の落ち方が桁違いに緩やかで、平均も分散も「数学的に存在しない」という性質を持ちます。中心極限定理の前提（有限分散）を破る分布の代表例であり、確率論を学ぶ上で必ず一度は通るべき「異常な」分布です。
この奇妙な分布は、しかし実応用の最前線でひっきりなしに登場します。素粒子物理の 共鳴ピーク（Breit-Wigner分布）として観測される散乱断面積の形、原子分光や NMRスペクトル の自然線幅（ローレンツ分布として）、そして外れ値の多い計測データに対する ロバスト回帰 の損失関数（コーシー損失）。さらに、中央値フィルタ が画像処理で頑健に動く理論的根拠も、コーシー型ノイズに対する最適性に結びついています。「平均が壊れる世界」を理解することは、現実の汚いデータを扱うすべての場面で実用的な意味を持つのです。
本記事の内容
- コーシー分布の確率密度関数の導出と、ローレンツ分布／Breit-Wigner分布／自由度1のt分布としての複数の顔
- 平均・分散が「存在しない」とはどういうことか — 積分が発散する具体的な計算
- 中央値とMADは有限であるという救い、そしてロバスト統計量がなぜ必要か
- 標本平均がn→∞でも収束しないというショッキングな事実のMonte Carlo検証
- 安定分布族としての位置づけと、2つの独立な標準正規の比としての構成
- M推定、Huber損失とコーシー損失を使った外れ値混入データへのロバスト線形回帰のPython実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
直感: 正規分布が破綻する場面
データ分析の教科書の最初のページに登場する正規分布は、あらゆる現象を支配しているように見えます。身長、体重、測定誤差、株価のリターン……「いざとなれば中心極限定理が助けてくれる」「サンプルを増やせば標本平均は真値に収束する」というのは、ほとんど自然法則のように感じられるでしょう。
ところが、そう感じてしまうのは、私たちが「裾が指数関数的に減衰する分布」しか普段見ていないからです。正規分布の確率密度は $\exp(-x^2/2)$ で減衰しますから、$x = 5$ の確率はもう $10^{-7}$ 程度しかありません。$x = 10$ なら $10^{-23}$ で、宇宙の年齢を全宇宙のサンプリングに使っても見つかりません。だから「外れ値はめったに来ない」「来ても少しずつしか平均を動かさない」と直感が育つわけです。
ところが、自然界には裾が多項式でしか減衰しない分布が存在します。コーシー分布の確率密度は大きな $x$ で $1/x^2$ のオーダーで減衰するだけ。これは何を意味するかというと、$|x| > 100$ となるサンプルが、宝くじではなく日常的にやってくる、ということです。さらに $|x| > 10^6$ も「滅多にない」程度には起こります。こうした極端値が混じるとき、標本平均は1つの巨大なサンプルによって簡単に乗っ取られ、いくらサンプル数を増やしても落ち着かないのです。
この事情を象徴するのが、19世紀のフランスの数学者ポアソンとコーシーの間で交わされた逸話です。ポアソンは「正規分布近似は普遍的だ」と主張しました。コーシーはそれに対し、「いや、こういう分布もある」と $f(x) \propto 1/(1+x^2)$ という反例を提示しました。これが今でいうコーシー分布で、当時の確率論の常識を覆したのです。本記事では、その「常識破り」の中身を一つひとつ数式と数値実験で確かめていきます。
左図は「光源から地面へランダムな角度で投影したとき、着地点がコーシー分布に従う」という幾何的な直感を示しています。角度が $\pm 90°$ に近づくほど着地点が遠くへ飛び、これが重い裾を生む仕組みです。右図では同じ事実を確率密度として表現しており、コーシー分布（赤）が中心ピーク付近で正規分布（青破線）より低く、しかし裾では遥かに高い値を保つことが一目で分かります。
まずはこの分布の数学的な姿——確率密度関数とその導出——から見ていきましょう。
コーシー分布の確率密度関数と導出
定義
コーシー分布は、位置パラメータ $x_0$（中央の位置）と尺度パラメータ $\gamma > 0$（半値半幅 HWHM）で特徴づけられ、確率密度関数は次で与えられます。
$$ f(x; x_0, \gamma) = \frac{1}{\pi \gamma}\, \frac{1}{1 + \left(\dfrac{x – x_0}{\gamma}\right)^2} $$
$x_0 = 0,\ \gamma = 1$ の場合を 標準コーシー分布 と呼び、
$$ f(x) = \frac{1}{\pi (1 + x^2)} $$
という非常にシンプルな形になります。グラフは正規分布と同じく左右対称・単峰の釣鐘型ですが、裾の落ち方が決定的に異なります。正規分布が $\exp(-x^2/2)$ で 指数的に 減衰するのに対し、コーシー分布は $1/x^2$ で 多項式的に しか減衰しません。これが「重い裾（heavy tail）」と呼ばれる所以です。
正規化定数の確認
定義式の係数 $1/(\pi\gamma)$ がきちんと確率密度になっていること、すなわち $\int_{-\infty}^{\infty} f(x)\, dx = 1$ を確かめておきます。標準コーシー（$x_0 = 0, \gamma = 1$）で示せば十分です。
$$ \int_{-\infty}^{\infty} \frac{1}{\pi(1 + x^2)}\, dx = \frac{1}{\pi}\Big[\arctan x\Big]_{-\infty}^{\infty} = \frac{1}{\pi}\left(\frac{\pi}{2} – \left(-\frac{\pi}{2}\right)\right) = 1 $$
$\arctan$ の積分公式 $\int dx/(1+x^2) = \arctan x$ を使い、両端で $\pm\pi/2$ になることから値が綺麗に $1$ になります。一般の $\gamma$ の場合も $u = (x – x_0)/\gamma$ と置換すれば同じ計算に帰着し、係数 $1/(\pi\gamma)$ が必然的に出てくることが確認できます。
幾何学的な由来（投影定理による導出）
コーシー分布は無味乾燥な数式として与えられるのではなく、非常に自然な幾何学的状況から「自動的に」生まれてきます。
地面（数直線 $y = 0$）から高さ $\gamma$ の位置 $(x_0,\gamma)$ に光源があり、その光源が一様にランダムな方向で光を放つとします。具体的には、光源から地面に向かう光線が鉛直下方となす角度 $\theta$ が、区間 $[-\pi/2, \pi/2]$ で一様分布に従う、と仮定します。このとき、光線が地面に当たる点 $X$ の分布が、まさにコーシー分布になります。
光線が地面に達する点の座標は、簡単な三角関数で
$$ X = x_0 + \gamma \tan\theta $$
と書けます。$\theta$ が一様分布 $U(-\pi/2, \pi/2)$ に従うので、その確率密度は
$$ f_\theta(\theta) = \frac{1}{\pi}, \qquad \theta \in \left(-\tfrac{\pi}{2}, \tfrac{\pi}{2}\right) $$
です。変数変換の公式 を使うと、$x = x_0 + \gamma\tan\theta$ の逆関数は $\theta = \arctan\!\big((x – x_0)/\gamma\big)$ であり、その微分は
$$ \frac{d\theta}{dx} = \frac{1}{\gamma}\, \frac{1}{1 + \big((x – x_0)/\gamma\big)^2} $$
となります。確率密度の変換則 $f_X(x) = f_\theta(\theta(x))\, |d\theta/dx|$ を適用すると、
$$ f_X(x) = \frac{1}{\pi} \cdot \frac{1}{\gamma}\, \frac{1}{1 + \big((x – x_0)/\gamma\big)^2} = \frac{1}{\pi\gamma}\, \frac{1}{1 + \big((x – x_0)/\gamma\big)^2} $$
と、まさにコーシー分布の式が現れました。$\tan\theta$ という非線形変換が、一様分布から重い裾を生み出しているのです。$\theta \to \pm\pi/2$ で $\tan\theta \to \pm\infty$ となり、角度がほんの少し $\pi/2$ に近づくだけで $X$ が爆発的に大きくなる——これが「外れ値が日常的にやってくる」直感的説明です。
2つの独立な標準正規の比としての構成
コーシー分布のもう一つの自然な構成法があります。$U, V$ が独立にそれぞれ標準正規分布 $\mathcal N(0,1)$ に従うとき、その比 $Z = U/V$ は標準コーシー分布に従う、という事実です。証明は変数変換 $(U, V) \to (Z, V)$ のヤコビアンを計算し、$V$ について周辺化することで得られます。
$$ f_Z(z) = \int_{-\infty}^{\infty} |v|\, f_U(zv)\, f_V(v)\, dv = \int_{-\infty}^{\infty} \frac{|v|}{2\pi}\, e^{-(z^2+1)v^2/2}\, dv $$
ここで $|v|$ の対称性から正の半区間の2倍に書き換え、$s = (z^2 + 1)v^2/2$ と置換すると、$ds = (z^2 + 1)v\, dv$ なので
$$ f_Z(z) = \frac{1}{\pi}\int_0^{\infty} e^{-(z^2+1)v^2/2}\, v\, dv = \frac{1}{\pi(z^2+1)}\int_0^\infty e^{-s}\, ds = \frac{1}{\pi(1 + z^2)} $$
と標準コーシー分布の密度が現れます。「2つの独立な正規ノイズの比」がコーシー分布を生むという事実は、たとえば角度推定（$\arctan(U/V)$ のような形）が一様分布になることや、正規分布の比が出てくるあらゆる場面でコーシー分布が現れることを意味します。
別名と物理での顔
同じ分布が、応用分野によって異なる名前で呼ばれます。
- ローレンツ分布（Lorentzian）: 分光学・NMR・原子物理で、共鳴線の自然線幅を記述する関数。励起状態の寿命 $\tau$ と不確定性関係から、線幅 $\gamma$ が $\hbar/\tau$ と結びつきます。
- Breit-Wigner分布: 素粒子・原子核物理で、共鳴粒子（$Z$ボソン、$\rho$中間子、励起原子核など）の散乱断面積のエネルギー依存性を記述する関数。共鳴エネルギー $E_0$ と幅 $\Gamma$ を使い、断面積 $\sigma(E) \propto 1/[(E-E_0)^2 + \Gamma^2/4]$ という形で現れます。
- 自由度1のt分布: スチューデントのt分布で自由度 $\nu = 1$ としたものが、ぴったり標準コーシー分布になります。自由度 $\nu \to \infty$ で標準正規に近づき、$\nu = 1$ では最も重い裾を持つ——コーシーはt分布族の「もう一方の端」なのです。
これら3つの「同じ分布の別名」を意識しておくと、物理の論文と統計の教科書を行き来するときに視点がつながります。
左の線形スケールでは3分布の外見の違いが分かりにくいですが、右の対数スケールに切り替えると裾の落ち方の差が20桁以上に及ぶことが鮮明に現れます。正規分布（青）が $|x|=10$ 付近で $10^{-23}$ 以下に消えるのに対し、コーシー分布（赤）は $10^{-3}$ 近くを維持しており、「重い裾」が単なる定性的な言葉でなく桁外れの量的差異であることが確認できます。t分布（緑, df=3）はちょうどその中間に位置します。
次に、この分布の最大の特徴——平均と分散が「存在しない」ことの中身を、積分計算で具体的に見ていきましょう。
なぜ平均・分散が存在しないか
「期待値が存在しない」の正確な意味
確率論で「平均（期待値）が存在しない」とは、単に「計算しにくい」とか「無限大になる」という曖昧な意味ではなく、明確にルベーグ積分の意味で定義できないということを指します。期待値 $E[X] = \int x f(x)\, dx$ が意味を持つためには、まず $\int |x| f(x)\, dx < \infty$（絶対可積分性）が要求されます。これが成り立たないとき、$\int x f(x)\, dx$ の値はそもそも一意に定まりません。
なぜ絶対可積分性が必要か。積分を $x > 0$ の部分と $x < 0$ の部分に分けると、
$$ \int x f(x)\, dx = \int_0^\infty x f(x)\, dx – \int_0^\infty x f(-x)\, dx $$
となります。コーシー分布のように両者がそれぞれ $+\infty$ になる場合、$+\infty – \infty$ という未定義の演算が生じ、和の順序や積分範囲の取り方によって値が変わってしまいます。これは数学的にきちんと定義された量とは呼べません。
標準コーシー分布の期待値が存在しないこと
標準コーシー分布で、絶対可積分性が破れていることを実際に計算します。
$$ \int_{-\infty}^{\infty} |x|\, f(x)\, dx = \int_{-\infty}^{\infty} \frac{|x|}{\pi(1 + x^2)}\, dx = \frac{2}{\pi}\int_0^{\infty} \frac{x}{1 + x^2}\, dx $$
対称性から $|x|$ の積分は正の半区間の2倍になっています。ここで $u = 1 + x^2$、$du = 2x\, dx$ と置換すると、
$$ \frac{2}{\pi}\int_0^{\infty} \frac{x}{1 + x^2}\, dx = \frac{1}{\pi}\int_1^{\infty} \frac{du}{u} = \frac{1}{\pi}\Big[\log u\Big]_1^{\infty} = \infty $$
$\log u$ は $u \to \infty$ で発散します。すなわち $\int |x| f(x)\, dx = \infty$ となり、絶対可積分性が破綻するので、期待値 $E[X]$ は存在しません。
左のグラフでは、被積分関数 $|x|f(x)$ がコーシー分布（赤）では正規分布（青）よりはるかにゆっくり減衰することが分かります。正規分布の $|x|f(x)$ は遠方で急速に0へ向かうのに対し、コーシー分布は $|x|/(\pi(1+x^2)) \approx 1/(\pi x)$ とごく緩やかにしか下がりません。右の累積積分では、その違いが決定的な形で現れています。正規分布の積分は $\sqrt{2/\pi} \approx 0.798$（点線）という有限値へ収束しますが、コーシー分布の積分は右端まで発散し続け、打ち止めが来ません。これが「コーシー分布に平均が存在しない」という数学的事実の視覚的証拠です。
コーシー主値による「対称性の誘惑」
「対称な分布なんだから平均はゼロでしょう？」と直感的に答えたくなります。実際、コーシー主値（Cauchy principal value）を取れば
$$ \mathrm{p.v.}\!\int_{-\infty}^{\infty} \frac{x}{\pi(1+x^2)}\, dx = \lim_{R\to\infty}\int_{-R}^{R} \frac{x}{\pi(1+x^2)}\, dx = \lim_{R\to\infty} 0 = 0 $$
となります。被積分関数が奇関数で、対称な区間で打ち消し合うからです。しかしこの「主値」は積分範囲を対称に取った場合に限った値であり、たとえば $[-R, 2R]$ のように非対称に広げると別の値（実際には $\frac{1}{\pi}\log 2$）に収束します。期待値が一意に定まるためには、積分範囲の取り方に依らずに同じ値になる必要があり、コーシー分布はその要件を満たさないため、「平均はゼロ」と言ってはいけないのです。
分散も存在しない
分散 $\mathrm{Var}(X) = E[X^2] – (E[X])^2$ は、まず $E[X^2]$ が定義される必要があります。標準コーシー分布で計算すると、
$$ \int_{-\infty}^{\infty} x^2 f(x)\, dx = \frac{1}{\pi}\int_{-\infty}^{\infty} \frac{x^2}{1 + x^2}\, dx $$
被積分関数を変形すると $\frac{x^2}{1+x^2} = 1 – \frac{1}{1+x^2}$ となり、
$$ \frac{1}{\pi}\int_{-\infty}^{\infty}\left(1 – \frac{1}{1 + x^2}\right) dx = \frac{1}{\pi}\Big[(x – \arctan x)\Big]_{-\infty}^{\infty} $$
$x$ の項が両端で $\pm\infty$ になり、$\arctan x$ は有界なので、この積分は発散します。よって $E[X^2] = \infty$ で、分散も存在しません。同様の議論で任意の正のべき乗 $E[|X|^p]$（$p \geq 1$）が存在しないことが示せます——平均だけでなく、すべての「整数モーメント」が壊れているのです。
しかし中央値とMADは存在する
平均と分散が壊れたなら、コーシー分布は手も足も出ない無秩序な分布なのでしょうか。実はそうではありません。中央値（median）と中央絶対偏差（MAD: Median Absolute Deviation）は、コーシー分布でも美しく定義可能です。
中央値は累積分布関数が $1/2$ となる点として定義されます。コーシー分布の累積分布関数は、密度を積分して
$$ F(x) = \int_{-\infty}^x \frac{1}{\pi\gamma(1 + ((t-x_0)/\gamma)^2)}\, dt = \frac{1}{\pi}\arctan\!\left(\frac{x – x_0}{\gamma}\right) + \frac{1}{2} $$
となります。$F(x) = 1/2$ を解くと $\arctan((x-x_0)/\gamma) = 0$ より $x = x_0$。すなわち中央値はちょうど位置パラメータ $x_0$ に一致します。
同様に、$|X – x_0|$ の中央値（MAD）も次のように計算できます。$|X – x_0| \leq m$ となる確率が $1/2$ になる $m$ を求めればよいので、
$$ P(|X – x_0| \leq m) = F(x_0 + m) – F(x_0 – m) = \frac{2}{\pi}\arctan(m/\gamma) = \frac{1}{2} $$
これを解くと $\arctan(m/\gamma) = \pi/4$ より $m = \gamma$。MADは尺度パラメータ $\gamma$ そのものになります。位置と尺度を表す「正しい」統計量が、コーシー分布では平均・標準偏差ではなく、中央値・MADなのです。
ここで重要な教訓が浮かび上がります。重い裾を持つ分布では、平均と分散を計算してはいけない（できない）。代わりに中央値・MADを使え、と。これがロバスト統計の出発点であり、後ほど詳しく扱います。
しかし、平均が存在しないという数学的事実が、実際のデータ分析にどんな衝撃をもたらすのか——その威力を最も強烈に示すのが「標本平均が収束しない」という現象です。次にそれを実演します。
2つの独立な標準正規の比としての構成
コーシー分布の理論的位置づけをさらに確かめるために、「$U, V \sim \mathcal{N}(0,1)$ 独立のとき $Z = U/V$ はコーシー分布に従う」という事実をシミュレーションで検証してみましょう。
50万サンプルの $U/V$ のヒストグラム（青）と、理論的なコーシー分布のPDF（赤線）が極めてよく重なっています。$|z| < 15$ の範囲に絞っても、全体の分布形状がコーシー分布に沿っていることが確認でき、「2正規の比」という非常にシンプルな構成からコーシー分布が生まれる事実が数値実証できました。この性質は、測定ノイズの比やアンテナの受信信号の比など、比の演算が出てくる場面でコーシー分布が頻繁に現れる理由を直接説明しています。
標本平均が収束しない驚き
大数の法則と中心極限定理が要求するもの
通常の確率分布で、私たちが何気なく頼っている2つの定理を確認しておきます。
- 大数の法則（弱／強）: $X_1, X_2, \dots$ が独立同分布で $E[|X|] < \infty$ ならば、標本平均 $\bar X_n = (X_1 + \cdots + X_n)/n$ は $n \to \infty$ で真の期待値 $\mu = E[X]$ に収束する。
- 中心極限定理: さらに $\mathrm{Var}(X) = \sigma^2 < \infty$ ならば、$\sqrt n (\bar X_n - \mu)$ は $\mathcal N(0, \sigma^2)$ に分布収束する。
どちらの定理も、その前提として「期待値が存在する」「分散が有限である」ことを要求しています。コーシー分布はその前提を満たしません。したがって大数の法則も中心極限定理も適用できない——前のセクションの議論は、こうして実用的な帰結に結びつくのです。
標本平均がまたコーシー分布になる
コーシー分布の最も衝撃的な性質は、次の事実です。
$X_1, X_2, \dots, X_n$ が独立に標準コーシー分布に従うとき、その標本平均 $\bar X_n = (X_1 + \cdots + X_n)/n$ もまた標準コーシー分布に従う。
すなわち $\bar X_n \overset{d}{=} X_1$。標本数を $n$ に増やしても、平均の分布は1サンプルの分布と全く同じ。$n$ を $10^6$ にしても $10^9$ にしても、$\bar X_n$ の散らばりは縮まりません。
この事実を視覚的に示す前に、まず大数の法則の成否を正規分布との比較で確認しておきましょう。
左の正規分布では5つのトライアルすべてが $n$ の増加とともにきれいに真の平均 $0$ へ収束しており、$n = 10^5$ ではほぼ重なっています——これが大数の法則の本来の姿です。右のコーシー分布は全く異なる景色を見せます。累積平均はサンプル中に紛れ込む巨大な値のたびに大きくジャンプし、$n$ を増やしても一向に安定しません。試行ごとに全く異なる軌跡を描いており、「大きくなれば収束する」という直感が完全に崩壊しています。
この事実は 特性関数 から最もエレガントに示せます。標準コーシー分布の特性関数は
$$ \varphi(t) = E[e^{itX}] = \int_{-\infty}^{\infty}\frac{e^{itx}}{\pi(1+x^2)}\, dx = e^{-|t|} $$
です（複素積分で証明される標準的な結果）。独立和の特性関数は積になるので、$S_n = X_1 + \cdots + X_n$ の特性関数は
$$ \varphi_{S_n}(t) = \big(e^{-|t|}\big)^n = e^{-n|t|} $$
となります。これは尺度パラメータ $n$ のコーシー分布の特性関数です。標本平均 $\bar X_n = S_n / n$ の特性関数は $\varphi_{\bar X_n}(t) = \varphi_{S_n}(t/n) = e^{-n|t/n|} = e^{-|t|}$ で、これは元の標準コーシー分布の特性関数とぴったり同じ。よって $\bar X_n$ は標準コーシー分布に従います。
直感的には、$n$ 個のサンプルの中にたまに紛れ込む巨大な値が、$n$ で割られても「平均を一気に動かすに十分な大きさ」を保ってしまう、ということです。指数的に減衰する裾を持つ正規分布なら巨大値はほぼ来ませんが、$1/x^2$ の裾は $n$ サンプル中の最大値が $n$ のオーダーで成長する（$\max|X_i| \sim n$）ため、平均を計算しても1サンプル分の散らばりが残るのです。
標本中央値は収束する
一方、標本中央値はちゃんと収束します。これは、中央値が累積分布関数 $F$ の $1/2$ 分位点として定まり、極端値の影響を受けないからです。コーシー分布の標本中央値 $\tilde X_n$ は、大標本では漸近的に正規分布
$$ \tilde X_n \overset{d}{\approx} \mathcal N\!\left(x_0,\, \frac{1}{4 n f(x_0)^2}\right) = \mathcal N\!\left(x_0,\, \frac{\pi^2 \gamma^2}{4n}\right) $$
に従います（中央値の漸近正規性、Bahadur表現）。分散は $1/n$ で減るので、サンプル数を増やせばちゃんと精度が上がります。平均は壊れるが中央値は使える——これはロバスト推定の核心です。
左の標本平均（コーシー）は試行ごとにまったく異なる軌跡を描き、$n=50000$ でも全く収束する気配がありません。一方、右の標本中央値は5つのトライアルすべてが $n$ の増加とともに着実に真の中央値 $0$ へ収束しています。漸近正規性の理論 $\tilde{X}_n \approx \mathcal{N}(0, \pi^2/(4n))$ が予測するとおり、$n$ が大きいほど散らばりが縮まっており、$n = 10^4$ 以上では $\pm 0.1$ 以内に収まっています。この対比が「コーシー分布の世界で使うべき統計量は中央値」であることを直感的に納得させてくれます。
安定分布族としての位置づけ
コーシー分布が「標本平均が再びコーシー分布になる」という性質を持つのは偶然ではありません。これは 安定分布（stable distribution） という大きな分布族の一員であることから来る性質です。安定分布とは、「独立同分布な確率変数の和が、適当にスケーリングするとまた同じ族に属する分布になる」という性質を満たす分布で、4つのパラメータ $(\alpha, \beta, \gamma, \delta)$ で特徴づけられます。$\alpha \in (0, 2]$ は安定パラメータと呼ばれ、裾の重さを決めます。
- $\alpha = 2$: 正規分布（裾は指数減衰、分散有限）
- $\alpha = 1$: コーシー分布（裾は $1/x^2$、分散も平均も非存在）
- $\alpha = 1/2$: レヴィ分布（片側分布、裾は $1/x^{3/2}$）
一般化中心極限定理（generalized CLT） はこの族を統一して、「重い裾を持つ独立同分布の和は、適切にスケーリングすれば安定分布に分布収束する」と述べます。通常の中心極限定理（正規分布への収束）はその特別ケース（$\alpha = 2$）にすぎず、コーシー分布もまた重要な「普遍極限」の一つなのです。
この理論的事実を特性関数の視点で可視化しておきましょう。
左図では、コーシー分布の特性関数 $e^{-|t|}$（赤）が $t=0$ に尖った頂点を持つのに対し、正規分布の $e^{-t^2/2}$（青）は滑らかな放物線を描きます。この「尖り」が重い裾と対応しています。右図の核心は、$n=3$ の標本平均を考えたとき、コーシーの特性関数（赤実線）が元の分布とまったく同じ形 $e^{-|t|}$ のまま変化しないことです。これが「$\bar{X}_n$ がサンプル数を増やしてもコーシー分布のまま」という事実の本質で、安定分布族に属する証拠が特性関数の不変性として現れています。一方、正規分布（青実線）は $n=3$ で $e^{-t^2/6}$ となり、分散が $1/n$ に縮んでいく（収束する）様子が見えます。
この理論的事実を、次のセクションでロバスト推定の枠組みに昇華させていきます。
ロバスト推定への応用
M推定の枠組み
正規分布の世界では、最尤推定が「残差の二乗和を最小化する」という形をとります。最小二乗法 $\hat\theta = \arg\min_\theta \sum_i (y_i – \mu(\theta; x_i))^2$ は、誤差が正規分布のときの最尤推定なのです。しかし誤差が重い裾を持つとき、二乗は外れ値の影響を爆発的に増幅させてしまいます。残差が10倍ずれたサンプル1つで、二乗和への寄与は100倍。これでは少数の外れ値が推定全体を支配してしまいます。
そこで、二乗の代わりにもっと「優しい」関数 $\rho(\cdot)$ を使った最小化問題
$$ \hat\theta = \arg\min_\theta \sum_i \rho\!\big(r_i(\theta)\big), \quad r_i(\theta) = y_i – \mu(\theta; x_i) $$
を考えます。これを M推定（M-estimator） と呼びます（Mは Maximum likelihood-like の意）。$\rho(r) = r^2$ なら最小二乗法、$\rho(r) = |r|$ なら最小絶対偏差（中央値回帰）、$\rho$ をうまく設計すれば外れ値に頑健な推定が得られます。
主要な損失関数
代表的なロバスト損失関数を3つ挙げます。
- 二乗損失（最小二乗）: $\rho(r) = \tfrac{1}{2} r^2$。導関数（影響関数）$\psi(r) = r$ が無制限に増加するため、外れ値の影響が線形に増える。
- 絶対値損失（$L_1$、最小絶対偏差）: $\rho(r) = |r|$。$\psi(r) = \mathrm{sgn}(r)$ で、外れ値の影響を一定値に抑える。誤差がラプラス分布のときの最尤推定。
- Huber損失: $\rho(r) = \tfrac{1}{2}r^2$（$|r| \leq k$）、$\rho(r) = k(|r| – \tfrac{k}{2})$（$|r| > k$）。$|r| \leq k$ では二乗、$|r| > k$ では絶対値に切り替わる「ハイブリッド」。$\psi(r) = \max(-k, \min(k, r))$ でクリップされる。
- コーシー損失: $\rho(r) = \tfrac{c^2}{2}\log\!\big(1 + (r/c)^2\big)$。$\psi(r) = \tfrac{r}{1 + (r/c)^2}$ で、$|r| \to \infty$ ではむしろゼロに近づく。誤差がコーシー分布のときの最尤推定。
影響関数 $\psi$ の比較
ロバスト性の本質は「1つの外れ値が推定にどれだけ影響を及ぼすか」を表す 影響関数（influence function） $\psi(r) = \rho'(r)$ にあります。各損失の影響関数を比べると、
| 損失 | $\psi(r)$ | $\lvert r \rvert \to \infty$ での挙動 |
|---|---|---|
| 二乗 | $r$ | $\pm\infty$（無制限） |
| 絶対値 | $\mathrm{sgn}(r)$ | $\pm 1$（有界） |
| Huber | $\max(-k, \min(k, r))$ | $\pm k$（有界） |
| コーシー | $r/(1 + (r/c)^2)$ | $0$（redescending） |
コーシー損失の影響関数は「ある程度大きな残差に対しては影響を持つが、極端な残差に対してはむしろ影響を打ち消す」というユニークな振る舞いをします。これは redescending M-estimator（戻り下がり型）と呼ばれ、強い外れ値を実質的に「無視」するメカニズムを内蔵しています。一方、Huberは外れ値の影響を頭打ちにするだけで、完全には無視しません。両者は哲学が異なる損失関数なのです。
左の損失関数グラフでは、二乗損失（青）だけが $|r|$ の増加に対して急激に立ち上がり続け、大きな残差への感度が桁違いに高いことが分かります。HuberとCauchyは遥かに緩やかです。右の影響関数（損失の微分）では、二乗損失 $\psi(r)=r$ が無制限に増加するのに対し、Huberは $\pm k$ でクリップ、コーシーはある値を超えるとむしろ減少して 0 に戻る（redescending） という本質的な違いが一目で確認できます。つまりコーシー損失は「非常に大きな残差のサンプルを推定から実質的に除外する」という自動メカニズムを持っているのです。
反復重み付き最小二乗法（IRLS）
M推定はそのままでは非線形最適化が必要ですが、実用的には 反復重み付き最小二乗法（IRLS: Iteratively Reweighted Least Squares） で解かれます。アイデアは、損失関数を各反復で「等価な重み付き最小二乗問題」に書き換えることです。
$\rho(r)$ の停留条件 $\sum_i \psi(r_i)\, \partial r_i/\partial \theta = 0$ を眺めると、これは重み $w_i = \psi(r_i)/r_i$ を持つ重み付き最小二乗の正規方程式と同じ形をしています。そこで、
- 初期推定 $\hat\theta^{(0)}$（通常はOLS）から残差 $r_i^{(0)}$ を計算
- 重み $w_i^{(t)} = \psi(r_i^{(t)})/r_i^{(t)}$ を計算
- 重み付き最小二乗 $\hat\theta^{(t+1)} = (X^\top W^{(t)} X)^{-1} X^\top W^{(t)} y$ で更新
- 収束するまで繰り返す
という反復計算で解けます。各反復で外れ値（残差が大きいサンプル）の重みが小さくなり、推定が頑健になっていきます。コーシー損失なら $w_i = 1/(1 + (r_i/c)^2)$、Huber損失なら $w_i = \min(1, k/|r_i|)$ という形になります。
コーシー分布の他の応用例
ロバスト回帰以外にも、コーシー分布は実用上の多様な場面で活躍します。
- MCMCの提案分布: メトロポリス・ヘイスティングス法で「広い裾」を持つ提案分布としてコーシー分布が使われることがあります。重い裾のおかげで、目標分布の多峰性に対しても遠くのモードに飛び移れる確率が高まります。
- 画像処理の中央値フィルタ: ソルト＆ペッパーノイズのような重い裾のノイズに対し、平均化フィルタが惨敗するのに対し、中央値フィルタは見事にノイズを除去します。これは中央値が「コーシー型ノイズに対する最尤推定」であるという理論に裏付けられています。
- 金融工学のリスク評価: 株価リターンは正規分布よりはるかに重い裾を持ち、コーシー分布や安定分布のフィッティングが使われます。VaR（Value at Risk）の正確な見積りに必要です。
- 物理学の共鳴フィット: 加速器実験でのピークフィッティングは、ほぼ常にローレンツ／Breit-Wigner関数（=コーシー分布）と検出器分解能（=ガウス分布）の畳み込みである フォークト関数 を使います。
理論はここまでです。M推定の実装と、コーシー分布の奇妙な性質をPythonで実際に目で見て確かめましょう。
Pythonでの可視化とロバスト回帰
確率密度と裾の比較
まず、標準コーシー分布と標準正規分布を重ねて、裾の重さの違いを直接見てみます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# x軸の範囲（裾の重さを見るため広く取る）
x = np.linspace(-10, 10, 1000)
# 確率密度
pdf_normal = stats.norm.pdf(x)
pdf_cauchy = stats.cauchy.pdf(x)
# 線形スケールと対数スケールで2枚並べる
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].plot(x, pdf_normal, 'b-', lw=2, label='Standard normal')
axes[0].plot(x, pdf_cauchy, 'r-', lw=2, label='Standard Cauchy')
axes[0].set_xlabel('x')
axes[0].set_ylabel('f(x)')
axes[0].set_title('Linear scale')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].semilogy(x, pdf_normal, 'b-', lw=2, label='Standard normal')
axes[1].semilogy(x, pdf_cauchy, 'r-', lw=2, label='Standard Cauchy')
axes[1].set_xlabel('x')
axes[1].set_ylabel('f(x) (log scale)')
axes[1].set_title('Log scale — tail comparison')
axes[1].set_ylim(1e-20, 1)
axes[1].legend()
axes[1].grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('cauchy_pdf_compare.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、2つの分布の決定的な違いが鮮明に読み取れます。線形スケールでは、コーシー分布の方が中心ピークがやや低く、肩が広く、裾が太い、という違いがわかります。しかし真の衝撃は対数スケールにあります。正規分布は $|x| = 10$ で確率密度が $10^{-23}$ 程度まで真っ逆さまに落ちていますが、コーシー分布は $|x| = 10$ でもまだ $10^{-3}$ 近くの値を保っています。実に 20桁 の差です。これが「$|x|$ が大きい値が日常的にやってくる」という直感の数学的根拠です。
標本平均が収束しないMonte Carlo
次に、コーシー分布のサンプルから累積平均を計算し、正規分布のそれと比較します。大数の法則が成り立たない様子を直接目で見るのが目的です。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
N = 100000 # 標本サイズ
# 5つのトライアル
n_trials = 5
fig, axes = plt.subplots(1, 2, figsize=(13, 5), sharey=False)
# 正規分布の累積平均
for k in range(n_trials):
samples = stats.norm.rvs(size=N, random_state=k + 100)
cum_mean = np.cumsum(samples) / np.arange(1, N + 1)
axes[0].plot(cum_mean, alpha=0.7, lw=0.8, label=f'trial {k+1}')
axes[0].axhline(0, color='black', ls='--', alpha=0.5, label='true mean = 0')
axes[0].set_xscale('log')
axes[0].set_xlabel('n (log scale)')
axes[0].set_ylabel('cumulative mean')
axes[0].set_title('Normal: cumulative mean converges to 0 (LLN holds)')
axes[0].set_ylim(-2, 2)
axes[0].legend(loc='upper right', fontsize=8)
axes[0].grid(True, which='both', alpha=0.3)
# コーシー分布の累積平均
for k in range(n_trials):
samples = stats.cauchy.rvs(size=N, random_state=k + 200)
cum_mean = np.cumsum(samples) / np.arange(1, N + 1)
axes[1].plot(cum_mean, alpha=0.7, lw=0.8, label=f'trial {k+1}')
axes[1].axhline(0, color='black', ls='--', alpha=0.5, label='naive "mean" = 0?')
axes[1].set_xscale('log')
axes[1].set_xlabel('n (log scale)')
axes[1].set_ylabel('cumulative mean')
axes[1].set_title('Cauchy: cumulative mean does NOT converge')
axes[1].legend(loc='upper right', fontsize=8)
axes[1].grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('cauchy_lln_fail.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフは大数の法則の崩壊をこれ以上ないほど雄弁に物語っています。左の正規分布では、5つのトライアルすべてが $n$ の増加とともにきれいに $0$ に収束していき、$n = 10^5$ ではすべての軌跡が $\pm 0.01$ の幅に収まっています。これが大数の法則の威力です。一方、右のコーシー分布は全く異なる景色を見せます。累積平均は $n$ をいくら増やしても安定せず、サンプル列の途中でたまに巨大な値が紛れ込むとジャンプします。$n = 10^5$ でも平均が $\pm 1$ や $\pm 10$ のスケールで暴れています。理論で示した「$\bar X_n$ はまた標準コーシー分布」という事実が、目に見える形で立ち現れているのです。
標本中央値は収束する
「平均は壊れる、では何を使えばよいか」の答えとして、同じシミュレーションを中央値で行います。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(123)
N = 100000
n_trials = 5
# n軸を間引いて高速化（全n点での中央値計算は重いため）
ns_eval = np.unique(np.round(np.logspace(0, np.log10(N), 200)).astype(int))
fig, ax = plt.subplots(figsize=(10, 5))
for k in range(n_trials):
samples = stats.cauchy.rvs(size=N, random_state=k + 300)
med = np.array([np.median(samples[:n]) for n in ns_eval])
ax.plot(ns_eval, med, alpha=0.7, lw=1, label=f'trial {k+1}')
ax.axhline(0, color='black', ls='--', alpha=0.5, label='true median = 0')
ax.set_xscale('log')
ax.set_xlabel('n (log scale)')
ax.set_ylabel('running median')
ax.set_title('Cauchy: sample MEDIAN converges (unlike sample mean)')
ax.set_ylim(-1, 1)
ax.legend(loc='upper right', fontsize=8)
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('cauchy_median_converges.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフは、前のグラフとは打って変わって美しい収束を見せます。5つのトライアルすべてが $n$ の増加とともに $0$（真の中央値 $x_0 = 0$）に収束していき、$n = 10^5$ では $\pm 0.02$ 程度の幅に収まっています。理論が予言した中央値の漸近正規性 $\tilde X_n \sim \mathcal N(0, \pi^2/(4n))$ では、$n = 10^5$ で標準偏差が約 $\pi/(2\sqrt{10^5}) \approx 0.005$ ですから、実測の散らばりとオーダーで整合しています。コーシー分布の真の位置を推定したければ、平均ではなく中央値を使え——これがロバスト統計の最も基本的な処方箋です。
Huber損失とコーシー損失の比較
ロバスト回帰の理論を実装する前に、二乗損失・Huber損失・コーシー損失の形状と影響関数を可視化しておきます。
import numpy as np
import matplotlib.pyplot as plt
r = np.linspace(-6, 6, 400)
k_huber = 1.345 # Huberの標準的なk（正規誤差で95%効率）
c_cauchy = 2.385 # コーシー損失の標準的なc
# 各損失関数
rho_l2 = 0.5 * r**2
rho_huber = np.where(np.abs(r) <= k_huber,
0.5 * r**2,
k_huber * (np.abs(r) - 0.5 * k_huber))
rho_cauchy = 0.5 * c_cauchy**2 * np.log(1 + (r / c_cauchy)**2)
# 影響関数（ρの導関数）
psi_l2 = r
psi_huber = np.clip(r, -k_huber, k_huber)
psi_cauchy = r / (1 + (r / c_cauchy)**2)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].plot(r, rho_l2, 'b-', lw=2, label='L2 (squared)')
axes[0].plot(r, rho_huber, 'g-', lw=2, label=f'Huber (k={k_huber})')
axes[0].plot(r, rho_cauchy, 'r-', lw=2, label=f'Cauchy (c={c_cauchy})')
axes[0].set_xlabel('residual r')
axes[0].set_ylabel('rho(r)')
axes[0].set_title('Loss function')
axes[0].set_ylim(0, 12)
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].plot(r, psi_l2, 'b-', lw=2, label='L2 (squared)')
axes[1].plot(r, psi_huber, 'g-', lw=2, label=f'Huber (k={k_huber})')
axes[1].plot(r, psi_cauchy, 'r-', lw=2, label=f'Cauchy (c={c_cauchy})')
axes[1].axhline(0, color='black', alpha=0.3)
axes[1].set_xlabel('residual r')
axes[1].set_ylabel('psi(r) = drho/dr')
axes[1].set_title('Influence function')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('robust_loss_compare.png', dpi=150, bbox_inches='tight')
plt.show()
左の損失関数のグラフから、二乗損失だけが残差の絶対値が大きいところで急激に立ち上がっているのが見えます。HuberとCauchyは大きな残差に対して緩やかにしか増えません。特にコーシー損失は対数オーダーでしか増加しないため、極端な残差にほとんど反応しません。右の影響関数を見ると違いがさらに明瞭です。二乗損失の $\psi(r) = r$ は際限なく増加し、Huberは $\pm k$ で頭打ちになります。最も特徴的なのがコーシー損失で、$|r|$ がある程度大きくなると $\psi$ が逆に減少して0に近づく——これがredescending M-estimatorの正体です。極端な外れ値は、コーシー損失にとっては「存在しないも同然」になるのです。
外れ値混入データへのロバスト線形回帰（IRLS実装）
最後に、外れ値を含む線形回帰データに対して、OLS／Huber／Cauchyの3つの推定を比較します。IRLSをスクラッチで実装します。
import numpy as np
def irls_regression(X, y, weight_fn, max_iter=50, tol=1e-6):
"""反復重み付き最小二乗法によるロバスト回帰
weight_fn(r, scale) は残差とスケールから重みを返す関数"""
n, p = X.shape
# 初期値はOLS
beta = np.linalg.lstsq(X, y, rcond=None)[0]
for it in range(max_iter):
r = y - X @ beta
# スケール推定はMAD（残差の中央絶対偏差を 0.6745 で割って σ 換算）
scale = np.median(np.abs(r - np.median(r))) / 0.6745
if scale < 1e-12:
scale = 1e-12
w = weight_fn(r, scale)
# 重み付き最小二乗で更新
W = np.diag(w)
beta_new = np.linalg.solve(X.T @ W @ X, X.T @ W @ y)
if np.linalg.norm(beta_new - beta) < tol * (np.linalg.norm(beta) + 1e-12):
beta = beta_new
break
beta = beta_new
return beta, it + 1
# Huberの重み関数（w_i = min(1, k/|z|)、z=r/scale）
def huber_weight(r, scale, k=1.345):
z = r / scale
return np.where(np.abs(z) <= k, 1.0, k / np.maximum(np.abs(z), 1e-12))
# コーシーの重み関数（w_i = 1/(1 + (z/c)^2)）
def cauchy_weight(r, scale, c=2.385):
z = r / scale
return 1.0 / (1.0 + (z / c)**2)
ここでスケール推定にMAD（残差の中央絶対偏差）を $0.6745$ で割ったものを使っているのは、正規分布のときにこの量がちょうど標準偏差の不偏推定量になるからです（$\mathrm{MAD}/0.6745 \approx \sigma$）。MADは中央値ベースの統計量なので、外れ値に頑健にスケールを推定できます。これがないと、外れ値がスケールを大きく見積もらせ、その結果すべてのサンプルが「内側」に見えてしまうという破綻が起こります。
import numpy as np
import matplotlib.pyplot as plt
# テストデータ: 真のモデル y = 2x + 1 + ノイズ + 外れ値
np.random.seed(7)
n = 100
x = np.linspace(0, 10, n)
y_true = 2.0 * x + 1.0
y = y_true + np.random.normal(0, 0.5, n)
# 全体の20%を強い外れ値で汚染
n_outlier = 20
outlier_idx = np.random.choice(n, n_outlier, replace=False)
y[outlier_idx] += np.random.normal(0, 20, n_outlier) + 15
# 計画行列
X = np.column_stack([np.ones(n), x])
# 3つの推定
beta_ols = np.linalg.lstsq(X, y, rcond=None)[0]
beta_huber, it_h = irls_regression(X, y, huber_weight)
beta_cauchy, it_c = irls_regression(X, y, cauchy_weight)
print(f"True params: [intercept, slope] = [1.0, 2.0]")
print(f"OLS: {beta_ols}")
print(f"Huber ({it_h} iters): {beta_huber}")
print(f"Cauchy ({it_c} iters): {beta_cauchy}")
# 可視化
xs = np.linspace(0, 10, 100)
fig, ax = plt.subplots(figsize=(10, 6))
clean_idx = np.setdiff1d(np.arange(n), outlier_idx)
ax.scatter(x[clean_idx], y[clean_idx], c='blue', s=30, alpha=0.6, label='clean')
ax.scatter(x[outlier_idx], y[outlier_idx], c='red', s=60, marker='x',
label='outlier (20%)')
ax.plot(xs, 2 * xs + 1, 'k--', lw=2, label='True (y=2x+1)')
ax.plot(xs, beta_ols[0] + beta_ols[1] * xs, 'orange', lw=2,
label=f'OLS (slope={beta_ols[1]:.2f})')
ax.plot(xs, beta_huber[0] + beta_huber[1] * xs, 'g-', lw=2,
label=f'Huber (slope={beta_huber[1]:.2f})')
ax.plot(xs, beta_cauchy[0] + beta_cauchy[1] * xs, 'm-', lw=2,
label=f'Cauchy (slope={beta_cauchy[1]:.2f})')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Robust regression: OLS vs Huber vs Cauchy (20% outliers)')
ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('robust_regression.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、3手法のロバスト性の差が一目瞭然に読み取れます。OLS（オレンジ）は外れ値（赤×）に強く引っ張られ、真の直線（黒破線、$y = 2x+1$）から明らかにずれた傾きを推定しています。一方、Huber（緑）とCauchy（マゼンタ）は外れ値をほぼ無視して、真の直線にぴたりと重なっています。出力された数値でも、OLSの傾きが真値2.0から大きくずれているのに対し、Huber/Cauchyは小数第2位まで真値に一致します。特にコーシー損失は影響関数がredescendingなため、極端な外れ値を完全に「忘れる」ことができ、より頑健な推定が得られています。これがM推定とIRLSの実力です。
実測値として OLS傾き≈1.66、Huber傾き≈2.00、コーシー損失傾き≈2.01 が得られており、真値2.0と比べてOLSは外れ値に引きずられて大きくずれているのに対し、HuberとCauchyはともに真値にほぼ一致しています。外れ値（赤×）が $y \approx 30〜60$ 付近に集中しているにもかかわらず、緑と紫の回帰線はクリーンなデータ点の傾向を正確にとらえています。HuberとCauchyのどちらも外れ値20%程度には十分頑健で、損失関数によってIRLSが自動的に外れ値の重みを絞ることでOLSの崩壊を回避できることが数値実証されています。
外れ値割合に対する頑健性のサーベイ
最後に、外れ値の割合を変えながら、各手法の推定精度を比較します。各手法が「どのくらい汚染に耐えられるか」を定量的に評価します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
n = 200
true_slope, true_intercept = 2.0, 1.0
x = np.linspace(0, 10, n)
X = np.column_stack([np.ones(n), x])
# 汚染割合を変える
contamination_rates = np.linspace(0, 0.5, 11)
n_repeat = 30 # ノイズ実現を変えて平均化
results = {'OLS': [], 'Huber': [], 'Cauchy': []}
for rate in contamination_rates:
errs = {'OLS': [], 'Huber': [], 'Cauchy': []}
for rep in range(n_repeat):
rng = np.random.default_rng(rep + int(rate * 1000))
y = true_intercept + true_slope * x + rng.normal(0, 0.5, n)
n_out = int(n * rate)
if n_out > 0:
idx = rng.choice(n, n_out, replace=False)
y[idx] += rng.normal(0, 30, n_out) + rng.choice([-30, 30], n_out)
beta_ols = np.linalg.lstsq(X, y, rcond=None)[0]
beta_huber, _ = irls_regression(X, y, huber_weight)
beta_cauchy, _ = irls_regression(X, y, cauchy_weight)
# 傾きの誤差
errs['OLS'].append(abs(beta_ols[1] - true_slope))
errs['Huber'].append(abs(beta_huber[1] - true_slope))
errs['Cauchy'].append(abs(beta_cauchy[1] - true_slope))
for k in results:
results[k].append(np.mean(errs[k]))
fig, ax = plt.subplots(figsize=(10, 5.5))
ax.plot(contamination_rates * 100, results['OLS'], 'o-', lw=2,
label='OLS', color='orange')
ax.plot(contamination_rates * 100, results['Huber'], 's-', lw=2,
label='Huber', color='green')
ax.plot(contamination_rates * 100, results['Cauchy'], '^-', lw=2,
label='Cauchy', color='magenta')
ax.set_xlabel('Outlier contamination rate [%]')
ax.set_ylabel('|estimated slope - true slope|')
ax.set_title('Slope estimation error vs contamination rate')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('robust_breakdown.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、3手法の頑健性の階層が明瞭に見えてきます。OLSは外れ値が増えるにつれて誤差が線形に増大し、汚染率20%ではすでに傾きの誤差が大きく崩れています。これが「OLSのbreakdown point は0」と言われる所以で、たった1サンプルの極端な外れ値でも理論上推定が無限大にずれうるのです。Huberは汚染率20〜30%程度までは持ちこたえますが、それを超えると徐々に崩れていきます。最も頑健なのがコーシー損失で、汚染率40%でも傾きを真値近くに保っています。redescendingな影響関数のおかげで、多くの外れ値が「重み0」相当として扱われているからです。実用では、外れ値が少ないと予想される場合はHuber（効率的）、多い／不明な場合はCauchy（頑健）と使い分けるのが定石です。
OLS（オレンジ）は汚染率0%でほぼ誤差ゼロからスタートするものの、汚染率5%を超えたあたりから急激に誤差が増え、40%では真値2.0からのずれが数倍に膨らんでいます。Huber（緑）は20%程度まで低誤差を維持しますが、それを超えると徐々に崩れます。コーシー損失（紫）は汚染率40〜50%の領域でも誤差が抑制されており、3手法の中で最も頑健であることが数値実験で確かめられました。外れ値の割合が不明な実務状況ではコーシー損失を選ぶ合理的な根拠がここにあります。
まとめ
本記事では、コーシー分布の数学的性質とロバスト統計への応用を、理論の導出からPython実装まで一気通貫で解説しました。
- コーシー分布の密度: $f(x; x_0, \gamma) = \dfrac{1}{\pi\gamma[1 + ((x-x_0)/\gamma)^2]}$。一様分布する角度の正接変換から自然に導出され、また2つの独立な標準正規の比 $U/V$ としても現れる。ローレンツ分布／Breit-Wigner分布／自由度1のt分布として複数の顔を持つ。
- 平均・分散が存在しない: $\int |x| f(x)\, dx = \infty$ となり、絶対可積分性が破綻するため期待値が定義できない。$E[X^2]$ も発散し分散も存在しない。一方、中央値とMAD は美しく定義され、それぞれ位置パラメータ $x_0$ と尺度パラメータ $\gamma$ に一致する。
- 標本平均が収束しない: 独立な標準コーシー和は再びコーシー分布、$\bar X_n$ は何度平均しても標準コーシー分布のまま。大数の法則も中心極限定理も適用できない。標本中央値は漸近的に $\mathcal N(x_0, \pi^2\gamma^2/(4n))$ に収束する。
- 安定分布族: コーシーは $\alpha = 1$ の安定分布で、$\alpha = 2$ の正規分布と並ぶ「もう一つの普遍極限」。一般化中心極限定理が、重い裾を持つ和の極限としてガウス分布の代わりに安定分布を与える。
- ロバスト推定: 二乗損失の代わりに $\rho(r)$ を最小化するM推定。Huber損失は影響関数を頭打ちにし、コーシー損失は redescending で極端な外れ値を実質的に無視する。IRLSで効率的に解ける。
- 実装での威力: 外れ値が20〜40%混入したデータでも、Huber／Cauchyは真の傾きを正確に推定。一方OLSは少数の外れ値で大きく崩れる。
コーシー分布が示すのは、「中心極限定理が万能ではない」という確率論の地平の広がりです。データに重い裾が疑われるとき、平均と分散をうかつに計算してはいけない——中央値・MAD・ロバスト推定を使うべきだ、という実践的教訓は、計測ノイズ・市場リターン・画像ノイズ・パルス性スパイクノイズといったあらゆる「汚いデータ」の場面で生きます。
次のステップとして、以下の記事も参考にしてください。
