はじめに

本記事で分かること

最小二乗法の概要
最小二乗法のMATLABを用いた実装プログラムと出力結果の例
最小距離二乗法の概要
最小距離二乗法を導出する数式
最小距離二乗法のMATLABを用いた実装プログラムと出力結果の例

参考文献について

最小距離二乗法という名は、おそらく正式名称ではありませんが、

“最小距離2乗法による回帰直線の求め方”（永島弘文，大館孝幸，荒井太紀雄，1986）
上記文献にて、この名称の記載があり、名称および解法についてこちらを参考にさせていただきました。

この方法は，Major Axis Regression (MA) と呼ばれるものです。訳としては「主成分回帰」とされることが多いと思います。

類似のものに，Reduces Major Axis Regression (RMA) というものもあります。

http://aoki2.si.gunma-u.ac.jp/R/MA.html
http://aoki2.si.gunma-u.ac.jp/R/RMA.html

本記事の概要

よく用いられる近似の手法として、最小二乗法があります。
簡単に言えば、XY平面上の適当なプロットに対して、近似される関数を求める手法になります。
（以下、すべて1次関数を用いて近似することにします。）

しかし、一般的な最小二乗法では、XかYどちらか一方のデータにのみ誤差が含まれるという仮定があります。
例えば、「Xは誤差を含まず、Yに誤差が含まれており、各点とY軸方向の距離が最小になる1次関数を求める」 というケースが多いでしょう。

それでは、XとYいずれのデータにも誤差が含まれる場合にはどうすれば良いのでしょうか？

ここで最小距離二乗法が1つの答えとなります。
「各点と、X軸方向とY軸方向の両方の距離（すなわち点と直線の距離）が最小になる1次関数を求める」ことができるように、最小二乗法を拡張したものが、最小距離二乗法になります。

これが，「主成分回帰」と言われる所以です。

最小二乗法とは

最小二乗法とは？（Wikipedia）
一般的な最小二乗法においては、Ｙ軸の値に誤差が含まれます。
例えば、Ｘ軸は時間や個数、人数などであり、Y軸は測定値であることが多いでしょう。
したがって、最小二乗法では、各点と1次関数のY軸方向の距離を最小化する１次関数を見つけます。

それでは、前章の最小二乗法の例について、MATLABを使って実践してみます。
まずは、適当なプロットを用意します。今回は以下のような8個の点を用意します。

x = [0 3 6 9 12 15 18 21]';
y = [32.5 35.5 30.5 40 41.9 47 39.6 42.5]';
plot(x, y, 'o');

次に、以下のような処理を施し、目標の1次関数の傾き$u(1)$とY切片$u(2)$を算出します。

A = ones(8,2);
A(:, 1) = x;
v = y;
u = inv(A'*A)*A'*v;

x2 = 0:21; %目標の1次関数のx軸範囲が、用意したxの0～21までであることを示す
y2 = u(1) .* x2 + u(2);
plot(x2, y2)

できました。

最小“距離”ニ乗法とは

さて、最小距離ニ乗法では、Ｙ軸の値のみならず、Ｘ軸の値に含まれる誤差も考慮します。
X軸も、Y軸も測定値であると考えれば良いでしょう。

（余談：ちなみに筆者は、センサから推定される歩行者の連続した位置情報から、進行方向を算出するために本手法を検討しました。）

したがって、最小距離二乗法では、各点と1次関数の距離を最小化する1次関数を見つけます。

本章では、はじめに最小距離二乗法の解の導出について考えてみます。
以下の図のように、点i $(x_i, y_i)$　(ただし、$i$は自然数で、$0<i\leq n$ を満たす) と
$y=ax+b$ の直線を考えます。
（一般的な場合を考えますので、今回は目盛りは無視してください。）

次に、点$i$から直線 $y=ax+b$ に向けて垂線を下ろし、距離$d_i$を求めます。
この計算は$d_i$と$y=ax+b$が直交していることを利用すると導出することができます。

$\displaystyle {d_i}^2 = \frac{1}{1+a^2} (y_i-ax_i-b)^2$　　　$- (1)$

さて、点の数は$n$個ありますから、
最小化したい値は、$n$個の点と直線の距離の2乗の総和になります。
この総和を$D$とすると、

$D = \displaystyle \Sigma{d_i}^2 = \Sigma{\frac{1}{1+a^2} (y_i-ax_i-b)^2}$　　　$-(2)$

$a$と$b$を変数とみなすと、$D$を最小にするためには以下を満たす必要があります。

$\displaystyle \frac{\Delta{D}}{\Delta{a}} = 0, \frac{\Delta{D}}{\Delta{b}} = 0$　　　$-(3)$

式(3)を$a$について解くと、

$I^2a+Ka-I = 0$　　　$-(4)$
$a = \displaystyle \frac{1}{2I}(-K\pm\sqrt{K^2+4I^2})$　　　$-(5)$

ただし、
$I = (\Sigma{x_i})(\Sigma{y_i}) - n\Sigma{x_i}{y_i} = n^2(\bar{x} \bar{y} - \bar{xy})$　　　$-(6)$
$K = n(\Sigma{y_i}^2 - \Sigma{x_i}^2) + (\Sigma{x_i})^2 - (\Sigma{y_i})^2 = n^2(\bar{y^2} - \bar{x^2} + \bar{x}^2 - \bar{y}^2)$　　　$-(7)$

また$b$については以下の通り。
$b = \displaystyle \frac{1}{n}(\Sigma{y_i} - a\Sigma{x_i}) = \bar{y} - a\bar{x}$

以上の計算を実装すれば、入力のプロットに対して、自動的にそれらの近似直線を出力することができます。

それでは、前々章の最小距離二乗法の例について、MATLABを使って実践してみます。
まずは、最小二乗法の実装と同様の8個の点を用意します。（コードおよび図は割愛）

次に、以下のような処理を施し、目標の1次関数$\space y = a \times x + b$の傾き$\space a$とY切片$\space b$を算出します。

n=8; % プロット数

avr_x = mean(x);
avr_y = mean(y);
d = dot(x,y);
avr_d = d/n;
avr_xx = dot(x,x)/n;
avr_yy = dot(y,y)/n;
avr_x2 = avr_x.^2;
avr_y2 = avr_y .^2;

% 計算式の簡略化のために置く
I = (n.^2) .* (avr_x.*avr_y - avr_d);
K = n.^2 .* (avr_yy - avr_xx + avr_x2 - avr_y2);

a = (-K - sqrt(K.^2 + 4.*(I.^2))) / (2.*I);
b = avr_y - a .* avr_x;
% (a, b) = (0.6833767859339025, 31.512043747694022)

x3 = 0:21;
y3 = a .* x3 + b;
plot(x3, y3)

できました。

最後に、両方の出力結果を比較してみましょう。
下図の黄色線が最小距離二乗法によって求められた1次関数で、
赤線色が最小二乗法によって求められた1次関数になります。

さすがに目視では定性的な違いはわかりませんが、
8個程度のデータでも、計算方法によって異なるグラフが出力されることがわかります。

おわりに

今回は、最小二乗法の拡張として、通称“最小距離二乗法”を紹介しました。
データのX軸とY軸の両方に誤差が含まれる場合に、それらのプロットを1次関数を用いて近似したいときに有効です。

最小“距離”二乗法とは？