🔖 キーワード索引

本ページで扱う中核キーワードを 「概念 / 推定量 / 評価指標 / 派生手法 / 実装ライブラリ」 の 5 系統で索引化した。 チップをクリックすると該当セクションへ移動する。 別名（Aliases）は丸括弧で併記。

別名 (Aliases): ロバスト推定、 抗外れ値統計、 R 推定、 高ブレイクダウン点推定、 などの呼び方が文献では混在する。 本ページでは 「ロバスト統計 = 外れ値や分布の歪みに対して結論が大きく変わらない推定・検定の総称」 として扱う。

🧭 現場での意思決定フロー

「ロバスト統計を使うべきか」は 4 つの判断軸で決まる。 単に外れ値が "見える" だけで M 推定に飛びつくと、 効率（標準誤差）を 1〜2 割犠牲にする可能性がある。 まず分布と業務要件を整理した上で、 以下のフローで段階的に判断する。

判断軸 1: データ生成過程（DGP）が正規からどの程度ずれているか

SSDSE-B-2026 の都道府県人口は、 東京・大阪・愛知が右側に長い裾を持つ重尾分布である。 シャピロ–ウィルク検定で p<0.001、 歪度=2.22、 尖度=7.95。 こうしたケースでは平均は東京 1 県に引きずられ「全国平均=265 万人 (＞ 中央値 155 万人の 1.71 倍)」となる。 ここで平均を使うと「平均より下=下位 35 県」という直感に反する読みになり、 報告の説得力を失う。 中央値 155 万人 + 第 1・3 四分位 103 万 / 264 万を併記すれば、 読者は「中位の県の規模感」と「ばらつき」を同時に得られる。

判断軸 2: 外れ値が「測定誤差」か「重要な少数派」か

外れ値=ゴミという思い込みは危険である。 たとえば医療費の右裾は 高額長期入院患者という臨床上重要なグループ を含む。 ここでロバスト推定で押し込めると、 政策議論で必要な「上位 5% の負担分布」が消える。 ロバスト統計は 「結論を平均値で要約したいが、 外れ値の存在は別途報告する」 という二段構えで使うのが正しい。

判断軸 3: 推定の用途（記述統計 / 推論 / 予測）

記述（県の典型を示す）→ 中央値 + IQR で十分。 推論（処置効果の検定）→ Wilcoxon 順位和検定など分布フリー手法、 もしくは Huber-White SE。 予測（回帰モデル）→ HuberRegressor / RANSAC / Theil–Sen。 用途に合わせて 3 系統を切り替える。

判断軸 4: サンプルサイズ n

n=47 県のような小標本では、 ロバスト推定の漸近理論（n→∞）が頼りにくい。 ブートストラップ 1,000 回で M 推定量の信頼区間を構築するのが現実的。 大標本（n>1,000）なら理論 SE をそのまま使ってよい。

この 4 軸表を最初に当てはめれば、 「とりあえず平均」も「とりあえず中央値」も避けられる。 結論を 1 つの統計量で要約するのではなく、 分布の形・外れ値・用途・n の 4 軸を全て言語化してから推定量を選ぶ 習慣がロバスト統計の本質である。

📈 感度分析プロトコル

ロバスト統計を「使う」のではなく、 結果が 外れ値に対してどの程度安定か を体系的に検証するのが感度分析（Sensitivity Analysis）である。 SSDSE-B-2026 の県データで 5 段階のプロトコルを示す。

ステップ 1: 全データで通常の推定量を計算

県人口 47 件で平均=265 万、 中央値=155 万、 標準偏差=280 万、 MAD=62 万。 平均と中央値の差 110 万、 SD と MAD の比 4.5 が 「分布が歪んでいる証拠」 となる。

ステップ 2: 上位 / 下位 5% を順次除外し再推定

トリム率 α を 0%, 5%, 10%, 20% と段階的に変えて平均（トリム平均）を計算。 α=0: 265 万 → α=5: 232 万 → α=10: 211 万 → α=20: 169 万。 単調減少が緩むのが「真の中心」の目安である。

ステップ 3: 外れ値除外前後で回帰係数の比較

総人口 X に対する延べ宿泊者数 Y（G7101）の OLS 回帰係数 β=4.11（宿泊/人）。 東京を除外すると β=2.94、 上位 3 県除外で β=2.61。 ここまで動くなら HuberRegressor で全データを使うのが正解。 β=3.97 (Huber) に落ち着く。

ステップ 4: ブートストラップで信頼区間

中央値の 95% CI を 10,000 回のブートストラップで構築。 [148 万, 178 万] が得られたら、 「県の中位的人口は 150〜180 万の範囲」と政策議論で使える表現が確定する。

ステップ 5: 結果のレポート

必ず 「平均 / 中央値 / トリム平均 / Huber 推定量」を 1 表に並べる。 4 つが概ね一致 → 結論は安定。 大きくズレる → 「分布の歪みと外れ値が結論に影響する」と明記する。 後者の場合、 ロバスト統計の真価が発揮される。

この感度分析を 毎回データを触る前に習慣化する ことが、 報告の信頼性を一段引き上げる。 ロバスト統計とは個別の関数や推定量を指す名前である以前に、 「平均・中央値・トリム平均・Huber 推定量を並べ、 結論が外れ値の有無で変わらないことを毎回確認する文化」 そのものなのである。 SSDSE-B-2026 で東京を含む 47 件を扱う限り、 この感度分析は必ず通る検問所となる。

🧪 SSDSE-B-2026 で見るロバスト統計の威力

前節までで理論と判断基準を整理した。 ここでは SSDSE-B-2026（47 都道府県 × 90 指標） の具体データで「平均と中央値の差」「ロバスト推定の効果」を 5 つの指標について検証する。 各ケースで「なぜロバスト統計が必要か」が腑に落ちる構成にした。

ケース 1: 総人口（A1101）— 極端な右裾分布

東京 1,409 万人、 神奈川 923 万人、 大阪 876 万人。 最下位は鳥取県 54 万人。 平均 265 万、 中央値 155 万、 比率 1.71。 SD=280 万、 MAD=62 万、 比率 4.5。 ヒストグラムは典型的な対数正規型に近い。 ここで「県の代表的な規模は？」と問われたら、 中央値 155 万人を提示するのが正しい。 平均 265 万人は「東京込みの算術平均」であり、 県の典型ではない。

ケース 2: 高齢化率（A1303）— 概ね対称、 ロバスト不要

高齢化率（A1303 高齢者人口 ÷ A1101 総人口）は 22.8%〜39.1% の範囲に 47 県が分布し、 平均 31.6%、 中央値 31.8% で差は 0.2 ポイント。 SD=3.3、 MAD=2.1、 比率 1.6。 総人口のような極端な右裾はなく、 ほぼ対称。 ここではロバスト推定の出番はほぼなく、 平均 + 標準偏差で十分。 「ロバスト統計を使うべきかは指標ごとに判定する」 という原則を象徴するケースである。

ケース 3: 延べ宿泊者数（G7101）— 重尾だが対数変換で改善

東京 8,027 万人泊、 大阪 4,401 万人泊。 最下位は徳島 170 万人泊。 平均 1,064 万、 中央値 585 万、 比率 1.82。 対数変換すると歪度が 3.32 → 0.70 に下がり、 ほぼ対称化。 ロバスト推定と対数変換は 代替手段 として位置づけられる。 「対数で対称になるなら通常の OLS でよい」 vs 「対数変換が解釈を難しくするならロバスト推定」のトレードオフで判断する。

ケース 4: 転入者数（A5101）— 外れ値は政策上重要

東京 40.7 万人、 神奈川 21.1 万人、 埼玉 16.1 万人、 大阪 16.0 万人、 千葉 14.0 万人が上位。 最下位は鳥取 0.8 万人。 平均 4.7 万、 中央値 2.4 万で右裾が長い。 問題は「上位 5 県だけで全国転入の 48%」という一極集中の政策論点。 ここでロバスト推定で「典型県」を出すと 議論の重要な部分（東京圏への人口移動）が消える。 中央値ではなく 分位点（q90=11.9 万, q99=31.7 万）+ 中央値の併記 が正解である。

ケース 5: 延べ宿泊者数（G7101）— 外れ値耐性回帰

総人口 X に対する延べ宿泊者数 Y を回帰。 OLS では β=4.11、 R²=0.71。 東京を除外すると β=2.94、 R²=0.53。 上位 3 県（東京・神奈川・大阪）を除外すると β=2.61、 R²=0.41。 ここで HuberRegressor (デフォルト ε=1.35) を全データに適用すると β=3.97。 「外れ値の影響を中和しつつ全データを使う」 ことができる。 RANSACRegressor では β=2.26（観光集中の 14 県を外れ値とみなして除外）。 用途次第で Huber か RANSAC を選ぶ。

この 5 ケースで言えるのは、 ロバスト統計は 「全データを保持しつつ外れ値の影響を中和する技術」 であり、 単純除外より優れること。 ただし指標ごとに必要性が異なり、 一律適用は推奨しない。 SSDSE-B のような公的データで分析する際、 最低でもケース 1〜5 のような 「指標ごとの分布診断 + 推定量の選定」 を行う習慣をつけたい。

📜 ロバスト統計の歴史と理論的背景

ロバスト統計の起源は古く、 1818 年に Laplace が中央値の漸近分散を導出 したのが学術的な原点とされる。 19 世紀後半には Edgeworth が「外れ値の影響を最小化する推定」を議論し、 20 世紀初頭には Fisher が最尤推定の枠組みで議論を進めた。 しかし「ロバスト統計」という名称と統一的理論を提唱したのは、 1964 年の Peter J. Huber の論文 "Robust Estimation of a Location Parameter" である。

Huber (1964) のブレークスルー

Huber は「真の分布が正規分布から ε の混合で歪んでいるとき、 推定量の最大漸近分散を最小化するミニマックス問題」を解いた。 結果として導かれたのが Huber loss であり、 これは中心では二乗誤差、 裾では絶対値誤差というハイブリッド構造を持つ。 1964 年から 60 年が経過した今も、 機械学習の損失関数（Smooth L1, Pseudo-Huber, Faster R-CNN の bounding box regression）で標準採用されている。 理論的優美さと実用性の両立がロバスト統計の核心である。

Hampel (1971) — 影響関数の導入

Frank Hampel は 「影響関数 (Influence Function, IF)」 を導入し、 推定量の頑健性を関数として可視化する枠組みを提供した。 IF(x; F, T) は「分布 F に微小な汚染 x を加えたとき、 推定量 T がどれだけ動くか」を表す。 平均の IF は線形（無制限に動く）、 中央値の IF は階段関数（±1/(2f(median))）、 Huber の IF はクリップ関数（c で打ち切り）。 これによりロバスト推定量の 「頑健性の度合いを定量化」 できるようになった。

Rousseeuw (1984) — 高ブレイクダウン点推定

Peter Rousseeuw は LMS (Least Median of Squares) 回帰 と MCD (Minimum Covariance Determinant) を提案。 これらは 50% のブレイクダウン点 を達成する最初の実用的アルゴリズムだった。 LMS は「残差平方の中央値を最小化する回帰」、 MCD は「h サブセットの共分散行列式が最小になる h 点を選ぶ」アルゴリズム。 計算量は当初 O(n²) で重かったが、 高速 MCD（FAST-MCD, 1999）で実用化された。

Yohai (1987) — MM 推定量

Victor Yohai は 「高ブレイクダウン点（50%）と高効率（95%）を両立する MM 推定量」 を提案した。 これは 2 段階で構成される: (1) LMS / S 推定でロバストな初期値とスケール推定、 (2) M 推定で効率の高い最終推定。 統計ソフト R の robustbase::lmrob() や Python の statsmodels.robust.linear_model.RLM で標準実装されている。

現代の発展 — 機械学習との融合

2010 年代以降、 ロバスト統計は深層学習の文脈で再注目された。 物体検出 (Fast R-CNN, 2015) は Smooth L1 loss を採用、 GAN の Wasserstein 距離 (2017) はロバスト推定の派生、 異常検知でも MCD と Isolation Forest が組み合わせて使われる。 さらに 差分プライバシー (Differential Privacy) や 連合学習 (Federated Learning) でも、 「外れ値≒攻撃データ」を中和するためにロバスト統計が使われている。 古典理論が現代的問題で生き続けている好例である。

この歴史的経緯を踏まえると、 ロバスト統計は 「単なる古典手法」ではなく、 現代の AI 時代で再評価されている深い理論体系」 と理解できる。 入門段階で中央値・MAD だけ覚えて終わるのではなく、 影響関数・BP・MM 推定量・そして深層学習との接続まで視野に入れて学ぶことを強く推奨する。

📖 推奨教科書と論文 5 選

独学・大学院演習用に最適な教材を 5 件紹介する。 (1) Huber & Ronchetti (2009) "Robust Statistics, 2nd ed." Wiley — 創始者本人による決定版。 (2) Maronna, Martin, Yohai, Salibián-Barrera (2019) "Robust Statistics: Theory and Methods, 2nd ed." Wiley — 最新の MM 推定量と R 実装。 (3) Hampel et al. (1986) "Robust Statistics: The Approach Based on Influence Functions" Wiley — 影響関数の古典的教科書。 (4) Wilcox (2017) "Introduction to Robust Estimation and Hypothesis Testing, 4th ed." Academic Press — 学部生向けの分かりやすい入門書。 (5) Rousseeuw & Leroy (1987) "Robust Regression and Outlier Detection" Wiley — LMS / MCD の原典。 これら 5 冊を順に読めば理論的にも実装的にも完全に習得できる。

🧰 Python / R のロバスト統計ライブラリ早見表

実務で使うときの定番ライブラリ。 Python: scipy.stats.median_abs_deviation, scipy.stats.trim_mean, sklearn.linear_model.HuberRegressor / RANSACRegressor / TheilSenRegressor, sklearn.covariance.MinCovDet, statsmodels.robust.scale.mad, statsmodels.robust.linear_model.RLM。 R: robustbase::lmrob(), MASS::rlm(), WRS2::trimmed.mean(), robust::lmRob()。 SSDSE-B-2026 を使った演習では Python の sklearn.linear_model.HuberRegressor から始めるのが学習効率が最も高い。 fit メソッド 1 行で OLS のロバスト版が手に入る。

❓ よくある質問

Q1: 平均値を全てロバスト統計に置き換えた方がよいか？
A1: いいえ。 正規分布に近いデータでは平均値は最も効率的（分散最小）であり、 中央値に置き換えると標準誤差が 1.25 倍に増える（漸近相対効率 ARE=0.637）。 分布の形と用途に応じて使い分ける のが正解で、 万能手法は存在しない。

Q2: Huber と Tukey biweight、 どちらを使うべきか？
A2: 経験則として、 外れ値が 連続的な裾（重尾分布）→ Huber、 外れ値が 離散的な異常値（測定エラー）→ Tukey biweight。 後者は閾値外を完全に切り捨てるため極端な外れ値に強いが、 微分可能でも凸ではないので局所最適に陥る可能性がある（初期値依存）。

Q3: ロバスト統計で検定の有意水準は変わるか？
A3: 変わる。 通常の t 検定の p 値はサンプル分布が正規であることを前提にする。 重尾分布で t 検定を使うと第一種過誤が 真の α より小さく なる（保守的）。 ノンパラメトリック検定（Wilcoxon、 順位和）かブートストラップ p 値を使うのが安全。

Q4: 機械学習モデルでも同じ考え方が通用するか？
A4: 通用する。 sklearn の HuberRegressor、 RANSACRegressor、 TheilSenRegressor はそれぞれ M 推定・コンセンサスサンプリング・順序ベースの 3 系統。 XGBoost / LightGBM では objective='reg:pseudohubererror' を選べば同様の効果が得られる。 深層学習でも SmoothL1Loss （Huber の派生）が物体検出 (R-CNN 系) で標準採用される。

Q5: 学習にどのくらいかかるか？
A5: 中央値・MAD・四分位数の意味を直感的に理解 → 1 日。 M 推定量と Huber loss を一つの例題で実装 → 2 日。 MM 推定量・MCD を sklearn / statsmodels で動かす → 4 日。 影響関数の図を自分で描く → 1 日。 計 1-2 週間で「ロバスト統計の文化」が身につく。 教科書は Maronna et al. (2006) Robust Statistics: Theory and Methods、 入門なら Wilcox (2017) Introduction to Robust Estimation。

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

日本の都道府県データを扱うと、 ほぼ必ず 「東京問題」に直面します。 総人口でも、 出生数でも、 消費支出でも、 東京は他 46 県と桁が違うレベルで大きい。 SSDSE-B-2026 の A1101（総人口）を例に、 47 県の 普通の平均 を計算すると：

SSDSE-B-2026 総人口（A1101）の代表値

「普通の県の総人口はだいたいどれくらい？」と問われたとき、 265 万人 と答えるのと 155 万人 と答えるのとでは、 意味がまったく違います。 平均は「東京を含めた全体の総量を 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% を超える推定量は原理的に作れません（半分以上が「壊れた値」なら、 もはやそれが多数派で真の値）。 中央値は最強です。

🎨 概念図で押さえるロバスト統計の核（補遺）

ロバスト統計の中核概念（外れ値の影響、 損失関数の形、 ブレイクダウン点）を 3 枚の図で再整理する補遺セクション。

🖼 図 1: 平均値 vs 中央値 — 外れ値 1 つで何が起きるか

古典統計の代表値「平均値」は 外れ値 1 つで任意に動かせる。 これを定量化したのがブレイクダウン点 0%。 一方の中央値はブレイクダウン点 50% で、 半数まで汚染されても揺るがない。

💬 読み方：上段 5 点では平均＝中央＝190 で一致。 下段で外れ値（赤）1 点を足すと平均は 268 へジャンプ（+78）するが、 中央値は 190 のまま。 これがブレイクダウン点 0% vs 50% の差である。

🖼 図 2: 二乗損失 vs Huber 損失 vs 絶対値損失

ロバスト推定の核は 損失関数の形を変えることにある。 二乗損失（L2）は外れ値に二次的に反応するため敏感、 絶対値損失（L1）は線形反応、 Huber 損失は中央が L2・外側で L1 に切り替わる「いいとこ取り」型。

💬 読み方：赤の L2（二乗）は外側で 放物線的に 急上昇 → 外れ値 1 点が全体推定を激しく揺らす。 緑 L1（絶対値）は線形なので外れ値の影響は有限。 青破線 Huber は |r|≤c で L2、 |r|>c で L1 に切り替わる連続関数で、 中心付近の効率と外れ値耐性を両立する。

🖼 図 3: ブレイクダウン点で見る推定量ランキング

ブレイクダウン点（BP）は 「データの何 % まで汚染されても推定が暴走しないか」を示す指標。 0% が「1 点で破綻」、 50% が「半数まで耐える」。 50% が理論上限である。

💬 読み方：赤（BP=0%）の平均・最小二乗は 1 点で破綻、 オレンジ（BP=10%）のトリム平均は α=20% カットで耐える、 緑（BP=50%）の中央値・MAD・LMS 回帰は 理論上限 に達する。 ロバスト性と効率のトレードオフを意識して使い分ける。

以上、 ロバスト統計の本質（外れ値感受性・損失関数・BP）を 3 枚の図で整理した。 詳細は本文参照。

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

① 中央値 と MAD

② トリム平均（trimmed mean）

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

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

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

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

⑤ Tukey の biweight（再降下型）

⑥ ロバスト回帰：LMS と LTS

📐 影響関数 — 感度の数学

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

📐 数学的詳細 — M 推定量の収束と漸近理論

M 推定量 $\hat{\theta}_n$ は $\sum_i \psi(x_i - \theta) = 0$ の解として定義されます。 ここで $\psi$ は影響関数の心臓部となる関数。 例えば Huber では:

$$\psi_c(u) = \begin{cases} u & |u| \le c \\ c \cdot \mathrm{sign}(u) & |u| > c \end{cases}$$

$c$ はクリッピング閾値。 $c \to \infty$ で OLS の正規方程式に一致、 $c \to 0$ で中央値方程式に近づきます。 漸近的には:

$$\sqrt{n}(\hat{\theta}_n - \theta_0) \xrightarrow{d} \mathcal{N}\!\left(0, \, \frac{\mathbb{E}[\psi^2(X-\theta_0)]}{\{\mathbb{E}[\psi'(X-\theta_0)]\}^2}\right)$$

分子は「変動」、 分母は「感度」を表します。 この比 $V(\psi, F)$ を最小化する $\psi$ を求めると Huber 関数が導かれる――これが Huber (1964) の minimax 定理の中身です。 一般化するとノンパラメトリック効率限界に到達でき、 ロバスト統計は単なる対症療法ではなく、 漸近最適理論として完成しています。

🔬 重み関数の比較

🌀 IRLS（Iteratively Reweighted Least Squares）

M 推定量を解く標準アルゴリズムが IRLS です:

1. 初期値 $\hat{\theta}^{(0)}$（例：中央値）
2. 残差 $r_i^{(k)} = x_i - \hat{\theta}^{(k)}$ を計算
3. 重み $w_i = \psi(r_i)/r_i$ を計算
4. 重み付き最小二乗 $\hat{\theta}^{(k+1)} = \sum_i w_i x_i / \sum_i w_i$
5. 収束まで繰り返し

凸 $\psi$（Huber 含む）なら大域収束保証あり。 再下降型（Tukey biweight）は局所最適に陥るので初期値依存性が強く、 通常は Huber を初期値として使う。 statsmodels の RLM(...).fit() は内部で IRLS を実装。

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

🔬 数式を言葉で読み解く（ロバスト統計を 4 ステップで比較）

ロバスト統計の核は「外れ値が混入しても結論が大きく崩れない」性質。 中央値・MAD・Huber 推定量を SSDSE-B-2026 で 4 ステップで比較する。

▶ ブロック 1: 平均 vs 中央値 (代表値の感度比較)

このコードでやること: SSDSE-B-2026 の県別総人口について平均と中央値を比較。 東京・大阪のような大都市が平均をどれだけ押し上げるかを見る。

📥 入力データ (抜粋):

1
2
3
4
5
6
7
8
9

import pandas as pd
import numpy as np

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', skiprows=[1], encoding='cp932')
pop = df[df['SSDSE-B-2026']==2023]['A1101'].values / 10000  # 総人口（万人）

print(f'平均: {pop.mean():.1f} (万人)')
print(f'中央値: {np.median(pop):.1f} (万人)')
print(f'差分: {pop.mean() - np.median(pop):.1f} (万人、 平均が押し上げられている)')

📤 実行例:

💬 平均 265 万 vs 中央値 155 万。 差 110 万は東京・大阪・神奈川などの大都市が平均を押し上げている証拠。 「典型的な県」を語るなら中央値の方が適切。

▶ ブロック 2: 標準偏差 vs MAD (ばらつきの感度比較)

このコードでやること: 標準偏差と MAD (中央絶対偏差) を比較。 標準偏差は二乗で平均するため外れ値で膨らみ、 MAD は中央値ベースで頑健。

1
2
3
4
5
6
7
8

from scipy.stats import median_abs_deviation

std = pop.std(ddof=1)
mad = median_abs_deviation(pop)
print(f'標準偏差: {std:.0f}')
print(f'MAD (中央絶対偏差): {mad:.0f}')
print(f'MAD × 1.4826 (正規分布で σ 相当): {mad*1.4826:.0f}')
print(f'std / (1.4826·MAD) = {std / (1.4826*mad):.2f}  (1 から離れる ≈ 外れ値あり)')

📤 実行例:

💬 std/MAD 比 3.03 → 標準偏差が正規化 MAD の 3 倍、 強い外れ値の影響。 ロバスト統計の世界では、 ばらつき指標として MAD を使うのが標準。

▶ ブロック 3: HuberRegressor で外れ値に強い回帰

このコードでやること: 総人口から延べ宿泊者数（G7101）を予測。 通常の OLS と Huber 回帰を比較し、 東京・観光集中県という外れ値が係数に与える影響を観察。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

from sklearn.linear_model import LinearRegression, HuberRegressor

dd = df[df['SSDSE-B-2026']==2023].copy()
X = dd[['A1101']].astype(float).values / 1e6  # 総人口（百万人）
y = dd['G7101'].astype(float).values / 1e6  # 延べ宿泊者数（百万人泊）

ols = LinearRegression().fit(X, y)
hr = HuberRegressor(epsilon=1.35).fit(X, y)

print(f'OLS 係数: {ols.coef_[0]:.5f}, 切片: {ols.intercept_:.1f}')
print(f'Huber 係数: {hr.coef_[0]:.5f}, 切片: {hr.intercept_:.1f}')
print(f'差: 東京・観光集中県の影響で OLS の傾きが上振れ')

📤 実行例:

💬 OLS と Huber で傾きが変わる。 東京 (延べ宿泊者数 8,027 万人泊) のような回帰直線から大きく上に外れる点が OLS の傾きを引っ張っている。 Huber は残差の大きい点の重みを下げて影響を抑えるため、 傾き 4.11 → 3.97 と地方部の「人口対宿泊者数」の関係に近づく。

▶ ブロック 4: 外れ値を人為的に追加して頑健性を比較

このコードでやること: 既存データに人為的な異常値 (例: 観測ミスで人口が 10 倍) を 1 件追加。 平均は爆発するが中央値は微動だにしない、 という実物実験。

1
2
3
4
5
6
7

pop_with_outlier = np.append(pop, [pop.max()*10])

print(f'元データの平均:     {pop.mean():.0f}')
print(f'外れ値追加後の平均: {pop_with_outlier.mean():.0f} (変化大)')
print()
print(f'元データの中央値:     {np.median(pop):.0f}')
print(f'外れ値追加後の中央値: {np.median(pop_with_outlier):.0f} (ほぼ不変)')

📤 実行例:

💬 平均が 2 倍以上に跳ねる一方、 中央値はほぼ動かない。 これがロバスト統計の威力。 中央値・外れ値検知・ブートストラップ と組み合わせて運用するのが定石。

🧮 実値で計算してみる — SSDSE-B 総人口

STEP 1：47 県の総人口（万人）の主要値

SSDSE-B-2026 から、 2023 年の都道府県別の総人口（A1101、 実値、 万人）を抜粋：

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

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

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

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

STEP 4：解釈の違い

• 「日本の都道府県の人口規模はどれくらい？」と問われ、 平均 265 万人と答えるのは正確だが誤解を招く。 47 県のうち 12 県程度しか 265 万人を超えない
• 中央値 155 万人なら、 「23 県が 155 万人以下、 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、 「長期累積収益」を語るなら算術平均、 と用途で使い分けるのが正解。

🧮 SSDSE-B-2026 で MM 推定量を試す（47 都道府県の総人口対延べ宿泊者数）

SSDSE-B-2026 の A1101（総人口）と G7101（延べ宿泊者数）を例に、 古典回帰 (OLS) とロバスト回帰 (MM 推定量) の差を確認します。 東京・大阪のような巨大県が回帰直線をどれだけ歪めるか、 ロバスト推定はどう違うかが実感できます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

# SSDSE-B-2026 MM推定量：総人口 → 延べ宿泊者数（東京が外れ値）
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.robust.robust_linear_model import RLM
from sklearn.linear_model import HuberRegressor, RANSACRegressor

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)  # 2023年度・47都道府県
x = df['A1101'].astype(float).values / 1e6   # 総人口（百万人）
y = df['G7101'].astype(float).values / 1e6   # 延べ宿泊者数（百万人泊）

X = sm.add_constant(x)
ols = sm.OLS(y, X).fit()
print('OLS      : slope=%.3f' % ols.params[1])

# statsmodels IRLS（Huber / Tukey biweight）
rlm_huber = RLM(y, X, M=sm.robust.norms.HuberT()).fit()
rlm_tukey = RLM(y, X, M=sm.robust.norms.TukeyBiweight()).fit()
print('Huber    : slope=%.3f' % rlm_huber.params[1])
print('Tukey    : slope=%.3f' % rlm_tukey.params[1])

# sklearn（内部で標準化）— 外れ値の影響がより顕著に出る
print('sk-Huber : slope=%.3f' % HuberRegressor(epsilon=1.35).fit(x.reshape(-1,1), y).coef_[0])
ran = RANSACRegressor(random_state=0).fit(x.reshape(-1,1), y)
print('RANSAC   : slope=%.3f (外れ %d 県)' % (ran.estimator_.coef_[0], (~ran.inlier_mask_).sum()))

# Huber で最も重みが下がった県（観光集中・大都市）
for k in np.argsort(rlm_huber.weights)[:5]:
    print(f'最低重み: {df.iloc[k]["Prefecture"]:8s} weight={rlm_huber.weights[k]:.3f}')

このコードを実行すると、 東京・埼玉・京都・神奈川・沖縄など観光集中・大都市の県が weight が 1 未満（東京 0.093）に押し下げられます。 sklearn の HuberRegressor（傾き 3.972）や RANSAC（2.257, 14 県を外れ値扱い）では OLS の 4.106 から傾きがはっきり下振れし、 古典 OLS が巨大都市・観光地に引きずられる現象 (leverage effect) を見抜けます。 一方 statsmodels の IRLS（Huber 3.982 / Tukey 3.910）は高レバレッジ点をある程度追随するため下がり方が緩やかで、 手法ごとの外れ値耐性の違いも観察できます。

🧮 数式に値を入れて手で計算する: 頑健統計量

SSDSE-B-2026 で「人口密度（百人/km²）」が極端な外れ値となる現象を再現する小例で、 5 道県（北海道, 岩手, 秋田, 高知, 鳥取）と東京都の人口密度値 [2, 4, 5, 6, 7, 100]（×100 人/km² の概数）から平均・中央・トリム平均を比較する。

Step 1: データ

Step 2: 頑健代表値

🐍 Python で再現

1
2
3
4
5
6

import numpy as np
from scipy import stats
x = np.array([2, 4, 5, 6, 7, 100])
print(f"平均: {x.mean():.2f}")
print(f"中央: {np.median(x)}")
print(f"20% trim: {stats.trim_mean(x, 0.2)}")

📤 実行結果

💬 手計算 (Step 2) 中央 5.5 / trim 5.5 と Python 出力が完全一致。

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

# 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=[1])
df = df.iloc[:, [0, 2, 3]].copy()  # 年度、都道府県、総人口（A1101 列）
df.columns = ['year', 'pref', 'pop']
df = df[df['year'] == 2023].dropna()
x = df['pop'].values / 10000  # 万人単位に換算（元データは人）

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

from statsmodels.robust.scale import huber
import numpy as np

# Huber M-estimator（位置・尺度を同時推定する Proposal 2）
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 出生数）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

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=[1])
df = df.iloc[:, [0, 2, 3, 18]].copy()
df.columns = ['year', 'pref', 'pop', 'births']
df = df[df['year'] == 2023].dropna()

x = df['pop'].values  # 総人口（人, A1101）
y = df['births'].values  # 出生数（人, A4101）
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

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

x = df['pop'].values.reshape(-1, 1)
y = df['births'].values  # 出生数（人）

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：可視化（外れ値の影響を視覚化）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

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 種類の頑健推定量を実装

SSDSE-B-2026 の総人口（A1101）で複数の頑健推定量を比較し、 総人口から出生数（A4101）へのロバスト回帰も試す：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

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])
df = df[df['SSDSE-B-2026'] == 2023]  # 2023年度・47都道府県
pop = df.iloc[:, 3].values  # 総人口（A1101）

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

# ロバスト回帰（総人口 → 出生数 A4101）
X = sm.add_constant(pop)  # 総人口
rlm = sm.RLM(df.iloc[:, 18].values, 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)：高速ロバスト共分散。

1
2
3
4
5
6
7
8
9

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 距離で外れ値として識別されます。

⚠️ ロバスト統計の 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 検定や中央値の使用を通じて医学・心理学の領域で広まり、 現在では「外れ値が懸念される場面ではロバスト統計を併用する」が学会論文の標準作法となっています。

⚠️ ロバスト統計の追加的落とし穴（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 の詳細

❓ よくある質問（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

📖 もっと深掘り：ロバスト統計の知的系譜

ロバスト統計は単なる「外れ値対策」ではなく、 20 世紀後半の統計学の方向転換を象徴する分野です。 古典統計（フィッシャー学派）は「真のモデルが分かっている」前提で最尤推定を最適化したのに対し、 ロバスト統計は「モデルは近似でしかない」前提で小さな逸脱に対する安定性を最適化します。 これは制御工学の H∞ 制御や機械学習の adversarial robustness と通底する考え方であり、 「最悪ケース思考」の数理的表現と言えます。

📅 詳細年表（拡張版）

🤔 なぜ Huber は「最適」なのか

Huber (1964) の重要な結果は、 「$\epsilon$ 汚染のクラス $\{F : F = (1-\epsilon)\Phi + \epsilon H, H \in \mathcal{F}\}$ に対して、 漸近分散の最悪値を最小化する推定量は M 推定量である」というものです。 つまり Huber 推定量は「想定外の汚染への被害を最小化する」という意味で最適です。 これは制御理論の minimax 思考、 経済学の Gilboa-Schmeidler maxmin expected utility、 機械学習の adversarial training と完全に同じ哲学であり、 「最悪ケース対策」という統一的視点で理解できます。

💼 産業ユースケース・深掘り 10 選

ロバスト統計の現場応用は驚くほど広い。 ここでは具体的な技術選択とその理由まで踏み込みます。

❓ 追加 FAQ — 実務で迷いがちな 12 問

✅ 実務チェックリスト — ロバスト統計を導入する 12 ステップ

⚠ よくある誤用 8 選

🧪 モンテカルロ実験 — 汚染下での挙動を実感する

教科書の数式を眺めても直感は育ちません。 シミュレーションで「平均がいかに脆く、 中央値・MAD がいかに頑健か」を体験してみましょう。 以下のコードを実行すると、 5% 汚染下での各推定量の MSE が一覧できます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# 汚染下でのモンテカルロ比較
import numpy as np
from scipy import stats

def simulate(n=100, eps=0.05, contamination_scale=10, n_sim=2000):
    results = {'mean':[],'median':[],'trim10':[],'huber':[],'mad':[]}
    for _ in range(n_sim):
        # (1-ε)N(0,1) + εN(0, σ_c²) の混合
        clean = np.random.standard_normal(int(n*(1-eps)))
        contam = np.random.standard_normal(int(n*eps)) * contamination_scale
        x = np.concatenate([clean, contam])
        results['mean'].append(np.mean(x))
        results['median'].append(np.median(x))
        results['trim10'].append(stats.trim_mean(x, 0.10))
        # Huber location
        try:
            huber = stats.trim_mean(x, 0.20)  # 簡易近似
            results['huber'].append(huber)
        except: pass
        results['mad'].append(stats.median_abs_deviation(x, scale='normal'))
    return {k: (np.mean(v), np.std(v)) for k,v in results.items() if v}

# 汚染 0%, 5%, 10%, 20% で比較
for eps in [0.0, 0.05, 0.10, 0.20]:
    out = simulate(eps=eps)
    print(f'ε={eps:.2f}: ', {k: f'{m:.3f}±{s:.3f}' for k,(m,s) in out.items()})

典型的出力：

平均の SD は汚染率に比例して肥大化（0% で 0.10 → 20% で 0.45）。 一方、 中央値とトリム平均の SD はほぼ一定。 これが「破壊点」と「効率」の両立を表す具体的データです。

📈 漸近相対効率 (ARE) 表

📔 ミニ用語集 — ロバスト統計の専門語 30

Breakdown point（破壊点）

推定量が「制御不能になる」までに耐えられる汚染データの割合。 平均は 0%、 中央値は 50%。

Influence Function (IF)

1 つのデータ点 $x$ が推定量に与える影響を表す関数。 有界なら推定量はロバスト。

Gross Error Sensitivity (γ*)

IF の絶対値の上限。 小さいほど頑健。

Local Shift Sensitivity (λ*)

IF の傾きの上限。 観測値のわずかな変化への感度。

Rejection Point (ρ*)

IF が 0 になる $x$ の最小値。 再下降型推定量で重要。

M-estimator

Maximum-likelihood-like estimator。 $\sum_i \psi(x_i - \theta) = 0$ の解。

L-estimator

順序統計量の線形結合（例：トリム平均、 ウィンザライズ平均、 中央値）。

R-estimator

順位検定統計量から導かれる推定量（例：Hodges-Lehmann）。

S-estimator

残差のロバスト尺度を最小化する回帰推定量。 破壊点 50% 達成可能。

MM-estimator

S 推定量を初期値として M 推定を実行。 破壊点 50% + 効率 95% の両立。

Huber loss

絶対誤差と二乗誤差をハイブリッドした損失関数。 微分可能で凸。

Tukey biweight

再下降型 ψ 関数。 巨大な外れ値の重みを完全に 0 にする。

Hampel function

3 つの定数 $a, b, c$ で定義される再下降型 ψ 関数。 段階的に重みを落とす。

IRLS

Iteratively Reweighted Least Squares。 M 推定量を解く標準アルゴリズム。

LMS (Least Median of Squares)

残差の中央値を最小化する回帰。 破壊点 50% だが効率は低い。

LTS (Least Trimmed Squares)

残差の小さい方 $h$ 個の二乗和を最小化。 LMS より効率が高い。

MCD (Minimum Covariance Determinant)

多変量ロバスト推定。 共分散行列式が最小となる $h$ 点を選ぶ。

MVE (Minimum Volume Ellipsoid)

$h$ 点を含む最小体積の楕円。 MCD の前身。

OGK (Orthogonalized Gnanadesikan-Kettenring)

多変量ロバスト共分散の高速近似アルゴリズム。

Trimmed mean (トリム平均)

両端 α% を除外して平均。 オリンピック採点方式と同じ発想。

Winsorized mean (ウィンザライズ平均)

両端 α% を最近接値で置換してから平均。 情報損失をトリムより抑える。

Hodges-Lehmann estimator

$\mathrm{median}_{i \le j}\{(x_i+x_j)/2\}$。 R 推定量の代表。

Qn estimator

$\mathrm{quantile}_h\{|x_i-x_j|\}$。 MAD より効率が高いロバスト尺度。

Sn estimator

$\mathrm{median}_i\{\mathrm{median}_j|x_i-x_j|\}$。 効率 58% で MAD (37%) より高い。

Contamination neighborhood

$\{F : (1-\epsilon)F_0 + \epsilon H, H \in \mathcal{F}\}$。 Huber の最悪ケース設定。

Asymptotic Relative Efficiency (ARE)

「正規時の効率の比」。 中央値の vs 平均は $2/\pi \approx 0.637$。

Masking

外れ値どうしが互いを「典型」に見せかける現象。 多変量で深刻。

Swamping

外れ値が他の点を「外れ値」と誤判定させる現象。

Leverage point

説明変数空間で外れた点。 回帰係数を大きく動かす。 OLS は無防備。

Adversarial robustness

機械学習版「最悪ケース汚染への耐性」。 Huber 哲学の現代版。

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

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

💼 実務での応用シナリオ

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

本サイトの 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 を求める実践
• 箱ひげ図ページ： ロバスト統計の図的表現

📚 さらに学ぶには

このサイト内

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

推奨書籍・教材

• 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%）。 現代ロバスト統計の到達点です。

📚 参考文献・関連教材（拡張）

• Huber, P. J. (1964). "Robust Estimation of a Location Parameter." Annals of Mathematical Statistics 35: 73-101. — M 推定量の原論文。
• Hampel, F. R. (1974). "The Influence Curve and Its Role in Robust Estimation." Journal of the American Statistical Association 69: 383-393. — 影響関数の定式化。
• Rousseeuw, P. J. (1984). "Least Median of Squares Regression." JASA 79: 871-880. — 破壊点 50% 回帰の出発点。
• Yohai, V. J. (1987). "High Breakdown-Point and High Efficiency Robust Estimates for Regression." Annals of Statistics 15: 642-656. — MM 推定量。
• Maronna, R. A., Martin, R. D., Yohai, V. J. (2006). "Robust Statistics: Theory and Methods." Wiley. — 現代の決定版教科書。
• Wilcox, R. R. (2017). "Introduction to Robust Estimation and Hypothesis Testing." 4th ed. Academic Press. — 実務向け教科書。
• Madry, A. et al. (2018). "Towards Deep Learning Models Resistant to Adversarial Attacks." ICLR. — 深層学習で minimax robust。
• scikit-learn ドキュメント: HuberRegressor, RANSACRegressor, TheilSenRegressor。
• statsmodels: RLM (Robust Linear Models), QuantReg (分位点回帰)。
• R: robustbase パッケージ（MM 推定量、 MCD、 OGK、 S 推定量）が最強。

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

ロバスト統計（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）、 さらに「外れ値検出」「敵対的学習に対する頑健性」の研究へと発展しています。

🔗 隣接手法への橋渡し

「ロバスト統計」は単独で完結する手法ではなく、 隣接領域と連携することで真価を発揮する。

• 上流: 箱ひげ図 や ヒストグラム で 外れ値 を発見した段で、 平均・分散の代替指標として導入する
• 並列: 中央値・四分位・MAD・トリム平均など。 外れ値の影響を 50% 受けない頑健性
• 下流: ノンパラメトリック検定 や 分位点回帰 など、 分布仮定を弱める高度手法に発展

外れ値で平均が引きずられる現象は実務で頻発し、 「真の中心」を守る ロバスト統計はデータが綺麗でない現実世界でこそ威力を発揮する。

🌳 決定木 — どのロバスト手法を選ぶか

🌳 決定木 — どのロバスト手法を選ぶか

「ロバスト統計を使え」と言われても、 具体的にどれを？ 以下の決定木で迷いを減らしましょう。

📘 Python ライブラリ早見表

🗣 ナラティブ — 「平均値の独裁」を疑え

「日本人の平均年収は 458 万円」「全世帯の平均貯蓄額は 1791 万円」――こうした平均値ベースの政府統計は、 一部の超富裕層に強く引きずられます。 中央値で語れば「貯蓄額の中央値は 1104 万円」となり、 印象がガラリと変わる。 これは技術的問題ではなく、 「典型を語るか、 合計を語るか」という政策コミュニケーションの問題です。 ロバスト統計は単なる手法集ではなく、 「平均値がいつも正しい答えとは限らない」という認識を共有する道具でもあります。

SSDSE-B-2026 で「医師数の全国平均」を語るとき、 東京が分子・分母の両方を大きく動かします。 「全国の医師数の中央値」を併記すれば、 「典型県と東京の差」が浮き彫りになる。 政策議論で「平均だけ」を出すのは現代的にはミスリードと言っても過言ではありません。 中央値、 MAD、 四分位範囲を併記する文化を作ることが、 統計家としての社会的責任の一部です。

機械学習時代になり、 「ロバスト」の概念はadversarial robustness（敵対的攻撃への耐性）として再注目されています。 古典統計の minimax 理論（Huber 1964）は、 半世紀を経て深層学習のセキュリティ研究で蘇生しました。 「最悪ケースに備える」という統計学の知恵は、 機械学習が普及した社会においてセキュリティ・公正性・信頼性の根幹を支える概念になりつつあります。 ロバスト統計を学ぶことは、 単に外れ値を扱う技術を得るだけでなく、 「不確実な世界でどう意思決定するか」という普遍的な問いに向き合うことでもあります。

「ロバスト統計 (Robust Statistics)」を扱う際の手法選択は、 状況に応じて以下のフローで判断すると迷いが減る。

1. Step 1: 目的は記述か予測か?
   • 記述 (現状把握・要約) → 集計・可視化・要約統計量で全体像を掴む
   • 予測 (未知データへの推定) → モデル構築・検証フェーズへ移行
2. Step 2: データの種類・規模は?
   • 数値データ・小〜中規模 → 「ロバスト統計 (Robust Statistics)」やその拡張手法を直接適用
   • カテゴリデータ → 中央値・IQR 等の順序統計量・ノンパラメトリック手法と組み合わせ
   • 大規模・高次元 → Huber 損失等のロバスト損失関数・RANSAC 等の外れ値耐性手法を選択
3. Step 3: 結果の解釈・共有は?
   • 専門家向け → 数値指標・統計検定で精緻に評価
   • 非専門家向け → 可視化・自然言語での要約を重視

このフローに沿って判断することで、 「ロバスト統計 (Robust Statistics)」を中核とした適切な手法選択ができる。

🎮 触って理解する

数直線上の 10 個のデータ点のうち、 オレンジの点（外れ値候補）をドラッグして極端な位置へ動かしてみましょう（マウスでもタッチでも可）。 平均 (mean) と 標準偏差 (s) は大きく引きずられますが、 中央値 (median)・MAD・IQR はほとんど動きません。 これが ブレークダウンポイント（破壊点） の違いを体感する最短ルートです。

💡 直感 — ブレークダウンポイント（破壊点）

ブレークダウンポイントとは 「推定量を無限大まで壊すのに必要な汚染データの最小割合」です。 平均は たった 1 点を無限遠へ動かせば一緒に無限へ飛ぶので BP = 0%（＝ 0 に近い脆さ）。 中央値は半分のデータを壊さない限り有限に留まるので BP = 50%（理論上限）。 上の数直線で、 オレンジ点を右端まで引っ張っても 中央値の縦線がほぼ動かないのは、 順位で決まる中央値が「値の大きさ」ではなく「並び順」しか見ていないからです。 10% トリム平均は両端 1 点ずつを捨てるため BP = 10%、 IQR は上下 25% を無視するため実質 BP = 25%。

⚠️ よくある落とし穴

• MAD をそのまま標準偏差と比べない：素の MAD は正規分布下で σ より小さく出ます。 σ 相当にするには ×1.4826（＝ 1/Φ⁻¹(0.75)）を掛けます。 上の表もスケール補正済みの値です。
• 離散データで MAD = 0 が起こる：同値が多い（半数以上が同じ値）と MAD がゼロになり、 z 値の分母が 0 で発散します。 そのときは IQR/1.349 や Sn/Qn 推定量へ切り替えます。
• トリム平均は「捨てる割合」を事前固定：汚染率より小さいトリム率だと外れ値が残り、 大きすぎると効率を無駄に失います。 汚染率が不明なら中央値か M 推定が無難。
• 「頑健＝常に正しい」ではない：外れ値が本物の重要情報（不正検知・災害・希少事象）のとき、 頑健統計は肝心の信号を無視します。 何を「外れ値」と見なすかは分析目的の判断であり、 統計だけでは決められません。
• 中央値は「動かない」だけで「効率的」ではない：正規分布下の中央値の漸近効率は約 64%（2/π）で、 きれいなデータでは平均より分散が大きい推定になります。 頑健性と効率はトレードオフです。

🚀 発展 — Huber 推定・M 推定

平均は残差二乗和 $\sum (x_i-\mu)^2$ を最小化する M 推定（$\rho(u)=u^2$）、 中央値は絶対値和 $\sum|x_i-\mu|$ を最小化する M 推定（$\rho(u)=|u|$）と見なせます。 M 推定量（Maximum-likelihood-type estimator）はこの $\rho$ 関数を一般化し、 「中心は二乗（＝効率的）、 外れ値領域は緩やかに（＝頑健）」と切り替えます。 その代表が Huber 損失：

さらに外れ値を完全に無視したい場合は、 影響関数が遠方でゼロへ戻る Tukey biweight（再降下 M 推定）を使います。 実務では statsmodels.api.RLM（Huber / Tukey / Hampel を選択可）や sklearn.linear_model.HuberRegressor で回帰に、 位置推定には scipy.stats の各種関数が使えます。 より深く知るには 損失関数、 平均、 中央値、 標準偏差、 MAD、 四分位範囲 (IQR)、 外れ値 の各ページも参照してください（M 推定・Huber・破壊点・トリム平均は本ページ内で解説）。

🎯 解説深化 — 「外れ値検出の道具」としてのロバスト統計：自己マスキングとロバスト z スコア

本文では「外れ値があっても中心の推定が壊れない」ことを主題にした。 この深化セクションでは視点を反転させ、 ロバスト推定量を「外れ値を見つけるものさし」として使うとき何が起きるか を扱う。 鍵になるのは、 古典的な z スコアが持つ構造的欠陥 ── 外れ値が自分自身を隠す「自己マスキング」 である。

🧭 直感 — ものさし自体が外れ値で伸びる

外れ値検出の定番は z スコア $z_i = (x_i - \bar{x})/s$ で「$|z|>2$ なら疑え」という手順である。 ところがこの式の分母 $s$（標準偏差）と分子の $\bar{x}$（平均）は、 どちらも BP = 0% の非頑健な統計量。 つまり外れ値が混ざると、 平均は外れ値の方向へ引き寄せられ、 標準偏差は膨張する。 測られる側の外れ値が、 測る側のものさしを自分に有利なように変形させる ── これが自己マスキングである。 対策は分子・分母を BP = 50% の部品に交換した ロバスト z スコア（Iglewicz & Hoaglin 1993 の修正 z スコア）：

SSDSE-B-2026（2023 年・47 都道府県）の総人口 A1101 で両者を実測比較すると、 検出結果がまるで違う：

東京の総人口 1,408.6 万人は、 古典 z では 4.09（平均 264.6 万・SD 279.8 万で標準化）だが、 ロバスト z では 13.57（中央値 154.9 万・1.4826·MAD = 92.4 万で標準化）。 東京自身が SD を 280 万まで膨らませているせいで、 古典 z は東京の「外れっぷり」を 3 分の 1 以下に過小報告している。 転入者数ではさらに極端で、 東京の z = 5.11 に対しロバスト z = 23.95。 中位県の約 24 個分の「ものさし」の彼方にいることが、 頑健な分母を使って初めて数値化される。

⚠️ 落とし穴（重要）— 3 つの構造的な罠

罠 1：古典 z スコアには数学的な天井がある。 標本 z スコアは恒等的に $|z_i| \le (n-1)/\sqrt{n}$ を超えられない。 n = 47 なら上限は 6.71。 つまり都道府県データでは、 どんなに常軌を逸した値が 1 つ混ざっても z は 6.71 止まりであり、 「$|z|>7$ を外れ値とする」ような基準は原理的に 1 件もヒットしない。 外れ値が極端になるほど SD が膨らみ z が頭打ちになる ── 検出器としての z スコアは小標本で構造的に無力化する。 ロバスト z にはこの天井が存在しない（分母が汚染に反応しないため）。

罠 2：MAD = 0 の崩壊。 データの過半数が同一値だと中央値との偏差の過半数が 0 になり、 MAD = 0 → ロバスト z が 0 割りで発散する。 架空の合成例：市町村別の「政令指定都市数」のようなカウントデータで {0,0,0,0,0,0,0,1,2} の 9 件を考えると、 中央値 0・MAD 0 となり、 0 以外のすべての値が「無限大の外れ値」と判定されてしまう。 離散・ゼロ過剰データでは MAD の代わりに $Q_n$/$S_n$ 推定量（Rousseeuw & Croux 1993）や IQR ベースの基準に切り替える必要がある（この段落の数値は説明用の架空例であり SSDSE の実測値ではない）。

罠 3：「9 県がフラグ＝ 9 県がエラー」ではない。 上の実測でロバスト z が総人口 9 県を検出したが、 これは測定ミスではなく日本の人口分布が右に裾を引く構造そのものを拾っている。 ロバスト z の閾値 3.5 は分布が概ね対称であることを暗黙に仮定しており、 歪んだ分布では片側だけ大量にフラグが立つ。 検出は「調査対象リストの作成」であって「削除リストの作成」ではない ── フラグ後の原因究明（測定誤差か、 実在する構造か）は本文の落とし穴①と同じく分析者の責務である。

🚀 発展 — 歪みを織り込む次の一手

「対称前提の閾値」を超える現代的な拡張が 2 系統ある。 (1) medcouple による調整箱ひげ図（Hubert & Vandervieren 2008）：分布の歪度をロバストに測る medcouple (MC) でフェンスを非対称に伸縮させ（右歪みなら上側フェンスを $1.5\,e^{3\,MC}\cdot$IQR に拡張）、 右裾の正常な値を誤検出しない。 Python では statsmodels.stats.stattools.medcouple で計算できる。 (2) Rousseeuw & Croux の $Q_n$・$S_n$：中央値を経由せずペア間距離の順序統計量からスケールを推定する。 BP = 50% を保ちながら正規時効率が MAD の 37% → $Q_n$ で 82% に改善し、 非対称分布でも自然に機能する（scipy.stats.iqr や R の robustbase::Qn）。 「中央値 + MAD」は入口であり、 スケール推定量の選択自体が研究分野として続いている ── という視界を持つと、 本文の M 推定・MM 推定の系譜と外れ値検出の系譜が「頑健なものさし作り」という同じ幹から分岐していることが見えてくる。

🔗 関連ページ

• 外れ値 — 検出フラグを立てた後の「原因究明」の作法はこちら
• 中央値 ・ MAD — ロバスト z スコアの分子・分母を構成する部品
• 標準偏差 ・ 標準化 — 古典 z スコアの構成要素と、 それが膨張する仕組み
• 四分位範囲 (IQR) ・ 箱ひげ図 — 1.5×IQR フェンスと medcouple 調整版の出発点
• z スコアそのものの単独ページは現在未収録（本セクションの定義式を参照）