🎨 直感で掴む — 「東京問題」を例に

日本の都道府県データを扱うと、 ほぼ必ず 「東京問題」に直面します。 県内総生産（GDP）でも、 人口でも、 大学進学率でも、 東京は他 46 県と桁が違うレベルで大きい。 47 県の 普通の平均 を計算すると：

SSDSE-B 県内総生産（2023年）の例

「普通の県の県内総生産はだいたいいくら？」と問われたとき、 12.3 兆円 と答えるのと 7.8 兆円 と答えるのとでは、 意味がまったく違います。 平均は「東京を含めた全体の総量を 47 で割った」値であり、 「平均的な県の像」ではありません。 中央値やトリム平均こそが「平均的な県」を表す指標です。

🎢 「破綻点」の感覚

ロバスト統計の核心概念が breakdown point（破綻点、 BP） です。 BP とは 「データの何 % が汚染されたら、 推定量が無限大（または不定）に発散するか」。 直感的には「どこまでの汚染に耐えられるか」のスコア：

• 平均：BP = 0%。 1 点でも +∞ にあれば平均も +∞ に発散。 最弱
• 標準偏差：BP = 0%。 同上
• 5% トリム平均：BP = 5%。 5% 以下なら吸収できる
• 10% トリム平均：BP = 10%
• Huber M-estimator：BP ≈ 5〜30%（チューニング次第）
• 中央値：BP = 50%。 理論的最大値。 半分のデータが壊れても生存
• MAD（中央絶対偏差）：BP = 50%。 中央値と並ぶ「最も頑健」な推定量

逆に、 BP = 50% を超える推定量は原理的に作れません（半分以上が「壊れた値」なら、 もはやそれが多数派で真の値）。 中央値は最強です。

📐 数式 — 代表的なロバスト推定量

① 中央値 と MAD

② トリム平均（trimmed mean）

③ Winsorized 平均（端を切り詰める）

④ Huber M-estimator（最尤系の頑健化）

M-estimator は、 損失関数 $\rho(\cdot)$ の和を最小化する位置パラメータ：

Huber の $\rho$ 関数は、 中心付近では二乗（最小二乗の効率性を保つ）、 外れ値領域では絶対値（線形にしか伸びない）：

⑤ Tukey の biweight（再降下型）

⑥ ロバスト回帰：LMS と LTS

🔬 数式を「言葉」で読み解く

🧮 実値で計算してみる — SSDSE-B 県内総生産

STEP 1：47 県の県内総生産（兆円）の主要値

SSDSE-B-2026 から、 2023 年の都道府県別の県内総生産（実値、 兆円）を抜粋：

※ 実際の中央値順位の県は年度・統計値により変動するため概数。 計算は実 CSV で確認のこと。

STEP 2：算術平均 vs 中央値 vs ロバスト

STEP 3：MAD で散らばりを測る

標準偏差は $s \approx$ 17.5 兆円（東京の影響で巨大）。 一方 MAD（中央絶対偏差）は：

STEP 4：解釈の違い

• 「日本の都道府県の経済規模はどれくらい？」と問われ、 平均 12.3 兆円と答えるのは正確だが誤解を招く。 47 県のうち 10 県程度しか 12 兆円を超えない
• 中央値 7.8 兆円なら、 「23 県が 7.8 兆円以下、 23 県が以上」と正確な分布の中央を表す
• 政策立案で「平均的な県の支援策」を考えるなら、 平均ではなく中央値を基準にすべき

📈 効率（efficiency）と頑健性（robustness）のトレードオフ

ロバスト統計の中心的な議論は 「効率」と「頑健性」のトレードオフです。 効率とは「正規分布下で、 推定量の分散がどれだけ小さいか（＝最尤推定との比較）」、 頑健性とは「外れ値や分布のズレに対して推定量がどれだけ揺らがないか」。 両者は逆方向に動くため、 折衷案が必要です。

選び方のフロー：(1) 外れ値の比率が不明・大きい → 中央値 / MAD、 (2) 外れ値はあるが少数 → Huber、 (3) BP を最大化したい → biweight / LMS / LTS、 (4) 多変量で BP も効率も両立 → MM-estimator。

影響関数（Influence Function, IF）の比較

影響関数は「1 点の汚染が推定量にどれだけ影響するか」を示す関数。 視覚的には「推定量の局所感度」のグラフです：

• 算術平均：IF($z$) = $z - \mu$ → 線形に無限大まで伸びる。 z = +∞ なら影響も +∞。 致命的弱点
• 中央値：IF($z$) = sign($z - \mu$) / (2$f(\mu)$) → ±定数で有界。 どれだけ極端な z でも影響は ±一定
• Huber M (k=1.345)：IF が中心は線形、 |z|> k で flat（±k）になる → 連続有界
• Tukey biweight：IF が中心はゼロから線形に上がり、 |z|>c でゼロに戻る（再降下）。 完全に外れ値を無視できる

影響関数が「無界 → 有界 → ゼロに収束」と進むにつれ、 頑健性は高まりますが、 効率は徐々に低下します。 Huber は「ほとんど効率を犠牲にせず、 IF を有界化」した秀作。

📊 外れ値混入シミュレーション — 推定量の挙動を可視化

正規分布 N(0, 1) からの n = 100 個のデータに、 「10% を外れ値 N(20, 1)」として混入させたとき、 各推定量の動きを比較します。

観察：算術平均は汚染比率に線形に動く（致命的）。 中央値は 50% を超えるまでほぼ動かない。 5% トリム平均は BP=5% を超えると破綻。 Huber は連続的に劣化するが、 50% 近くまで耐える。 「BP を超えたら推定量はもはや信用ならない」という事実が体感できる例。

収益データでの応用例

金融分野では、 「日次収益の平均は外れ値（暴落・暴騰）に致命的に弱い」ことが知られており、 中央値や Huber 推定量での要約が業界標準です：

• 日経平均日次収益の単純平均（10 年） ≈ 0.04% / 日 ← リーマンや 2020/3 の暴落で歪む
• 同じデータの中央値 ≈ 0.07% / 日 ← 「典型的な 1 日」を表す
• Huber 平均 ≈ 0.06% / 日 ← 中央値と平均の折衷

「典型的な収益」を語るなら中央値・Huber、 「長期累積収益」を語るなら算術平均、 と用途で使い分けるのが正解。

🐍 Python 実装 — 4 通りのロバスト推定量を比較

方法 A：手早く scipy で（推奨）

# SSDSE-B 県内総生産での ロバスト推定量の比較
import numpy as np
import pandas as pd
from scipy import stats

# データ読込
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=2)
df = df.iloc[:, [0, 2, 31]].copy()  # 年度、都道府県、県内総生産（B4101 列）
df.columns = ['year', 'pref', 'gdp']
df = df[df['year'] == 2023].dropna()
x = df['gdp'].values / 1e6  # 兆円単位に換算（元データは百万円）

print(f'件数 n = {len(x)}')
print(f'算術平均: {np.mean(x):.2f}')
print(f'中央値: {np.median(x):.2f}')
print(f'5% トリム平均: {stats.trim_mean(x, 0.05):.2f}')
print(f'10% トリム平均: {stats.trim_mean(x, 0.10):.2f}')
print(f'標準偏差: {np.std(x, ddof=1):.2f}')
print(f'MAD (raw): {stats.median_abs_deviation(x):.2f}')
print(f'MAD scaled (~σ): {stats.median_abs_deviation(x, scale="normal"):.2f}')
print(f'IQR: {stats.iqr(x):.2f}')

方法 B：statsmodels で Huber M-estimator

from statsmodels.robust.scale import Huber
from statsmodels.robust.norms import HuberT
import numpy as np

# Huber M-estimator（位置パラメータ）
huber = Huber(c=1.345)
mu_huber, s_huber = huber(x)
print(f'Huber μ̂ = {mu_huber:.2f} 兆円  (s = {s_huber:.2f})')
print(f'比較：算術平均 {np.mean(x):.2f}, 中央値 {np.median(x):.2f}')

# Tukey biweight も試す
from statsmodels.robust.norms import TukeyBiweight
# biweight は scale=Hampel など別途、 詳細は公式 doc

方法 C：ロバスト回帰（県内総生産 vs 人口）

import numpy as np
import pandas as pd
import statsmodels.api as sm

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=2)
df = df.iloc[:, [0, 2, 3, 31]].copy()
df.columns = ['year', 'pref', 'pop', 'gdp']
df = df[df['year'] == 2023].dropna()

x = df['pop'].values  # 千人
y = df['gdp'].values / 1e6  # 兆円
X = sm.add_constant(x)

# 通常の OLS（東京で歪む）
ols = sm.OLS(y, X).fit()
print(f'OLS:    intercept={ols.params[0]:.3f}, slope={ols.params[1]:.6f}')

# Huber 回帰（IRLS, iteratively reweighted least squares）
huber = sm.RLM(y, X, M=sm.robust.norms.HuberT(t=1.345)).fit()
print(f'Huber:  intercept={huber.params[0]:.3f}, slope={huber.params[1]:.6f}')

# Tukey biweight 回帰
biw = sm.RLM(y, X, M=sm.robust.norms.TukeyBiweight()).fit()
print(f'Tukey:  intercept={biw.params[0]:.3f}, slope={biw.params[1]:.6f}')

方法 D：sklearn HuberRegressor

from sklearn.linear_model import HuberRegressor, LinearRegression
import numpy as np

x = df['pop'].values.reshape(-1, 1)
y = df['gdp'].values / 1e6

ols = LinearRegression().fit(x, y)
huber = HuberRegressor(epsilon=1.35).fit(x, y)

print(f'OLS 傾き:   {ols.coef_[0]:.6f}')
print(f'Huber 傾き: {huber.coef_[0]:.6f}')
print(f'Huber が外れ値扱いした県数: {(~huber.outliers_).sum()} 正常 / {huber.outliers_.sum()} 外れ')

方法 E：可視化（外れ値の影響を視覚化）

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(x, y, alpha=0.7, s=40)
xs = np.linspace(x.min(), x.max(), 200).reshape(-1, 1)
ax.plot(xs, ols.predict(xs), 'r-', lw=2, label=f'OLS（東京に引っ張られる）')
ax.plot(xs, huber.predict(xs), 'b-', lw=2, label=f'Huber（東京を緩和）')
ax.set_xlabel('総人口（千人）'); ax.set_ylabel('県内総生産（兆円）')
ax.legend(); ax.set_title('47都道府県の人口 vs 県内総生産：OLS vs Huber')
plt.tight_layout(); plt.savefig('robust_regression.png', dpi=120)

⚠️ ロバスト統計の 5 つの落とし穴

🏛 歴史的背景 — Tukey から Huber へ

ロバスト統計の起源は、 1953 年の John W. Tukey のテクニカルレポート "The Future of Data Analysis" にさかのぼります。 Tukey は「データ解析は確率論や数理統計から独立した実証科学である」と宣言し、 探索的データ解析（EDA）の創始者として、 平均ではなく中央値・四分位を主要な記述統計として推奨しました。

その後 Peter J. Huber が 1964 年の "Robust Estimation of a Location Parameter"（Annals of Mathematical Statistics）で数学的厳密性を伴う理論体系を構築。 M-estimator、 効率の最小化、 minimax 戦略といった現代ロバスト統計の中核概念がここで定式化されました。

主要マイルストーン

日本における普及は 1980 年代後半から、 Wilcoxon 検定や中央値の使用を通じて医学・心理学の領域で広まり、 現在では「外れ値が懸念される場面ではロバスト統計を併用する」が学会論文の標準作法となっています。

🌐 ロバスト統計の派生・関連手法

ロバスト統計は単一の手法ではなく、 分野横断のフレームワークです。 代表的な道具を分類して整理：

💼 実務での応用シナリオ

① 都道府県データ・市区町村データ（東京問題）

本サイトの SSDSE-B のような、 少数の超巨大行政単位（東京・大阪・神奈川・愛知）と多数の中小行政単位が混在するデータでは、 単純平均が東京に強く引きずられます。 こうした場合、 中央値・トリム平均・対数変換 + 平均（GDP のような桁差データ）を併用し、 「どの代表値で見ても同じ結論か」をチェックするのが基本です。

② 賃金・所得データ（右に長い裾）

個人所得は典型的な右裾分布（log-normal に近い）で、 上位 1% の高所得者が平均を引き上げます。 「平均年収 700 万円」は実態と乖離する典型例。 中央値所得 400 万円のほうが「中間層の典型」を示し、 政策議論ではこちらが標準。

③ レスポンスタイム・処理時間（重い裾）

Web サービスのレスポンスタイム、 API 呼び出し時間、 トランザクション処理時間など、 「ほとんどは 10ms 以下だが、 たまに 5 秒かかる」というデータでは、 単純平均は数秒単位に跳ね上がります。 業界標準は 中央値（p50）、 p95、 p99 のパーセンタイルで報告する SRE プラクティス。

④ センサーデータ・IoT

センサー誤動作で時折ありえない値（−999、 +99999）が混入。 リアルタイム監視では、 ストリーミング中央値（streaming median）や Huber フィルタを使い、 「外れ値を自動的に減衰させつつ平均トレンドを追う」のが定番。

⑤ アンケート調査

5 段階リッカート尺度の集計で、 「いたずら回答」「セット回答（全部 3 を選ぶ）」が混入。 平均評価より中央値・トリム平均で集計、 また外れ回答を Mahalanobis 距離などで検出して除外するのが先進的アプローチ。

⑥ 機械学習の損失関数

回帰モデルの損失関数として、 MSE は外れ値に二乗のペナルティで弱い。 Huber 損失は中心は二乗（収束が速い）、 外側は線形（外れ値に頑健）。 sklearn の HuberRegressor、 PyTorch の nn.SmoothL1Loss として実装され、 物体検出（Fast R-CNN）でも採用。

⑦ A/B テストの効果推定

クリック率や購入金額の A/B テストで、 まれな大口購入者が平均を歪める。 中央値での比較、 順位検定（Mann-Whitney U）、 Bootstrap + トリム平均などが頑健化策。 Netflix や Airbnb の A/B テストプラットフォームではロバスト推定が組込み済み。

🔗 関連用語（前提・並列・発展）

📚 前提（先に押さえたい）

🤝 並列（同じレイヤーの仲間）

🚀 発展（その先へ）

🧪 古典統計 vs ロバスト統計 — 早見表

運用ルール：論文・レポートを書くときは、 古典推定量とロバスト推定量を両方計算し、 「両者で結論が一致するか」を必ず確認。 一致すれば結論は強固。 大きく異なれば、 データの中に外れ値や分布の歪みが存在し、 そこに重要な情報があるサインです。

🛠 トラブルシューティング集

Q1：平均と中央値の値が大きく違う。 どちらを報告すべき？

回答：両方を併記するのが正解。 「平均は X、 中央値は Y、 差の原因は上位 10% に集中している大口値」のように説明文を付ける。 読者が用途に応じて使い分けられる。 平均を報告するなら必ず標準偏差も併記し、 中央値を報告するなら MAD または IQR を併記。

Q2：Huber と Tukey biweight、 どちらを使うべき？

回答：「Huber でまず試す」が原則。 Huber は連続的に重み付けし、 効率と頑健性のバランスが良く、 IRLS で安定収束する。 Tukey biweight は外れ値を完全にゼロ重みにする再降下型なので、 BP が高いが収束が難しいことがある。 多変量や複雑な状況では MM-estimator（biweight + Huber 初期値）が現代的な選択肢。

Q3：ロバスト回帰で SE をどう報告すべき？

回答：解析的な SE 公式は手法ごとに異なり複雑。 もっとも汎用な解は Bootstrap で SE を推定。 scipy.stats.bootstrap や sklearn.utils.resample で 1,000〜10,000 回再標本化して、 推定量の分布から CI を出す。

Q4：n が小さい（n < 20）。 ロバスト統計は使える？

回答：使えるが効果は限定的。 n=10 で「外れ値の比率」を見積もるのが困難で、 BP 50% を保証する手法は事実上 LMS や中央値だけ。 n が小さいときは「外れ値の候補を 1 つずつ手で確認」「Bootstrap で SE を推定」「複数推定量で結論の一貫性をチェック」が王道。

Q5：外れ値を除外したほうが結果がきれいになる。 除外していい？

回答：基本的に NG。 外れ値の原因を調査するのが先。 (a) 測定ミス・データ入力ミスなら除外可（ただしその旨を論文に明記）、 (b) 真の極端値（東京、 GAFAM 等）なら除外せずに対数変換、 ロバスト推定、 階層モデルなどで対応。 「結果がきれいになるから除外」は p-hacking と同じ問題行動。

Q6：機械学習で外れ値があるとき、 損失関数として何を選ぶ？

回答：回帰タスクなら Huber 損失（PyTorch の SmoothL1Loss、 sklearn の HuberRegressor）が定番。 中心は二乗で学習が速く、 外れ値領域は線形で頑健。 quantile regression（pinball loss）でメディアン回帰も選択肢。 物体検出の Fast R-CNN、 強化学習の DQN でも Huber 損失が事実上の標準。

Q7：箱ひげ図の「外れ値（whiskers の外の点）」は除外すべき？

回答：箱ひげ図の外れ値判定（Q1 − 1.5×IQR より下、 Q3 + 1.5×IQR より上）は 視覚補助のための慣習的閾値であり、 統計的判定ではない。 正規分布データでも 0.7% は外れ値判定される。 「視覚的に検出された外れ値の個数や位置に意味があれば調査」「集計上は除外せずロバスト推定で対応」が標準作法。

Q8：A/B テストでロバスト統計を使うときの注意点

回答：(1) サンプルサイズ設計は中央値・トリム平均の SE で計算（平均より大きめに必要）、 (2) Mann-Whitney U 検定または permutation test を使う、 (3) winsorized 平均なら winsorize の閾値を実験開始前に固定（事後決定は p-hacking）、 (4) 効果量も中央値の差で報告。

📚 関連グループ教材

• 論文一覧トップへ戻る — 都道府県データで「東京の影響をどう扱うか」が問題になる再現論文集
• 平均ページ： 普通の平均が東京問題でなぜ歪むかを実感
• 中央値ページ： 最強のロバスト推定量を深掘り
• 外れ値ページ： 検出・除外・重み付けの判断基準
• ブートストラップページ： ロバスト推定量の SE を求める実践
• 箱ひげ図ページ： ロバスト統計の図的表現

🗺 概念マップ — ロバスト統計の家系図

ロバスト統計（Robust Statistics, Huber 1964）
├── 位置パラメータの頑健化
│   ├── 順位ベース → 中央値、 トリム平均、 ウィンザライズ平均
│   └── M-estimator
│       ├── Huber (区分線形)
│       ├── Tukey biweight (再降下)
│       ├── Hampel (3 段階)
│       └── L1 = LAD = median
├── スケール（散らばり）の頑健化
│   ├── MAD（中央絶対偏差、 BP=50%）
│   ├── IQR（BP=25%）
│   └── Q_n / S_n 推定量（Rousseeuw & Croux）
├── 回帰の頑健化
│   ├── LAD（L1）回帰
│   ├── Huber 回帰、 RLM
│   ├── LMS、 LTS（高 BP 回帰）
│   ├── MM-estimator
│   └── RANSAC（CV 由来）
├── 多変量手法の頑健化
│   ├── MCD（Min Covariance Determinant）
│   ├── ロバスト PCA
│   └── ロバスト Mahalanobis 距離
└── 検定の頑健化
    ├── Wilcoxon、 Mann-Whitney
    ├── Kruskal-Wallis
    └── permutation tests
  

系譜の祖先：Huber (1964) "Robust Estimation of a Location Parameter" でロバスト統計の数学的基盤が築かれ、 その後 Hampel、 Rousseeuw、 Tukey らによって理論と実用が拡張されました。 現代では機械学習にも組み込まれ（Huber 損失、 RANSAC）、 さらに「外れ値検出」「敵対的学習に対する頑健性」の研究へと発展しています。

📚 さらに学ぶには

このサイト内

• 論文一覧に戻る — ロバスト統計を実際に使った再現論文をハンズオン形式で読む
• 関連用語ページ — このページの「🔗 関連用語」から派生

推奨書籍・教材

• Huber, P. J. & Ronchetti, E. M. (2009) Robust Statistics (2nd ed.), Wiley. ── ロバスト統計の聖典。 数学的に厳密
• Hampel, Ronchetti, Rousseeuw & Stahel (1986) Robust Statistics: The Approach Based on Influence Functions, Wiley
• Maronna, Martin, Yohai & Salibián-Barrera (2019) Robust Statistics: Theory and Methods (with R) (2nd ed.). R 実装込みの実用書
• 『統計的データ解析の基礎』（小西貞則、 共立）── 日本語でロバスト推定の基礎を学べる

オンライン教材

• statsmodels.robust 公式ドキュメント── RLM、 Huber、 Tukey の Python 実装
• scikit-learn HuberRegressor / RANSACRegressor ── 機械学習側からの実装
• R パッケージ `robustbase` ── LMS、 LTS、 MCD などの定番
• 『現代統計学ハンドブック』（朝倉書店）── ロバスト統計章で日本語の体系整理

困ったときは

「Huber と Tukey のどちらが良いか」「k や c のチューニング値はどうするか」「ロバスト回帰の SE をどう報告するか」など、 ロバスト統計は判断項目が多い分野です。 (1) まずデータ可視化で外れ値の構造を把握、 (2) いくつかの推定量を並列に計算して比較、 (3) ブートストラップで SE を補強、 (4) 最終的に「平均と中央値で結論が違わないか」を確認、 という流れが標準的。

🧠 破壊点 — ロバスト性の数学的定義

ロバスト統計の核心は 「破壊点 (breakdown point)」 という概念です。 破壊点とは、 推定量の値を任意に大きくするのに必要な汚染データの最小割合。 標本平均の破壊点は $1/n$（1 個の外れ値で台無し）、 中央値は $50\%$ で最強。 SSDSE-B-2026 で「東京」を含めるかどうかで平均人口が約 2 倍変わる現象は、 標本平均の破壊点が極めて低いことを示しています。

$$\epsilon^*(T_n) = \min\left\{\frac{m}{n} : \sup_{x_1, \ldots, x_m} |T_n(x_1, \ldots, x_n)| = \infty\right\}$$

Hampel (1968) はこの概念を導入し、 「破壊点が高い推定量」を体系的に探す枠組みを作りました。 中央値、 トリム平均、 M 推定量、 S 推定量、 MM 推定量へと発展していきます。

📊 ロバスト推定量の総合比較表

「破壊点が高く、 かつ正規時の効率が高い」を両立するのが MM 推定量（破壊点 50%、 効率 95%）。 現代ロバスト統計の到達点です。

📐 影響関数 — 感度の数学

影響関数 (Influence Function) は「1 サンプル増やしたら推定量がどう動くか」の感度を表す関数：

$$\mathrm{IF}(x; T, F) = \lim_{\epsilon \to 0} \frac{T((1-\epsilon)F + \epsilon \delta_x) - T(F)}{\epsilon}$$

標本平均の影響関数は $\mathrm{IF}(x; \bar X, F) = x - \mu$、 つまり「外れ値ほど大きく影響」。 中央値の影響関数は有界（$\pm c$ の階段）、 これがロバスト性の源です。

Huber 関数のロバスト推定では、 影響関数が「中央で線形、 外側で一定」（クリップされた線形）。 これにより外れ値の影響を限定しつつ、 中心では効率を保ちます。

🐍 5 種類の頑健推定量を実装

SSDSE-B-2026 の県内総生産で複数の頑健推定量を比較：

import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.api as sm

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=1)
gdp = df.iloc[:, 11].values

print(f'平均: {gdp.mean():,.0f}')
print(f'中央値: {np.median(gdp):,.0f}')
print(f'10% トリム平均: {stats.trim_mean(gdp, 0.1):,.0f}')
print(f'Huber: {sm.robust.scale.huber(gdp)[0]:,.0f}')

# ロバスト回帰
X = sm.add_constant(df.iloc[:, 3].values)  # 人口
rlm = sm.RLM(gdp, X, M=sm.robust.norms.HuberT()).fit()
print(rlm.params)

東京を含む 47 都道府県データで平均値は中央値より大きく出ます。 トリム平均と Huber 推定量は両者の中間で、 「東京の影響を抑えつつ全データを使う」バランスを取ります。

📈 ロバスト回帰 — 7 手法

回帰での頑健化：

• M 推定 (Huber)：残差の二乗ではなく Huber 損失を最小化。
• LMS (Least Median of Squares)：残差二乗の中央値を最小化。 破壊点 50%。
• LTS (Least Trimmed Squares)：上位 $h$ 個の残差二乗を最小化。
• S 推定量：残差のロバストスケール推定量を最小化。
• MM 推定量：S 推定量 → M 推定量の二段階。 高破壊点 + 高効率。
• RANSAC：ランダム標本でモデル試行、 インライアー多数のものを採用。 コンピュータビジョンで標準。
• Theil-Sen 推定量：全ペアの傾きの中央値。 単回帰で破壊点 29.3%。

🛡 多変量ロバスト共分散行列

多変量データでのロバスト共分散行列推定：

• MCD (Minimum Covariance Determinant)：行列式最小のサブセットから共分散を計算。 Rousseeuw (1984)。
• MVE (Minimum Volume Ellipsoid)：最小体積楕円体。
• S 推定量：M 推定量の多変量拡張。
• OGK (Orthogonalized Gnanadesikan-Kettenring)：高速ロバスト共分散。

from sklearn.covariance import MinCovDet

X = df.iloc[:, 3:11].values
mcd = MinCovDet(random_state=0).fit(X)
print(f'ロバスト平均ベクトル: {mcd.location_}')
# Mahalanobis 距離で外れ値検出
d = mcd.mahalanobis(X)
outliers = np.where(d > stats.chi2.ppf(0.975, X.shape[1]))[0]
print(f'外れ値県: {df.iloc[outliers, 1].tolist()}')

SSDSE-B-2026 を MCD で分析すると、 東京・大阪・愛知が頑健 Mahalanobis 距離で外れ値として識別されます。

⚠️ ロバスト統計の追加的落とし穴（8 個）

• 「ロバスト = 外れ値を捨てる」と誤解：ロバスト推定は「重みを下げる」だけで、 完全除外ではない。
• 効率の損失を過小評価：MAD は正規時 SD より 37% の効率しかない。 大標本でカバー可能だが、 小標本では深刻。
• 初期値依存：M 推定の反復は初期値次第で別解に。 中央値スタートが定石。
• 破壊点の誤理解：「破壊点 50%」は「データの 50% が外れ値でも壊れない」意味であり、 「常に正しい答え」ではない。
• 外れ値の構造を見ない：ロバスト推定は外れ値を弾くが、 なぜ外れ値かは教えてくれない。 必ず可視化と原因究明。
• 非対称分布での適用：ロバスト推定は対称分布前提。 歪んだ分布では別途処理が必要。
• 多変量での落とし穴：1 次元で頑健でも多変量で「外れ値の塊」（masking）に騙される。 MCD など多変量ロバスト法を使う。
• 計算コスト過小評価：MM 推定量や MCD は OLS より 10-100 倍重い。

🏛 ロバスト統計の歴史 — Tukey から adversarial training まで

Tukey は 1960 年に「Conventional statistics have not been particularly successful in dealing with the data we have」と書き、 ロバスト統計の必要性を訴えました。 Huber (1964) が M 推定量を提案、 Hampel (1968) が破壊点を導入、 Rousseeuw & Leroy (1987) が LMS/LTS を体系化、 Yohai (1987) が MM 推定量で破壊点と効率を両立――半世紀でロバスト統計は完成した分野になりました。

• 1960：Tukey がロバスト統計の必要性を提唱。
• 1964：Huber が M 推定量を提案。
• 1968：Hampel が破壊点・影響関数の概念を導入。
• 1981：Rousseeuw が LMS を提案、 破壊点 50% を達成。
• 1984：Rousseeuw が MCD を提案、 多変量へ拡張。
• 1987：Yohai が MM 推定量で「破壊点 50% + 効率 95%」を実現。
• 2000 年代：機械学習で robust loss が再注目（Huber loss、 Tukey biweight loss）。
• 2020 年代：深層学習で adversarial training と関連、 「分布外ロバスト性」が新潮流。

🎯 影響関数と RANSAC の詳細

🔗 関連用語（拡張）

• 中央値：最もシンプルなロバスト推定量。
• MAD：ロバスト分散の標準。
• 外れ値：ロバスト統計の主敵。
• 外れ値処理：ロバスト法と相補的。
• 分位点：ロバスト統計量の基本要素。
• ブートストラップ：ロバスト推定の SE 計算に。
• 回帰分析：ロバスト回帰の応用先。
• Huber 損失：M 推定の代表的損失関数。

❓ よくある質問（FAQ）

🏛 ロバスト統計 の歴史

• 1822：Gauss: 誤差論で「最小二乗法はガウス分布前提」と明言、 既にロバスト性問題を認識。
• 1960：Tukey: A Survey of Sampling from Contaminated Distributions で問題提起。
• 1964：Huber: Robust Estimation of a Location Parameter で M 推定量を提案。
• 1968：Hampel: Contributions to the Theory of Robust Estimation 博論で破壊点・影響関数を導入。
• 1981：Huber: Robust Statistics 教科書、 分野の確立。
• 1984：Rousseeuw: LMS 推定量で破壊点 50% を達成。
• 1984：Rousseeuw & Yohai: S 推定量で破壊点と効率の同時改善。
• 1987：Yohai: MM 推定量で「破壊点 50% + 効率 95%」を実現。
• 1990：Hampel et al.: Robust Statistics: The Approach Based on Influence Functions 教科書。
• 2006：Maronna et al.: Robust Statistics: Theory and Methods 現代の決定版教科書。
• 2019：Hampel-Yohai 後継: 機械学習でのロバスト損失関数、 adversarial training と接続。

💼 産業応用事例

🛡 影響関数の図解

影響関数 (IF) は「1 つのデータ点を $x$ に置いたとき、 推定量がどれだけ変化するか」を表します。 各推定量の IF：

「再下降型」（redescending）は最も強力なロバスト性。 Tukey biweight は外れ値の影響を完全に 0 にできる一方、 多重最適解の問題を持ちます。 これが MM 推定量で「S 推定量で初期値 → biweight で最適化」する理由です。

ロバスト統計のフロー

データに外れ値ありか？
├── ない（正規） → 通常の SD / 平均で OK
├── 軽い汚染（< 10%） → トリム平均、 Huber M
├── 中程度（10-25%） → MAD、 M 推定量
└── 重度（> 25%）
    ├── 単変量 → 中央値、 LMS
    ├── 回帰 → MM 推定量、 LTS
    └── 多変量 → MCD、 OGK