参考図書
多変量解析法入門 (ライブラリ新数学大系) 単行本 – 2001/4/10
永田 靖 (著), 棟近 雅彦 (著)

プログラミングコード、および次数を増やした時の最小二乗法はこちら

確信をもって言えることですが、
以下は間違いを含みます。
正確な記述は参考図書をご覧ください。

単回帰

誤差、あるいは残差$${\epsilon_i}$$
観測データ$${y_i}$$
観測地点$${x_i}$$は、これは例えば時間だったりするであろう。
当てはめたい関数を$${f(x_i)}$$として

観測データとの誤差、あるいは残差を

$$
\epsilon_i=y_i-f(x_i)
$$

とすると、この誤差（残差）$${\epsilon_i}$$はゼロになるのが理想であって、その時観測データ$${y_i}$$は関数$${y_i=f(x_i)}$$上にあることになる。基本的には誤差ゼロにならん場合に、誤差ゼロにならんまでも誤差の最小を目指すのが最小二乗法である。

1個の観測データの1変数(x)単回帰

$$
y=f(x)+\epsilon=a+bx+\epsilon
$$

$$
\begin{bmatrix}
y\\
\end{bmatrix}
=
\begin{bmatrix}
1 & x \\
\end{bmatrix}
\begin{bmatrix}
a \\
b \\
\end{bmatrix}
+
\begin{bmatrix}
\epsilon \\
\end{bmatrix}
$$

最小二乗法、あるいは線形回帰、単回帰は
この誤差の二乗の総和を最小化する。
すなわちn個のデータ点に対し、回帰したい関数（この場合直線）との誤差を求め、その総和をとったうえでその式を係数（上記の場合$${a,b}$$）で微分=0とし、連立方程式を解く。

誤差は

$$
\epsilon = y - (a + bx)
$$

２乗して

$$
\epsilon^2 = \lbrace y - (a + bx) \rbrace^2\\
=y^2-ay-bxy-ay+a^2+abx-bxy+abx+b^2x^2\\
=y^2-2ay-2bxy+2abx+a^2+b^2x^2
$$

n=1なので、総和は勘弁しておくとして
係数$${a,b}$$に関して微分して=0をとると

$$
\frac{\partial \epsilon^2}{\partial a}=2a+2bx-2y=0\\
\frac{\partial \epsilon^2}{\partial b}=2ax+2bx^2-2xy=0
$$

連立方程式は

$$
a+bx=y\\
ax+bx^2=xy
$$

行列にすると

$$
\begin{bmatrix}
1 & x\\
x & x^2
\end{bmatrix}
\begin{bmatrix}
a\\
b
\end{bmatrix}
=
\begin{bmatrix}
y\\
xy
\end{bmatrix}
$$

この時

$${\begin{bmatrix}1 & x\\x & x^2\end{bmatrix}}$$の1行目をx倍して2行目から減算すると$${\begin{bmatrix}1 & x\\0 & 0\end{bmatrix}}$$となる。
ガウスの消去法を用いて対角に0が現れる場合、その連立方程式は解けない。

今回のケースは解多すぎてとけないのパターン。
すなわち１つの点を通過する直線などいくらでもあるのパターン。

2個の観測データの1変数(x)単回帰

$$
\begin{bmatrix}
y_1\\
y_2\\
\end{bmatrix}
=
\begin{bmatrix}
1 & x _1\\
1 & x _2\\
\end{bmatrix}
\begin{bmatrix}
a \\
b \\
\end{bmatrix}
+
\begin{bmatrix}
\epsilon \\
\end{bmatrix}
$$

誤差の２乗は

$$
\epsilon^2 = \lbrace y-(a+bx) \rbrace^2
$$

n個の観測データの誤差２乗の総和を関数$${S_e}$$とすると。(S:sum, e:error)

$$
S_e=\sum\limits_{i=1}^n \epsilon_i^2=\sum\limits_{i=1}^n\lbrace y_i-(a+bx_i) \rbrace^2
$$

今、n=2。
係数$${\lbrace a, b \rbrace}$$で微分して０ならば最小値（放物線の底、おわんの底）
２乗の効果で、係数に関して、誤差の２乗の総和関数は負値をとらない。変な形をしない。係数に関して微分してゼロならば最小値である。
ただし最小値を満たす点が多すぎることはあり得る。逆に言うと微分してゼロの点は最大値ではないし変曲点でもないし鞍点でもないと思われる。多分。

$$
\frac{\partial S_e}{\partial a} = -2\sum\limits_{i=1}^n (y_i-(a+bx_i)) = 0\\
\frac{\partial S_e}{\partial b} = -2\sum\limits_{i=1}^n x_i(y_i-(a+bx_i)) = 0
$$

$$
\frac{\partial S_e}{\partial a}=2(a+bx_1-y_1)+2(a+bx_2-y_2)=0\\
\frac{\partial S_e}{\partial b}=2(ax_1+bx_1^2-x_1y_1)+2(ax_2+bx_2^2-x_2y_2)=0
$$

連立方程式は

$$
2a+b(x_1+x_2)=(y_1+y_2)\\
a(x_1+x_2)+b(x_1^2+x_2^2)=(x_1y_1+x_2y_2)
$$

行列にすると

$$
\begin{bmatrix}
n & \sum\limits_{i=1}^n x_i\\
\sum\limits_{i=1}^n x_i & \sum\limits_{i=1}^n x_i^2
\end{bmatrix}
\begin{bmatrix}
a\\
b
\end{bmatrix}
=
\begin{bmatrix}
\sum\limits_{i=1}^n y_i\\
\sum\limits_{i=1}^n x_iy_i
\end{bmatrix}
$$

今、n=2。
$${(x_1,y_1)=(x_2,y_2)}$$とかでない限り解ける。
ただし結果はデータ点を通るただの直線。

n個の観測データの1変数(x)単回帰

$$
y_i=f(x_i)+\epsilon_i=a+bx_i+\epsilon_i \quad i=1,2,…,n
$$

$$
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_n\\
\end{bmatrix}
=
\begin{bmatrix}
1 & x_{1} \\
1 & x_{2} \\
\vdots & \vdots\\
1 & x_{n} \\
\end{bmatrix}
\begin{bmatrix}
a \\
b \\
\end{bmatrix}
+
\begin{bmatrix}
\epsilon_1 \\
\epsilon_2 \\
\vdots\\
\epsilon_n\\
\end{bmatrix}
$$

誤差を２乗して総和

$$
S_e=\sum\limits_{i=1}^n \epsilon_i^2=\sum\limits_{i=1}^n\lbrace y_i-(a+bx_i) \rbrace^2
$$

連立方程式は

$$
na+b\sum\limits_{i=1}^nx_i=\sum\limits_{i=1}^n y_i\\
a\sum\limits_{i=1}^n x_i+b \sum\limits_{i=1}^n x_i^2=\sum\limits_{i=1}^n x_iy_i
$$

nは$${\sum\limits_{i=1}^n1}$$でもよい。

行列にすると

$$
\begin{bmatrix}
n & \sum\limits_{i=1}^n x_i\\
\sum\limits_{i=1}^n x_i & \sum\limits_{i=1}^n x_i^2
\end{bmatrix}
\begin{bmatrix}
a\\
b
\end{bmatrix}
=
\begin{bmatrix}
\sum\limits_{i=1}^n y_i\\
\sum\limits_{i=1}^n x_iy_i
\end{bmatrix}
$$

正規方程式

$$
\mathbf y=
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_n\\
\end{bmatrix},
\mathbf X =
\begin{bmatrix}
1 & x_{1} \\
1 & x_{2} \\
\vdots & \vdots\\
1 & x_{n} \\
\end{bmatrix},
\bm \beta=
\begin{bmatrix}
\beta_0 \\
\beta_1 \\
\end{bmatrix},
\bm \epsilon=
\begin{bmatrix}
\epsilon_1 \\
\epsilon_2 \\
\vdots\\
\epsilon_n\\
\end{bmatrix}
$$

とすると、

$$
\bm \epsilon = \mathbf y - \bm X\bm \beta
$$

誤差の２乗の総和は$${\bm \epsilon}$$同士の内積で表現できると考えてよく、行列なら

$$
\bm \epsilon^T\bm \epsilon
$$

が誤差の２乗の総和を表す。

$$
\bm \epsilon^T\bm \epsilon= (\mathbf y - \bm X\bm \beta)^T(\mathbf y - \bm X\bm \beta)\\
=(\mathbf y^T-\bm \beta^T\bm X^T)(\mathbf y - \bm X\bm \beta)\\
=\mathbf y^T\mathbf y-\mathbf y^T\bm X\bm \beta-\bm \beta^T\bm X^T\mathbf y+\bm \beta^T\bm X^T\bm X\bm \beta
$$

ここで$${\mathbf y^T\bm X\bm \beta}$$は
(1*N)行列と
(N*2)行列と
(2*1)行列の積であるからして、結果出力はスカラー
スカラーなので転置しても変わらず
転置の性質$${(AB)^T=B^TA^T}$$を用いて

$${\mathbf y^T\bm X\bm \beta=(\mathbf y^T\bm X\bm \beta)^T=(\bm X\bm \beta)^T\mathbf y=\bm \beta^T\bm X^T\mathbf y}$$

よって誤差の２乗の総和は

$$
S_e=\bm \epsilon^T\bm \epsilon=\mathbf y^T\mathbf y-2\bm \beta^T\bm X^T\mathbf y+\bm \beta^T\bm X^T\bm X\bm \beta
$$

これを微分して=0とすると

$$
\frac{\partial S_e}{\partial \bm \beta}=-2\bm X^T\mathbf y+2\bm X^T\bm X\bm \beta = 0
$$

ここで$${\mathbf y^T\mathbf y}$$は$${\partial \bm \beta}$$に応じないからそのまま消える。
$${-2\bm \beta^T\bm X^T\mathbf y}$$を$${-2(\bm \beta^T)(\bm X^T\mathbf y)}$$とみると、これは(1*2)と(2*N)(N*1)の積であり、これは内積である。$${\bm \beta^T\bm X^T\bm X\bm \beta}$$は二次形式。

内積のベクトル微分は

$$
\frac{\partial \bm x^T\bm a}{\partial \bm x}=\bm a
$$

二次形式のベクトル微分は

$$
\frac{\partial \bm x^TA\bm x}{\partial \bm x}=2A\bm x
$$

よって

$$
\frac{\partial S_e}{\partial \bm \beta}=-2\bm X^T\mathbf y+2\bm X^T\bm X\bm \beta = 0
$$

から

$$
\bm \beta =(\bm X^T\bm X)^{-1}\bm X^T\mathbf y
$$

この形は重回帰でも同じ。

正規方程式にはQR分解がハマり、その場合

$$
\bm X=QR\\
R\bm\beta=\bm X^T\mathbf y
$$

を解く。

正規方程式が解ける場合

変数側が非線形であっても、
係数（ウェイト）側が線形なら
正規方程式が立つ。
あとは連立方程式が解ける解けない、
すなわちX行列が逆行列を持つかどうかの問題になる。

変数側が関数であっても成り立つ。
総じて線形回帰モデルと呼ばれる。

正規方程式が解けない場合

座標が重複している場合。
（X行列の列ベクトルの重複）
座標の定数倍が含まれている場合
（X行列の列ベクトルに、他の列ベクトル×定数がおる場合）

重回帰

1個の観測データの2変数(x,y)重回帰

$$
z=f(x, y)+\epsilon=a+bx+cy+\epsilon
$$

$$
\begin{bmatrix}
z
\end{bmatrix}
=
\begin{bmatrix}
1 & x & y \\
\end{bmatrix}
\begin{bmatrix}
a \\
b \\
c \\
\end{bmatrix}
+
\begin{bmatrix}
\epsilon
\end{bmatrix}
$$

誤差の２乗は

$$
\epsilon^2=\lbrace z-(a+bx+cy) \rbrace^2
$$

$$
\epsilon^2
=(z-a-bx-cy)(z-a-bx-cy)\\
=z(z-a-bx-cy)\\
-a(z-a-bx-cy)\\
-bx(z-a-bx-cy)\\
-cy(z-a-bx-cy)\\
=z^2-az-bxz-cyz\\
-az+a^2+abx+acy\\
-bxz+abx+b^2x^2+bcxy\\
-cyz+acy+bcxy+c^2y^2\\
=z^2+a^2+b^2x^2+c^2y^2-2az-2bxz-2cyz+2abx+2acy+2bcxy
$$

係数が3つなので、
微分=0からなる連立方程式は3つ。

$$
\frac{\partial \epsilon^2}{\partial a}=a+bx+cy-z=0\\
\frac{\partial \epsilon^2}{\partial b}=ax+bx^2+cxy-xz=0\\
\frac{\partial \epsilon^2}{\partial c}=ay+bxy+cy^2-yz=0
$$

行列にすると

$$
\begin{bmatrix}
1 & x & y\\
x & x^2 & xy\\
y & xy & y^2
\end{bmatrix}
\begin{bmatrix}
a\\
b\\
c
\end{bmatrix}
=
\begin{bmatrix}
z\\
xz\\
yz
\end{bmatrix}
$$

$$
\begin{bmatrix}
1\\
x\\
y
\end{bmatrix}
\begin{bmatrix}
1 & x & y
\end{bmatrix}
\begin{bmatrix}
a\\
b\\
c
\end{bmatrix}
=
\begin{bmatrix}
1\\
x\\
y
\end{bmatrix}
z
$$

p変数に一般化すると

$$
\begin{bmatrix}
1\\
x_1\\
x_2\\
\vdots\\
x_p
\end{bmatrix}
\begin{bmatrix}
1 & x_1 & x_2 & \cdots & x_p
\end{bmatrix}
\begin{bmatrix}
\beta_0\\
\beta_1\\
\beta_2\\
\vdots\\
\beta_p
\end{bmatrix}
=
\begin{bmatrix}
1\\
x_1\\
x_2\\
\vdots\\
x_p
\end{bmatrix}
z
$$

これは正規方程式

$$
\bm \beta =(\bm X^T\bm X)^{-1}\bm X^T\mathbf y
$$

と同じ。
ただし、xはデータの個数分行方向に伸びるため、
xは標準が横ベクトルで転置したら縦ベクトル。

n個の観測データの2変数(x,y)重回帰

$$
z_i=f(x_i, y_i)+\epsilon_i=a+bx_i+cy_i+\epsilon_i \quad i=1,2,…,n
$$

$$
\begin{bmatrix}
z_1 \\
z_2 \\
\vdots\\
z_n\\
\end{bmatrix}
=
\begin{bmatrix}
1 & x_{1} & y_{1} \\
1 & x_{2} & y_{2} \\
\vdots& \vdots& \vdots\\
1 & x_{n} & y_{n} \\
\end{bmatrix}
\begin{bmatrix}
a \\
b \\
c \\
\end{bmatrix}
+
\begin{bmatrix}
\epsilon_1 \\
\epsilon_2 \\
\vdots\\
\epsilon_n\\
\end{bmatrix}
$$

誤差を２乗して総和

$$
S_e=\sum\limits_{i=1}^n \epsilon_i^2=\sum\limits_{i=1}^n\lbrace z_i-(a+bx_i+cy_i) \rbrace^2
$$

微分=0からなる連立方程式は

$$
\frac{\partial S_e}{\partial a}=(a+bx_1+cy_1)+(a+bx_2+cy_2)-z_1-z_2=0\\
\frac{\partial S_e}{\partial b}=(ax_1+bx_1^2+cx_1y_1)+(ax_2+bx_2^2+cx_2y_2)-x_1z_1-x_2z_2=0\\
\frac{\partial S_e}{\partial c}=(ay_1+bx_1y_1+cy_1^2)+(ay_2+bx_2y_2+cy_2^2)-y_1z_1-y_2z_2=0
$$

行列にすると

$$
\begin{bmatrix}
n & \sum\limits_{i=1}^n x_i & \sum\limits_{i=1}^n y_i\\
\sum\limits_{i=1}^n x_i & \sum\limits_{i=1}^n x_i^2 & \sum\limits_{i=1}^n x_iy_i\\
\sum\limits_{i=1}^n y_i & \sum\limits_{i=1}^n x_iy_i & \sum\limits_{i=1}^n y_i^2
\end{bmatrix}
\begin{bmatrix}
a\\
b\\
c
\end{bmatrix}
=
\begin{bmatrix}
\sum\limits_{i=1}^n z_i\\
\sum\limits_{i=1}^n x_iz_i\\
\sum\limits_{i=1}^n y_iz_i
\end{bmatrix}
$$

nは$${\sum\limits_{i=1}^n1}$$でもよい。

n個の観測データの1変数単回帰が

$$
\begin{bmatrix}
n & \sum\limits_{i=1}^n x_i\\
\sum\limits_{i=1}^n x_i & \sum\limits_{i=1}^n x_i^2
\end{bmatrix}
\begin{bmatrix}
a\\
b
\end{bmatrix}
=
\begin{bmatrix}
\sum\limits_{i=1}^n y_i\\
\sum\limits_{i=1}^n x_iy_i
\end{bmatrix}
$$

だったので、
見ただけでn個データのp変数重回帰に一般化できれば
苦労はないけれども、
それができる人はこのページ見とらんと考えて、
ここでは正規方程式にのっとって

$$
\bm X^T\bm X\bm \beta=\bm X^T\mathbf y
$$

の$${\bm X^T\bm X}$$部分は

$$
\begin{bmatrix}
1 & 1 & 1\\
x_1 & x_2 & x_3\\
y_1 & y_2 & y_3
\end{bmatrix}
\begin{bmatrix}
1 & x_1 & y_1\\
1 & x_2 & y_2\\
1 & x_3 & y_3
\end{bmatrix}\\
=
\begin{bmatrix}
1+1+1 & x_1+x_2+x_3 & y_1+y_2+y_3\\
x_1+x_2+x_3 & x_1^2+x_2^2+x_3^2 & x_1y_1+x_2y_2+x_3y_3\\
y_1+y_2+y_3 & x_1y_1+x_2y_2+x_3y_3 & y_1^2+y_2^2+y_3^2
\end{bmatrix}\\
$$

変数は2つのまま、n個のデータに一般化すると

$$
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1\\
x_1 & x_2 & x_3 & \cdots & x_n\\
y_1 & y_2 & y_3 & \cdots & y_n
\end{bmatrix}
\begin{bmatrix}
1 & x_1 & y_1\\
1 & x_2 & y_2\\
1 & x_3 & y_3\\
\vdots &\vdots &\vdots\\
1 & x_n & y_n
\end{bmatrix}
=
\begin{bmatrix}
n & \sum\limits_{i=1}^n x_i & \sum\limits_{i=1}^n y_i\\
\sum\limits_{i=1}^n x_i & \sum\limits_{i=1}^n x_i^2 & \sum\limits_{i=1}^n x_iy_i\\
\sum\limits_{i=1}^n y_i & \sum\limits_{i=1}^n x_iy_i & \sum\limits_{i=1}^n y_i^2
\end{bmatrix}
$$

$$
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1\\
x_1 & x_2 & x_3 & \cdots & x_n\\
y_1 & y_2 & y_3 & \cdots & y_n
\end{bmatrix}
\begin{bmatrix}
1 & x_1 & y_1\\
1 & x_2 & y_2\\
1 & x_3 & y_3\\
\vdots &\vdots &\vdots\\
1 & x_n & y_n
\end{bmatrix}
\begin{bmatrix}
a\\
b\\
c
\end{bmatrix}
=
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1\\
x_1 & x_2 & x_3 & \cdots & x_n\\
y_1 & y_2 & y_3 & \cdots & y_n
\end{bmatrix}
\begin{bmatrix}
z_1\\
z_2\\
z_3\\
\vdots\\
z_n
\end{bmatrix}
$$

1個の観測データのp変数重回帰

$$
Y=f(X_1,X_2,…,X_p)=\beta_0+\beta_1X_1+\beta_2X_2+…\beta_pX_p+\epsilon
$$

$$
\begin{bmatrix}
y_1 \\
\end{bmatrix}
=
\begin{bmatrix}
1 & x_{11} & x_{12} & \cdots & x_{1p}\\
\end{bmatrix}
\begin{bmatrix}
\beta_0 \\
\beta_1 \\
\vdots\\
\beta_p\\
\end{bmatrix}
+
\begin{bmatrix}
\epsilon_1 \\
\end{bmatrix}
$$

連立方程式は前述の通り

$$
\begin{bmatrix}
1\\
x_1\\
x_2\\
\vdots\\
x_p
\end{bmatrix}
\begin{bmatrix}
1 & x_1 & x_2 & \cdots & x_p
\end{bmatrix}
\begin{bmatrix}
\beta_0\\
\beta_1\\
\beta_2\\
\vdots\\
\beta_p
\end{bmatrix}
=
\begin{bmatrix}
1\\
x_1\\
x_2\\
\vdots\\
x_p
\end{bmatrix}
z
$$

この式は前述の通り正規方程式

$$
\bm \beta =(\bm X^T\bm X)^{-1}\bm X^T\mathbf y
$$

n個の観測データのp変数重回帰

$$
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots\\
y_n\\
\end{bmatrix}
=
\begin{bmatrix}
1        & x_{11} & x_{12} & \cdots & x_{1p}\\
1        & x_{21} & x_{22} & \cdots & x_{2p}\\
\vdots& \vdots& \vdots& \vdots & \vdots\\
1        & x_{n1} & x_{n2} & \cdots & x_{np} \\
\end{bmatrix}
\begin{bmatrix}
\beta_0 \\
\beta_1 \\
\vdots\\
\beta_p\\
\end{bmatrix}
+
\begin{bmatrix}
\epsilon_1 \\
\epsilon_2 \\
\vdots\\
\epsilon_n\\
\end{bmatrix}
$$

連立方程式というか正規方程式は単回帰の時と変わらず

$$
\bm \beta =(\bm X^T\bm X)^{-1}\bm X^T\mathbf y
$$

$$
\begin{bmatrix}
1 & 1 & \cdots & 1\\
x_{11} & x_{21} & \cdots & x_{n1}\\
x_{12} & x_{22} & \cdots & x_{n2}\\
\vdots & \vdots & \ddots & \vdots\\
x_{1p} &  x_{2p} &  \cdots & x_{np}
\end{bmatrix}
\begin{bmatrix}
1 & x_{11} & x_{12} & \cdots & x_{1p}\\
1 & x_{21} & x_{22} & \cdots & x_{2p}\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
1 & x_{n1} & x_{n2} & \cdots & x_{np}\\
\end{bmatrix}
\begin{bmatrix}
\beta_0\\
\beta_1\\
\beta_2\\
\vdots\\
\beta_p
\end{bmatrix}
=
\begin{bmatrix}
1 & 1 & \cdots & 1\\
x_{11} & x_{21} & \cdots & x_{n1}\\
x_{12} & x_{22} & \cdots & x_{n2}\\
\vdots & \vdots & \ddots & \vdots\\
x_{1p} &  x_{2p} &  \cdots & x_{np}
\end{bmatrix}
\begin{bmatrix}
y_1\\
y_2\\
y_3\\
\vdots\\
y_n
\end{bmatrix}
$$