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

# はじめに
#### 本記事で分かること
* 最小二乗法の概要
* 最小二乗法のMATLABを用いた実装プログラムと出力結果の例
* 最小距離二乗法の概要
* 最小距離二乗法を導出する数式
* 最小距離二乗法のMATLABを用いた実装プログラムと出力結果の例

#### 参考文献について
最小距離二乗法という名は、おそらく正式名称ではありませんが、

 [“最小距離2乗法による回帰直線の求め方”（永島弘文，大館孝幸，荒井太紀雄，1986）](https://www.jstage.jst.go.jp/article/jject1975/12/1/12_1_51/_article/-char/ja/)
上記文献にて、この名称の記載があり、名称および解法についてこちらを参考にさせていただきました。

この方法は，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次関数を求める」** というケースが多いでしょう。

![最小二乗法の例.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/326d7b5d-a5b4-b9aa-7ad5-80dbde716f12.jpeg "最小二乗法の例")

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

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

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

![最小距離二乗法の例.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/92aa3a73-a4d8-c8be-a2cb-331568d26be5.jpeg "最小距離二乗法の例")

# 最小二乗法とは
[最小二乗法とは？（Wikipedia）](https://ja.wikipedia.org/wiki/%E6%9C%80%E5%B0%8F%E4%BA%8C%E4%B9%97%E6%B3%95)
一般的な最小二乗法においては、Ｙ軸の値に誤差が含まれます。
例えば、Ｘ軸は時間や個数、人数などであり、Y軸は測定値であることが多いでしょう。
したがって、**最小二乗法では、各点と1次関数のY軸方向の距離を最小化する１次関数**を見つけます。


それでは、前章の最小二乗法の例について、MATLABを使って実践してみます。
まずは、適当なプロットを用意します。今回は以下のような8個の点を用意します。
```matlab
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');
```
![plots.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/ab164918-6c60-943b-a85b-0d409c45ac19.jpeg "プロット")

次に、以下のような処理を施し、目標の1次関数の傾き$u(1)$とY切片$u(2)$を算出します。
```matlab
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)
```
![最小二乗法の例 補助線なし.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/c5f40409-ddad-4597-f751-72009955d42d.jpeg "最小二乗法の出力結果")

できました。

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

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

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


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

![点i.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/a5ea8f66-5b6e-1b87-4023-17ca23e916bb.png)

次に、点$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)$

![di.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/bbdd65ad-ab8c-4a3f-7ac1-3daaca2bdafd.png)

さて、点の数は$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$を算出します。

```matlab
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)
```

![最小距離二乗法の例 補助線なし.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/af0a8c4f-a874-2065-f774-5b360e2180a4.jpeg "最小距離二乗法の出力結果")

できました。


最後に、両方の出力結果を比較してみましょう。
下図の黄色線が最小距離二乗法によって求められた1次関数で、
赤線色が最小二乗法によって求められた1次関数になります。
![出力結果.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/598885/748893f2-4db9-96be-4b7b-2d701d2c17ea.jpeg "出力結果 比較")

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

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