NumPyで行列　回帰分析その7 最小2乗法｜Pythonで数学を学ぼう！第28回

最小2乗法｜なぜ2乗なのか

データの点と回帰直線上の点の差を2乗した和の最小化を考えるのが最小2乗法ですが、そもそもなぜ2乗なのでしょうか？絶対値や和や積ではダメなのでしょうか？

これまで見てきたように、2乗した結果がシンプルな回帰係数の公式です。回帰直線を図示してみれば、それが近似直線であることも疑いありません。ここまでの説明では「差」と表現してきました。これを誤差とはしませんでした。

身長と体重の例では、各データと回帰直線との差は誤差ではないからです。しかし真の値がある現象に対して行われる測定の場合には、差は誤差となります。「正確に測定できれば100であるが、測定によって98が得られる」ということです。真の値がある場合には「真の値はわからない」けれど、測定誤差を含んだデータは得られます。

18世紀の数学者ガウスは、ゲッティンゲン大学の天文台長を務めるほど天文学者としても活躍しました。1794年に最小2乗法はガウスによって発見されました。惑星ケレスは1801年1月1日にピアッツィによって発見されるも、見失われてしまいました。

そこでガウスは最小2乗法を駆使して惑星ケレスの軌道計算を行い、それをもとに同年12月31日にオルバースが再発見するという出来事がありました。ガウスは天文学以外にも測地学にも通じており、測定につきものの測定誤差を考察していく中で数学的に真の値を推測する研究（誤差論）を始めました。

ガウスは確率論を応用することで、誤差が正規分布（ガウス分布とも呼ばれる）している事実を突き止めました。説明変数xと目的変数yがy=f(x)というモデルが想定した場合、最小2乗法によって得られるモデルと誤差が正規分布することから得られるモデルが一致することが証明されます。

最小2乗法は、残差（誤差）が正規分布で乗っているという仮定から導かれます。ここに2乗する合理性があります。誤差が正規分布することをもとにした推定法を最尤法（さいゆうほう：最も尤もらしい法）といいます。

最尤法とは、モデル上にデータがのる確率が最大化を考える推定法です。現代のデータサイエンスのおおもとは、地球や宇宙の観測を実際に行っていた数学者ガウスによって発見された数学にあるのです。

最小2乗法＝残差（誤差）が正規分布するときの最尤法

ガウス生誕225周年を記念して2002年にガウス賞が設立されました。社会の技術的発展と日常生活に対して優れた数学的貢献をなした研究者に贈られる賞で、栄えある第1回受賞者は伊藤清（確率微分方程式）です。

ガウス賞の受賞メダルの表面にはガウスの肖像、裏面には惑星セレスの軌道を表す曲線と円（惑星セレス）そして正方形が描かれています。この正方形こそ最小2乗法を表したものです。

重回帰分析における最小2乗法

線形回帰モデルの設定のところで、n個のデータおよびm個の説明変数の関係式を立てました。この式はベクトルと行列を用いて次のように表すことができます。

y=Xθ+e　… (1)

ここで、yは目的変数ベクトル、Xは説明変数行列、θは回帰係数ベクトル、そしてeは差ベクトルです。Xはn行m+1列行列、列ベクトルθはm+1行1列行列なので積Xθはn行1列行列すなわちn次元列ベクトルです。

差の2乗和Qをベクトルと行列で表すところからはじめます。差ベクトルe（列ベクトル）に対してその転置行列e'（'で転置行列を表します）は行ベクトルになります。したがって、
Q = e'e
と表すことができます。

次に(1)式より、e=y-Xθとなるので、差の2乗和Qは次のように変形できます。

Q = e'e
= (y-Xθ)'(y-Xθ)
= y'y -2y'Xθ＋θ'X'Xθ

Qをθの関数とみて、θで微分（偏微分）した式を0とすればQが最小となります。
∂Q/∂θ＝-2X'y + 2X'Xθ= 0 … (2)

正規方程式

(2)式から次の式が得られます。これは正規方程式と呼ばれます。
X'Xθ = X'y　… (3)

正規方程式(3)において、行列X'Xに逆行列が存在すれば、差の2乗和Qを最小にする回帰係数ベクトルθは次のように求められます。
θ = (X'X)^-1X'y　… (4)

転置行列

重回帰分析では回帰係数ベクトルを求める際に転置行列の計算が何度も現れます。Pythonで確かめてみましょう。

転置行列の公式

転置行列の公式の正規方程式の導出にも使われます。

(1) (A')=A
(2) (AB)' = B'A'
(3) αとβがスカラーならば
α' = α　(4) (αA+βB)' = αA' + βB'

3×2行列Aと2×3行列Bを用いて、公式(1)と(2)を確かめてみます。

>>>　 A = np.array([[1, 2, 3],
>>>　 [4, 5, 6]])
>>>
>>>　 B = np.array([[-1, 0],
>>>　 [2, 3],
>>>　 [4, 6]])
>>>
>>>　 print((A.T).T)
>>>　 [[1 2 3]
>>>　 [4 5 6]]
>>>
>>>　 print((A@B).T)
>>>　 [[15 30]
>>>　 [24 51]]
>>>
>>>　 print(B.T @ A.T)
>>>　 [[15 30]
>>>　 [24 51]]

偏微分

転置行列とともに正規方程式導出に使われるのが偏微分です。関数の極値を求める場合に用いるのが微分です。例えば、f(x)=x²のときf(x)をxで微分するとf'(x)=2xとなります。この微分はdf/dxとも書かれます。

2つ以上の変数を持つ多変数関数に対する微分が偏微分と呼ばれます。例えば、2変数関数f(x, y) = x²+y³に対して、f(x, y)をxで微分することをxで偏微分するといいます。これを∂f/∂x = 2xと表します。

xで偏微分する場合、変数x以外の変数は定数とみなして微分します。f(x, y)をyで偏微分すると、∂f/∂y = 3y²となります。

最小2乗法で扱う「差の2乗和Q」は2つ以上の回帰係数の多変数関数です。したがって、Qの極値を求める場合には、Qを説明変数ごとに偏微分して0とおきます。

はたして、連立方程式を解けばQを最小化する回帰係数が求まります。偏微分と聞くと難しそうに思われるかもしれませんが、結局1つだけの変数に着目し、それ以外の変数を定数とみて微分すればいいわけです。

次回は重回帰分析の回帰係数の導出プログラムと結果の可視化をみていきます。