17. ロバスト回帰：外れ値に強い推定

OLSの外れ値脆弱性とロバスト回帰の動機

OLS（通常最小二乗法）の目的関数は残差の二乗和の最小化として定義されます。

$$\hat{\beta}_{\mathrm{OLS}} = \arg\min_{\beta} \sum_{i=1}^{n} (y_i – x_i^\top \beta)^2$$

二乗損失は残差の絶対値に対して二次的に増加するため、大きな残差をもつ観測値が総損失を支配します。この性質により、OLS推定量は単一の外れ値によって大きく歪む可能性があります。

単一観測値が推定量に与える影響はCook距離によって定量化されます。

$$D_j = \frac{(\hat{\beta} – \hat{\beta}_{(j)})^\top (X^\top X)(\hat{\beta} – \hat{\beta}_{(j)})}{p\,\hat{\sigma}^2}$$

ここで$\hat{\beta}_{(j)}$は第$j$観測値を除いた推定量、$p$はパラメータ数です。$D_j$が大きいほど観測値$j$がOLS推定に強く影響します。

外れ値には複数の類型があります。縦方向外れ値は$y$方向に真の関係から乖離した観測値です。高レバレッジ点は説明変数空間で孤立した位置にある観測値であり、回帰直線を強く引き寄せますが、残差が小さく現れることがあるため通常の残差診断では検出が困難です。損失関数の形状が外れ値感度を決定し、二乗損失の非有界な増加がOLS推定量の非ロバスト性の根本原因です。

外れ値の定義は分析文脈に依存します。測定誤差・異常観測・分布の裾の重さによって解釈が異なるため、外れ値の機械的な除外は科学的に不適切です。

（Fig1. 外れ値存在下でのOLS回帰直線とHuber M推定回帰直線の比較）

影響関数と崩壊点

推定量のロバスト性を定量評価するための理論的枠組みとして、影響関数と崩壊点があります。

影響関数は推定量$T$の分布$F$に対する局所感度を測定します。確率的汚染$F_\epsilon = (1-\epsilon)F + \epsilon\delta_x$（$\delta_x$は点$x$にマスを置くディラック分布）に対して次式で定義されます。

$$\mathrm{IF}(x;\,T,\,F) = \lim_{\epsilon \to 0} \frac{T(F_\epsilon) – T(F)}{\epsilon}$$

OLS推定量の影響関数は$x$に対して非有界であり、単一の外れ値が推定値を無制限に歪める可能性を示します。

有限標本崩壊点は、推定量を任意に大きな値へ移動させることのできる汚染観測値の最大割合として定義されます。

$$\epsilon_n^* = \min\left\{\frac{m}{n} \,:\, \sup_{\mathcal{Z}_m} \|\hat{\theta}(\mathcal{Z}_m)\| = \infty\right\}$$

ここで$\mathcal{Z}_m$は元のデータから$m$個の観測値を任意の値に置換したデータセットです。OLS推定量の崩壊点は$1/n$であり、標本サイズの増加とともに0に収束します。これは単一の外れ値によってOLS推定量が任意に歪められることを意味します。理論的に達成可能な崩壊点の上限は$\lfloor n/2 \rfloor / n$（概ね50%）です。

masking効果とは、複数の外れ値が集団を形成することで互いの検出を妨害する現象です。外れ値群がクラスタを形成する場合、個別診断では集団的な影響を捕捉できません。

影響関数は漸近論的概念であり、有限標本では近似にとどまります。実際の有限標本感度は影響関数による予測から乖離する場合があります。

M推定：一般化損失関数による定式化

M推定は二乗損失を一般化した損失関数$\rho$による最小化として定式化されます。

$$\hat{\beta}_M = \arg\min_{\beta} \sum_{i=1}^{n} \rho\!\left(\frac{e_i}{\hat{\sigma}}\right)$$

ここで$e_i = y_i – x_i^\top\beta$は残差、$\hat{\sigma}$はスケール推定量です。最小化の一階条件は$\psi = \rho’$（スコア関数）を用いて次式で表されます。

$$\sum_{i=1}^{n} \psi\!\left(\frac{e_i}{\hat{\sigma}}\right) x_i = 0$$

重み関数$w(e) = \psi(e)/e$を定義すると、推定方程式は加重最小二乗問題と等価になります。$\rho$・$\psi$・$w$の三者は損失・スコア・重みのそれぞれに対応し、この関係がIRLSによる数値計算の基礎となります。

Huber損失関数はチューニング定数$k$を用いて次式で定義されます。

$$\rho_H(e) = \begin{cases} \dfrac{1}{2}e^2 & |e| \leq k \\[6pt] k\,|e| – \dfrac{1}{2}k^2 & |e| > k \end{cases}$$

対応するスコア関数は$\psi_H(e) = \min(k,\,\max(-k,\,e))$であり、大きな残差には定数スコアを与えます。Bisquare（Tukey）損失関数はチューニング定数$c$を用いて次式で定義されます。

$$\rho_B(e) = \begin{cases} \dfrac{c^2}{6}\!\left[1 – \!\left(1 – \dfrac{e^2}{c^2}\right)^{\!3}\right] & |e| \leq c \\[6pt] \dfrac{c^2}{6} & |e| > c \end{cases}$$

Bisquare損失は$|e| > c$の外れ値に定数ペナルティを与え、推定への影響を完全に遮断します。

スケール補正にはMAD（中央絶対偏差）が使用されます。

$$\hat{\sigma} = \frac{\mathrm{median}(|e_i – \mathrm{median}(e_j)|)}{0.6745}$$

除数0.6745は正規分布における一致推定のための補正定数です。

チューニング定数の設定はガウス誤差での漸近効率によって決まります。Huber損失での$k=1.345$およびBisquare損失での$c=4.685$は、正規誤差下でOLSの95%の漸近効率を維持するように設定されています。スケール推定量との同時推定が必要であり、MAD以外のスケール推定量を使用する場合には結果への感度確認が求められます。

M推定が偏りのない推定を行うためには、誤差分布の対称性が必要です。また、誤差のスケールが一致推定可能であることも前提となります。高レバレッジ外れ値に対しては標準的なM推定が有効でなく、GM推定（一般化M推定）の適用が必要になります。スケール推定量の選択が最終的な推定結果に影響するため、感度確認が求められます。

（Fig2. OLS・Huber・Bisquare損失関数の形状比較）

IRLS：反復加重最小二乗による数値計算

M推定量の数値解法として、IRLS（反復重み付き最小二乗法）が使用されます。各反復で重み関数

$$w(e) = \frac{\psi(e)}{e}$$

によって重み行列$W^{(t)} = \mathrm{diag}(w(e_1^{(t)}/\hat{\sigma}),\,\ldots,\,w(e_n^{(t)}/\hat{\sigma}))$を更新し、加重OLS更新式

$$\hat{\beta}^{(t+1)} = (X^\top W^{(t)} X)^{-1} X^\top W^{(t)} y$$

によって係数を更新します。収束判定は係数変化量が閾値$\delta$を下回ることで行われます。

$$\max_j \left|\hat{\beta}_j^{(t+1)} – \hat{\beta}_j^{(t)}\right| < \delta$$

IRLSの反復手順は次のとおりです。まず初期推定（LAD推定や中央値回帰が推奨されます）から始め、残差を計算し、残差に基づいて重みを更新した後、加重OLSを解きます。この操作を収束まで繰り返します。初期値にOLSを用いると外れ値汚染の影響を引き継ぐため、外れ値感度の低い推定量を初期値として選択します。各反復において外れ値の重みは自動的に縮小され、最終的に外れ値には最小の重みが割り当てられます。

Huber損失の場合のみ目的関数が凸であり、IRLSによる大域的収束が保証されます。Bisquare損失は非凸であるため局所解に陥る可能性があり、初期値依存性が高くなります。計算コストはOLSの$O(p^2 n)$に反復回数が乗算されます。

（Fig3. IRLS収束時のHuber重みと標準化残差の関係）

高崩壊点推定：LTSとRANSAC

崩壊点50%を達成する手法として、LTS（最小トリム二乗法）とRANSAC（ランダムサンプルコンセンサス）があります。

LTSは残差の順序統計量$e_{(1)}^2 \leq e_{(2)}^2 \leq \cdots \leq e_{(n)}^2$のうち最小の$h$個の二乗和を最小化します。

$$\hat{\beta}_{\mathrm{LTS}} = \arg\min_{\beta} \sum_{i=1}^{h} e_{(i)}^2(\beta)$$

$h \approx 0.75n$の設定が典型的であり、崩壊点50%を達成します。$h$の選択が崩壊点と統計的効率のバランスを決定し、$h$が小さいほど崩壊点は高くなりますが推定効率は低下します。

RANSACの必要反復数$N$は外れ値混入率$\varepsilon$、最小標本サイズ$s$、目標成功確率$p$から次式で決定されます。

$$N \geq \frac{\log(1-p)}{\log(1-(1-\varepsilon)^s)}$$

RANSACはランダムに$s$点を選択してモデルを適合し、残差がしきい値以下の観測値をインライアとして判定します。この操作を$N$回繰り返し、インライア数が最大のモデルを採用します。FastLTSアルゴリズムは確率的探索によって計算量を削減し、実用的な実行時間を実現します。

LTSは組み合わせ最適化問題であり、厳密解の計算量は$O(n^s)$に達するため、FastLTSなどの近似アルゴリズムが実用的です。RANSACはランダム性を含むため結果が再現不可であり、乱数シードの固定が必要です。RANSACでは外れ値混入率$\varepsilon$の事前知識が必要であり、この推定が不正確な場合には必要反復数が大幅に増加します。LTSでは$h$の選択が崩壊点と効率のバランスを決定します。

仮定・診断・推論

ロバスト回帰が有効に機能するためには、クリーンデータが観測値全体の過半数を占めること（汚染率が崩壊点未満であること）、および真のモデルが線形であることが前提となります。

推定後の統計的推論にはサンドウィッチ型分散推定量が使用されます。スケール補正済み残差を$r_i = e_i / \hat{\sigma}$と定義します。$r_i$はMADスケールで正規化された残差であり、以下の分散推定量の構成に使用されます。

$$\hat{V}(\hat{\beta}) = (X^\top \hat{W} X)^{-1}\!\left(\sum_{i=1}^{n} \psi(r_i)^2\, x_i x_i^\top\right)\!(X^\top \hat{W} X)^{-1}$$

一方、外れ値診断のための標準化ロバスト残差スコアはレバレッジ値$h_{ii}$（ハット行列の対角要素）を用いて次式で定義されます。

$$r_i^* = \frac{e_i}{\hat{\sigma}\sqrt{1 – h_{ii}}}$$

$r_i$（スケール補正済み残差）と$r_i^*$（レバレッジ補正済み診断スコア）は異なる量です。前者は分散推定において外れ値のスコアを評価するために用いられ、後者はハット行列によるレバレッジ効果を補正した外れ値診断のための指標です。

外れ値は四類型に分類されます。縦方向外れ値は$y$方向の乖離、高レバレッジ点は$x$空間での孤立、影響点は推定量を実際に大きく移動させる観測値、masking点は複数の外れ値が相互に検出を妨害する観測値です。swamping効果は正常観測値が外れ値と誤認される現象であり、外れ値集団によって回帰直線が歪んだ場合に発生します。ロバスト回帰後の残差診断によってinlier（正常観測）とoutlier（外れ値）を事後的に特定することが推奨されます。

ロバスト推定量は漸近正規性をもつため、信頼区間はサンドウィッチ型分散に基づいて構成されます。ただし漸近的性質は有限標本では近似にとどまります。

ロバスト回帰は外れ値の影響を抑制するものであり、外れ値の発生原因の調査を代替しません。外れ値が真の現象を反映している場合、影響の縮小より発生原因の調査が優先されます。また、汚染率が50%を超えると崩壊点条件が満たされず、すべての手法が破綻します。

品質管理への応用と代替手法との比較

製造ラインの寸法計測データには、測定機器の誤差・設備異常・工程変化点に起因する外れ値が混入します。こうしたデータにOLSでCpk（工程能力指数）を推定すると、外れ値が回帰係数を歪め、工程能力の誤判定につながります。ロバスト回帰を適用することで、外れ値を明示的に除外せずCpkの安定推定が実現します。

ただし、工程変化点（シフト）による系統的な外れ値はロバスト回帰によって抑制・無視されます。このため変化点の検出には、Shewhart管理図やCUSUM（累積和管理図）との併用が不可欠です。

代替手法との適用場面の使い分けを次に整理します。OLSは誤差が正規分布に従い外れ値のない場合に最大効率を発揮しますが外れ値に脆弱です。WLS（加重最小二乗法）は既知の分散不均一性への対応に適しますが、外れ値の重みは分析者が事前に指定する必要があり外れ値に自動適応しません。分位点回帰は条件付き分位数を推定対象とするため外れ値感度が低くなりますが、条件付き平均の推定には対応しません。外れ値除外とロバスト回帰の比較では、測定ミスなど真に無効な観測値に対して除外後のOLSが有効ですが、根拠のない機械的削除は情報損失を招きます。ロバスト回帰は観測値を保持したまま外れ値の影響を縮小します。