外れ値がある場合にもいい感じにしてくれる反復再重み付け最小二乗法(Iteratively Reweighted Least Squares method)がどういう原理で何をしているのかを勉強した痕跡です。

普通の最小二乗法ではいけないのか

外れ値に弱いから。

通常の最小二乗法でも，データが精密であればそれなりの結果は得られます。
(下図の場合，真の値$y=2x+1$に対して$y=1.94x+1.39$となります)

ソースコード

import matplotlib.pyplot as plt
import numpy as np
import random as rnd

@np.vectorize
def func(x, a, b):
  return a*x + b

X = np.linspace(0, 10, 20)
Y = [func(x, 2, 1) + rnd.normalvariate(0, 1) for x in X]
plt.scatter(X, Y)

coeffs = np.polyfit(X, Y, 1)
print(coeffs)
t = [0, 10]
plt.plot(t, func(t, *coeffs))

plt.xlim(xmin=0, xmax=10)
plt.show()

この方法は簡単でいいのですが，外れ値が入ってしまうとそれに引きずられやすいという重大な欠点があります。
(上と同条件で$y=1.59x+3.60$)

一応ソースコード

import matplotlib.pyplot as plt
import numpy as np
import random as rnd

@np.vectorize
def func(x, a, b):
  return a*x + b

X = np.linspace(0, 10, 20)
Y = [func(x, 2, 1) + rnd.normalvariate(0, 1) for x in X]
X = np.append(X, .5)
Y.append(20)
plt.scatter(X, Y)

coeffs = np.polyfit(X, Y, 1)
print(coeffs)
t = [0, 10]
plt.plot(t, func(t, *coeffs))

plt.xlim(xmin=0, xmax=10)
plt.show()

ではどうしましょうかと言うと，外れ値を手動で取り除くということが考えられます。
上の例のように明らかであればいいのですが，際どい点になってくると恣意的な作業になってしまい，客観性が損なわれます。
実際，例では左上の$(0.5, 20)$のみが意図的に加えた外れ値ですが，図のように線が引いてあると$(9, 20)$付近の線より上にある点も外れ値として除外する人がいるかもしれません。

反復再重み付け最小二乗法

データに重みを付ける

極端な外れ値の影響を受けてしまう原因として，データに対する仮定が挙げられます。
通常の最小二乗法(OLS)ではすべてのデータの誤差は同一の正規分布$N(0, \sigma^2)$に従うと仮定します。
すなわち，正常なデータと外れ値の誤差を対等な重みで評価しているため，外れ値が結果に影響を与えることになります。

このような影響を排除するには，各データに重みを付ける必要があります。
OLSでは残差平方和を
$$
s=\sum_i(y_i-\hat y_i)^2
$$
と評価しますが，これに対して
$$
s_w = \sum_iw_i(y_i-\hat y_i)^2
$$
というように重みを付けることで，外れ値の影響を抑えることができます。
ここで，$y_i$はデータ，$\hat y_i$は推定値を表しています。

重みの決定方法としては，測定誤差の分散$\sigma^2$に対して$w_i=1/\sigma_i^2$とする方法が考えられます。
また，各データが$k_i$回の測定結果の平均値である場合に，$w_i=k_i$とすることも考えられます。
(より多くの測定に基づくデータをより信頼できるデータとして評価する)

ロバスト最小二乗法

しかし，実際にデータを見て適切に重みの付け方を決めることができるとは限りませんし，重みの決め方に恣意性が生じては元も子もありません。
そこで，この過程も機械的に処理する方法として，ロバストな最小二乗法，反復再重み付け最小二乗法(以下「IRLS」)が利用されます。

ここでは，説明変数$x_{ij}$によってあらわされるデザイン行列
$$
X=\begin{bmatrix}
1 & x_{11} & x_{12} & \dots & x_{1k} \\
1 & x_{21} & x_{22} & \dots & x_{2k} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n1}　& x_{n2} & \dots & x_{nk}
\end{bmatrix}
$$
と観測値ベクトル$\boldsymbol{y}$について，
$$
\boldsymbol{y}=X\boldsymbol\beta+\boldsymbol\varepsilon
$$
となるパラメータベクトル$\boldsymbol\beta$を求めることとします。
ただし，$\boldsymbol\varepsilon$は誤差ベクトルを表します。

パラメータの推定値$\hat{\boldsymbol\beta}$に対して応答変数の推定値を$\hat{\boldsymbol{y}}=X\hat{\boldsymbol\beta}$とするとき，残差ベクトルは$\boldsymbol{e}=\boldsymbol{y}-X\hat{\boldsymbol\beta}$と定義されます。
誤差ベクトルとデザイン行列が独立であるとき，$X^\mathsf{T}{\boldsymbol e}=\boldsymbol0$より
$$
X^\mathsf{T}(\boldsymbol{y}-X\hat{\boldsymbol\beta})=\boldsymbol0
$$
と書けます。
これより
$$
X^\mathsf{T}X\hat{\boldsymbol\beta}=X^\mathsf{T}\boldsymbol{y}
$$
となり，$(X^\mathsf{T}X)^{-1}$が存在するとき
$$
\begin{array}{l}
\hat{\boldsymbol\beta}=(X^\mathsf{T}X)^{-1}X^\mathsf{T}\boldsymbol{y}\\
\hat{\boldsymbol{y}}=X\hat{\boldsymbol\beta}=X(X^\mathsf{T}X)^{-1}X^\mathsf{T}\boldsymbol{y}
\equiv H\boldsymbol{y}
\end{array}
$$
を得ます。
ここで定義した$H=X(X^\mathsf{T}X)^{-1}X^\mathsf{T}$はハット行列ないし射影行列と呼ばれます¹。

データに重みを付ける場合には，対角成分に各データの重みを持つ行列$W$を用いて，$\boldsymbol\beta$の推定値は
$$
\hat{\boldsymbol\beta}=(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W\boldsymbol{y}
$$
として求められます。
ここで，$\hat{\boldsymbol\beta}$の計算には$W$が必要ですが，$W$を求めるためには$\hat{\boldsymbol\beta}$が既知である必要があります。
当然のことながらこれでは問題が解決しないので，反復計算により$\hat{\boldsymbol\beta}$を計算することになります。
そのためのパラメータ(及び重み)を計算するためのアルゴリズムがIRLSです。

IRLSでは，以下のアルゴリズムでパラメータを決定します。

モデルによる近似を行い，パラメータを得る (初回はOLSの結果などを用いる)
調整済み残差を計算し，標準化する
各データの重みを決定する
1に戻る

近似が収束したら終了になります。

調整済み残差

調整済み残差は，
$$
r_{adj}=\frac{r_i}{\sqrt{1-h_i}}=\frac{(y_i-\hat y_i)^2}{\sqrt{1-h_i}}
$$
として計算されます。
ここで，$h_i$はてこ比です。
てこ比は，ハット行列の対角成分$h_{ii}$であり，
$$
\left\{\begin{array}{l}
0\leq h_{ii}\leq1\\
\sum_{i}h_{ii}=p
\end{array}\right.
$$
という条件を満たします。
ただし，$p$は回帰モデルの係数の数を表します。
てこ比は入力空間で特定の観測位置の回帰予測に対する影響を評価します。
入力空間の中心から離れるほどてこ比は大きくなるため，てこ比が大きい点は外れ値であると考えられます。
$r_{adj}$の計算ではこれを反映して$h_i$が大きい点の影響を小さくするようにしています。

残差の標準化

標準化された残差は
$$
u=\frac{r_{adj}}{Ks}
$$
と定義されます。
ここで，$K$は調整定数，$s=MAD/0.6745$は誤差項の標準偏差の推定値を表します。
調整定数については後述します。

誤差項の標準偏差の推定値

$$
s=\frac{MAD}{0.6745}
$$
において，$MAD$はmedian absolute deviation(中央絶対偏差)であり²，
$$
MAD = \mathrm{median}\left(|x_i - \tilde{x}|\right)
$$
と定義されます($\tilde{x}$はデータの中央値)。
MADは標準偏差の推定に利用され，(標準偏差と比べて)比較的ロバストであるとされています。
標準偏差との関係は，累積分布関数$\Phi$を用いて
$$
MAD=\Phi^{-1}(0.75)\sigma
$$
であり³，正規分布を仮定したうえでこれを解いて標準偏差の推定値$\hat\sigma\approx MAD/0.6745$を得ます。

重みの決定

標準化した残差を基にして各データの重みを決定します。
重みを決定する関数として様々なものが提案されていますが，代表的なものにTukeyのbiweight関数
$$
w(u)=\begin{cases}
(1-u^2)^2 & |u|\leq1 \\
0 & |u|>0
\end{cases}
$$
があります。
このときにこの関数のロバスト性を決定するパラメータが上述の調整定数$K$であり，$K$を大きな値にすると$u$は小さくなるため各データの重みが大きくなり，結果としてロバスト性が低下することになります。
MATLABのrobustfitなどの既定値は4.685(漸近有効性が0.95になる値)となっています。

以上でめでたく重み$w_i$を決定できたため，対角成分にこれらの重みを持つ行列$W$を用いてパラメータベクトルの推定値
$$
\hat{\boldsymbol\beta}=(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W\boldsymbol{y}
$$
を計算することができます。
得られた推定値を用いて再度重みを計算し，推定値を計算する……ということを繰り返し行うことで適当な値に収束し，パラメータベクトルの決定は収束となります。

参考文献

推定値にハットが付いているのでハット行列，観測値ベクトルを推定値ベクトルに射影するので射影行列と呼ばれます。 ↩
meanもmedianもmodeもすべてmであるため，中央値であることを分かりやすくするためにMedADと表記している文献もあります。 ↩
確率誤差の定義より，平均に対して対称かつ面積が2分の1となる範囲は左側25%点～右側25%となるので，0.75を使用します。 ↩