はじめに　ー　週末研究ノートとは？

個人的に研究的な活動をやるにあたり、オープンにしてみたら面白いかもと思い、自分が興味を持っているざっくりテーマについて、これから、ゆるい週末研究を公開していこうと思います。（有識者の方のアドバイスも、ちょっとというかかなり期待してます！笑）

どこかの権威的な学会やジャーナルなどで発表する予定はないため、万が一、私の記事を利用する際には自己責任でお願いします。そんな人はいないと思いますが、念のため。

今回のサマリ（TL; DR）

偏相関にかかわるデータの相関する成分の（$X'$ への直交射影上の）系列と無相関成分を抜き出した（直交補空間上の射影上の）系列と、重回帰モデルの予測系列と残差系列との間にそれぞれどのような関係があるか、を明確にするため、理論を整理して実機検証しました。

結論、以下のような関係になることを確認しました

今回は、理論的な側面だけではなく、実際に実機で確かめた自分を褒め称えたいと思います！笑

また、重回帰モデルの一般的でない一般解を求めることもでき、非常に有意義でした。

環境

Colab
- コードは、こちら
調査協力
- ChatGPT GPTs : Mathtermind

今回の週末研究ノート

重回帰モデルの一般的ではない一般解

$X$ を $(n, m)$ 行列、$a \in \mathbb{R}^m$, $b \in \mathbb{R}$ としたときに、最小二乗法で解くと $a, b$ は以下のように表せます

\begin{align}
a &= (X'^{\mathsf{T}} X')^{-1} X'^{\mathsf{T}} y' \\
b &= \bar{y} - \bar{X} a \\
\end{align}

ただし、

\begin{align}
y' &:= y - \bar{y} \\
X' &:= X - \bar{X} \\
\bar{y} &:= \frac{1}{n} \sum _{i=1}^n y_i \\
\bar{X} &:= [\bar{X_{*, 1}} \cdots \bar{X_{*, m}}] = ([\frac{1}{n} \sum _{i=1}^n X_{i,j}]_{j=1}^m)^{\mathsf{T}} \\
\end{align}

とします

順を追って証明をしていきましょう

重回帰モデルの切片を求める

まずは、$\hat{y} := Xa + b$ を最小二乗法で解いたとき $b$ が
$$b = \bar{y} - \bar{X} a$$
であることを示します

証明

\begin{align}
\chi &:=
	\begin{bmatrix}
	 X & \mathbb{1} \\
	\end{bmatrix} \\
\mathsf{a} &:=
	\begin{bmatrix}
	 a \\
	 b \\
	\end{bmatrix} \\
\end{align}

と置きます

ここで、 $\mathbb{1}$ は、 $\mathbb{1} := [1, ..., 1]^\mathsf{T} \in \mathbb{R}^n$ です

正規方程式は、

\begin{align}
\chi^{\mathsf{T}} \chi \mathsf{a} &= \chi^{\mathsf{T}} y \\
\end{align}

であるので、分割行列を用いて以下のように表せます

\begin{align}
	\begin{bmatrix}
	 X^{\mathsf{T}} \\
	 \mathbb{1}^{\mathsf{T}} \\
	\end{bmatrix}
	\begin{bmatrix}
	 X & \mathbb{1} \\
	\end{bmatrix}
	\begin{bmatrix}
	 a \\
	 b \\
	\end{bmatrix}
	&=
	\begin{bmatrix}
	 X^{\mathsf{T}} \\
	 \mathbb{1}^{\mathsf{T}} \\
	\end{bmatrix}
	y \\
\\ 
	\begin{bmatrix}
	 X^{\mathsf{T}} X & \sum X \\
	 (\sum X)^{\mathsf{T}} & n \\
	\end{bmatrix}
	\begin{bmatrix}
	 a \\
	 b \\
	\end{bmatrix}
	&=
	\begin{bmatrix}
	 X^{\mathsf{T}} y \\
	 \sum y \\
	\end{bmatrix}
\end{align}

2行目を展開すると

\begin{align}
(\sum X)^{\mathsf{T}} a + n b &= \sum y \\
\frac{1}{n} (\sum X)^{\mathsf{T}} a + b &= \frac{1}{n} \sum y \\
\bar{X} a + b &= \bar{y} \\
\therefore b &= \bar{y} - \bar{X} a \\
\end{align}

重回帰モデルの係数ベクトルを求める

直前の結果 $b = \bar{y} - \bar{X} a$ を用いて、係数ベクトル $a$ の一般解を求めます

ここで、

\begin{align}
\hat{y} &:= X a + b \\
\hat{y'} &:= X' a \\
\end{align}

とするとき、 $y$ の予測値 $\hat{y} := X a + b$ を最小二乗法で解いたときの解 $a$ と $y'$ の予測値 $\hat{y'} := X' a$ を最小二乗法で求めた解 $a$ が一致することを示します

式で表すと以下のようになります

\begin{align}
\min_{a ~|~ b=\bar{y} - \bar{X} a} \| y - \hat{y} \|^2 =  \min_{a} \| y' - \hat{y'} \|^2 \\
\end{align}

証明

\begin{align}
\min_{a ~|~ b=\bar{y} - \bar{X} a} \| y - \hat{y} \|^2
&=  \min_{a ~|~ b=\bar{y} - \bar{X} a} \| y - (X a + b) \|^2 \\
&=  \min_{a ~|~ b=\bar{y} - \bar{X} a} \| y - \bar{y} - (X a + b - \bar{y} ) \|^2 \\
&=  \min_{a} \| y - \bar{y} - (X a + \bar{y} - \bar{X} a - \bar{y} ) \|^2 \\
&=  \min_{a} \| y - \bar{y} - (X a - \bar{X} a) \|^2 \\
&=  \min_{a} \| y - \bar{y} - (X - \bar{X}) a \|^2 \\
&=  \min_{a} \| y' - X' a \|^2 \\
&=  \min_{a} \| y' - \hat{y'} \|^2 \\
\therefore  \min_{a ~|~ b=\bar{y} - \bar{X} a} \| y - \hat{y} \|^2 &=  \min_{a} \| y' - \hat{y'} \|^2 \\
\end{align}

つまり、$y$ の予測値 $\hat{y} := X a + b$ を最小二乗法で解いたときの解 $a$ と $y'$ の予測値 $\hat{y'} := X' a$ を最小二乗法で求めた解 $a$ が一致します

以上から、平均値を引いた $y', X'$ を用いて最小二乗法で解けばよいとわかります

正規方程式を考えると、

\begin{align}
X'^{\mathsf{T}} X' a &= X'^{\mathsf{T}} y'
\end{align}

であるので、$X'^{\mathsf{T}} X'$ が正則であるとすると、

\begin{align}
a &= (X'^{\mathsf{T}} X')^{-1} X'^{\mathsf{T}} y'
\end{align}

以上より、

\begin{align}
\min_{a, b} \| y - \hat{y} \|^2
= \min_{a, b} \| y - (X a + b) \|^2
\end{align}

を解くと、

\begin{align}
a &= (X'^{\mathsf{T}} X')^{-1} X'^{\mathsf{T}} y' \\
b &= \bar{y} - \bar{X} a \\
\end{align}

です

偏相関行列の直交射影と重回帰モデルの一般解との関係

さて、前回の記事で触れた通り、偏相関行列を求める際に導入した直交射影 $P, Q$ は、

\begin{align}
P &:= X'(X'^{\mathsf{T}}X')^{-1}X'^{\mathsf{T}} \\
Q &:= I - P
\end{align}

と表せます (注：$X'$ は、前回と今回とでまったく同じ定義です)

つまり、

\begin{align}
P y' 
&= X'(X'^{\mathsf{T}}X')^{-1}X'^{\mathsf{T}} y' \\
&= X' a \\
&= \hat{y'} \\
\end{align}

です

また、 $Q$ の定義 $Q := I - P$ より、$Q + P = I$ であることを踏まえると

\begin{align}
y' 
&= I y' \\
&= (P + Q) y' \\
&= P y' + Q y' \\
&= \hat{y'} + Q y'
\end{align}

ですので、

$y'$の予測残差を $e_{y'}$ とするとき、

\begin{align}
e_{y'} 
& = y' - \hat{y'} \\
&= Q y' \\
\end{align}

と表せます

これは、$y'$の予測残差と$y'$を $X'$ の直交補空間上に射影$Q$ した像が一致することを示しています

さらに、$e_{y} = Q y' = e_{y'}$ になることも見ていきます

\begin{align}
y'
&= P y' + Q y' \\
&= X' a + Q y' \\
&= (X - \bar{X}) a + Q y' \\
\therefore y - \bar{y} &= X a - \bar{X} a + Q y' \\
\therefore y &= X a + \bar{y} - \bar{X} a + Q y' \\
&= X a + b + Q y' \\
&= \hat{y} + Q y' \\
\therefore y - \hat{y} &= Q y'\\

\therefore e_{y} &:= y - \hat{y} = Q y' = e_{y'} \\
\end{align}

$y'$を $X'$ の直交補空間上の射影$Q$ の像は、$y$, $y'$ の $X$ の重回帰モデルの予測残差 $e_{y}, e_{y'}$ に一致する、ということがわかります

さらにさらに、

\begin{align}
\hat{y}
&= X a + b \\
&= X a + \bar{y} - \bar{X} a \\
&= (X - \bar{X}) a + \bar{y} \\
&= X' a + \bar{y} \\
&= P y' + \bar{y} \\
\end{align}

と、予測値 $\hat{y}$ を、直行射影 $P$ と平均値（ベクトル） $\bar{y}$ で表すことができます

言い換えると、

\begin{align}
y
&= \hat{y} + Q y' \\
&= P y' + \bar{y} + Q y' \\
\end{align}

と $y$ を直交射影 $P, Q$ で分解できることを意味しています

図で表すと以下のようになります

ポイント

特に、残差と直交補空間上の射影の関係を踏まえると、
系列$y_1, y_2$, と共変量$X$に対する偏相関係数は、$Q y_1, Q y_2$ の相関係数であるため、$y_1, y_2$ それぞれを共変量 $X$ による重回帰モデルの予測の残差 $e_{y_1}, e_{y_2}$ の相関係数を求めることで代用ができることを意味しています

実機検証

以上で確認してきた、等式関係を実際に数値計算でも成立することを確認していきます
全てのコードは、こちらを参照ください

import numpy as np
import pandas as pd
import pingouin as pg

直交射影 $P$ 、著効補空間上への射影 $Q$ を定義しておきます

def P(X: np.ndarray):
    if len(X.shape) == 1:
        X = X.reshape(-1, 1)
    return X.dot(np.linalg.inv(X.T.dot(X))).dot(X.T)


def Q(X: np.ndarray):
    return np.eye(X.shape[0]) - P(X)

以下のデータ系列 $X, y_1, y_2$ (コード上では、X, Y1, Y2) を使います

$X$ を共変量として、$y_1, y_2$ を生成していることにご留意ください

import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression


# 仮想データの生成（実際には実データを使用）
np.random.seed(42)

X = np.random.randn(100, 3) * 1.5 - 3   # 既知
alpha1 = [1.1, -0.7, 0.3]   # 未知
alpha2 = [-1.3, 0.2, 0.5]   # 未知
beta1 = 4   # 未知
beta2 = 1   # 未知

U1 = 3 * (np.sin(np.random.randn(100)) ** 2) - 10   # 未知の変動
U2 = 2 * (np.cos(np.random.randn(100)) ** 3) - 5    # 未知の変動
Y1 = X.dot(alpha1) + U1 + beta1     # 既知
Y2 = X.dot(alpha2) + U2 + beta2     # 既知

# Xに基づいてY1の線形回帰モデルを構築
model_Y1 = LinearRegression().fit(X, Y1)
# Y1の予測値を計算し、残差を求める
predicted_Y1 = model_Y1.predict(X)
residuals_Y1 = Y1 - predicted_Y1

# Xに基づいてY2の線形回帰モデルを構築
model_Y2 = LinearRegression().fit(X, Y2)
# Y2の予測値を計算し、残差を求める
predicted_Y2 = model_Y2.predict(X)
residuals_Y2 = Y2 - predicted_Y2

※ 作成協力 ChatGPT GPTs Mathtermind

推定したパラメータを変数として保持しておきます

a1 = model_Y1.coef_
b1 = model_Y1.intercept_
a2 = model_Y2.coef_
b2 = model_Y2.intercept_

残差間の相関係数と直交補空間上の射影の相関係数

df = pd.DataFrame({"X0": X[:, 0], "X1": X[:, 1], "X2": X[:, 2], "Y1": Y1, "Y2": Y2})
df.head(3)

index	X0	X1	X2	Y1	Y2
0	-2.254928770483151	-3.207396451756777	-2.0284671928489613	-5.2131594113529065	-2.724331871134173
1	-0.7154552153879621	-3.351230062085004	-3.351205435423771	-4.599536087596906	-4.289912336279545
2	-0.6311807767389128	-1.8488479062706369	-3.704211578902428	-5.125571848256945	-3.4014228487818787

ライブラリを使った偏相関係数

pcorr_df = pg.partial_corr(data=df, x="Y1", y="Y2", covar=["X0", "X1", "X2"])
pg_pcorr = pcorr_df["r"][0]
pg_pcorr

-0.2059147359980154

$X'$ の直交補空間上に射影した像を使った偏相関係数

# Y1, Y2 | X
pcorr = np.corrcoef(
    Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)),
    Q(X - X.mean(axis=0)).dot(Y2 - Y2.mean(axis=0)),
)[0, 1]
pcorr

-0.20591473599801607

残差間の相関係数（→　偏相関係数）

# np.corrcoef(residuals_Y1, residuals_Y2)
residual_pcorr = np.corrcoef(residuals_Y1, residuals_Y2)[0, 1]
residual_pcorr

-0.20591473599801643

この3つの偏相関係数が一致することを確認

# 3種類の方法によって算出した偏相関係数が一致することを確認
assert np.allclose(pcorr, pg_pcorr)
assert np.allclose(pcorr, residual_pcorr)

以上から、$y_1, y_2$ の $X$ に対する偏相関係数は、$y_1, y_2$ を $X$ で重回帰モデルで予測した残差 $e_{y_1}, e_{y_2}$ の相関係数と一致することがわかります

予測系列と直交射影上の像の等値性

\hat{y} = P y' + \bar{y}

重回帰モデルの予測系列

predicted_Y1[:5]

array([-5.60473505, -3.9955015 , -5.19655728, -5.60487606, -4.94716265])

$X'$ への直交射影の像系列

(P(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + Y1.mean(axis=0) )[:5]

array([-5.60473505, -3.9955015 , -5.19655728, -5.60487606, -4.94716265])

等値性の確認

assert np.allclose(
    predicted_Y1,
    P(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + Y1.mean(axis=0),
)

以上から、$y$ を $X$ による重回帰モデルで予測した系列は、$y$ を $X'$ 上の直交射影 $P$ の像に、$y$ の平均値（ベクトル）を加えたベクトル（系列）に一致することがわかります

残差と直交補空間上の像の等値性

e_{y} = Q y'

重回帰モデルの残差系列

residuals_Y1[:5]

array([ 0.39157564, -0.60403459,  0.07098544, -0.33731501, -1.51955793])

直交補空間上の像系列

Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0))[:5]

array([ 0.39157564, -0.60403459,  0.07098544, -0.33731501, -1.51955793])

等値性の確認

assert np.allclose(
    residuals_Y1,
    Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0))
)

以上から、$y$ を $X$ による重回帰モデルで予測した系列の残差 $e_{y}$は、$y$ を $X'$ に対する直補空間上の射影 $Q$ の像ベクトル（系列）に一致することがわかります

続いて、$y, y'$ の予測残差 $e_{y}, e_{y'}$ が一致することを確認します

e_{y'} = e_{y}

$y'$ を $X'$ で予測する重回帰モデルを構築

# X' = X - bar{X} に基づいて Y1' = Y1 - bar{Y1} の線形回帰モデルを構築
X_prime = X - X.mean(axis=0)
Y1_prime = Y1 - Y1.mean(axis=0)
model_Y1_prime = LinearRegression().fit(X_prime, Y1_prime)

# Y1'の 予測値を計算し、残差を求める
predicted_Y1_prime = model_Y1_prime.predict(X_prime)
residuals_Y1_prime = Y1_prime - predicted_Y1_prime

$y'$ の残差系列

residuals_Y1_prime[:5]

array([ 0.39157564, -0.60403459,  0.07098544, -0.33731501, -1.51955793])

等値性の確認

assert np.allclose(
    residuals_Y1_prime,
    residuals_Y1,
)

以上から、 $y'$ を $X$ による重回帰モデルで予測した系列の残差 $e_{y'}$は、 $y$ を $X$ による重回帰モデルで予測した系列の残差 $e_{y}$ と一致することがわかります

系列分解

以上をまとめると、以下のように系列 $y$ を分解できることがわかります

y =  P y' + \bar{y} + Q y'

$y$ 系列

Y1[:5]

array([-5.21315941, -4.59953609, -5.12557185, -5.94219107, -6.46672058])

分解系列

(P(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + Y1.mean(axis=0) + Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)))[:5]

array([-5.21315941, -4.59953609, -5.12557185, -5.94219107, -6.46672058])

等値性の確認

assert np.allclose(
    Y1,
    P(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + Y1.mean(axis=0) + Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)),
)

以上から、系列 $y$ を、直交射影 $P$ とその補空間上の射影 $Q$ を用いて分解できることがわかります

系列分解の可視化

ここからは、実際に系列を可視化して確認していきます

実際に X の成分を抽出できているか？

$X'$ へ直交射影した像系列 $P y_1'$ と $y_1$ の実際の $X$ 成分系列 $X \alpha_1$ を比較します

df = pd.DataFrame({r"$P y_1'$": P(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)), r"$X \alpha_1$":  X.dot(alpha1)})
df.plot()
plt.show()

散布図

sns.pairplot(df)
plt.show()

相関係数

np.corrcoef(P(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)), X.dot(alpha1)).round(3)

array([[1.   , 0.999],
       [0.999, 1.   ]])

散布図、相関係数から $X$ の成分を抽出できているとわかります

また、$y_1'$ の $X'$ への直交射影上の像 $P y_1'$ (変動成分) に、$X$ 成分の推定値の平均 $\bar{X a_1}$ を加えると、$y_1$に含まれる $X$ 成分 $X \alpha_1$ とほぼ一致します

df = pd.DataFrame({r"$P y_1' + \bar{X a_1}$": P(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + X.dot(a1).mean(axis=0), r"$X \alpha_1$":  X.dot(alpha1)})
df.plot()
plt.show()

つまり、$y_1$ から $X$ の成分を抽出した系列 $X \alpha_1$ は、$P y_1' + X a_1$ で近似できることを意味しています

実際に X 成分を除去できているか？

$X'$ の直交補空間上に射影した像系列 $Q y_1'$ と $y_1$ の実際の $X$ 成分系列 $X \alpha_1$ を除いた系列 $y_1 - X \alpha_1$ を比較します

df = pd.DataFrame({r"$Q y_1'$": Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)), r"$y_1-X \alpha_1$": Y1 - X.dot(alpha1)})
df.plot()
plt.show()

散布図

sns.pairplot(df)

相関係数

np.corrcoef(Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)), Y1 - X.dot(alpha1)).round(3)

array([[1.   , 0.977],
       [0.977, 1.   ]])

散布図と相関係数から、$X$ 成分を除去した系列と十分に相関が高いことから、$X$ の成分を除去しているといえます

さて、ここで、実際に $y_1$ から $X$ の成分 $X \alpha_1$ を除いた系列について考えましょう

実は、切片 $b1$ を加えるとよい近似になっていることがわかります

$y_1 = X a_1 + b_1 + Q y_1'$ であることから $y_1$ から $X$ 成分の近似 $X a_1$ を除くと、$y_1 - X a_1 = b_1 + Q y_1' ~~~(= \bar{y_1} - \bar{X_1} a_1 + Q y_1)$

assert np.allclose(Y1 - X.dot(a1), Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + b1)
assert np.allclose(Y1 - X.dot(a1), Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + Y1.mean(axis=0) - X.mean(axis=0).dot(a1))
assert np.allclose(b1, Y1.mean(axis=0) - X.mean(axis=0).dot(a1))

$a_1 = \hat{\alpha_1}$ より、$a_1 \approx \alpha_1$ であるから、 $ b_1 + Q y_1' \approx y_1 - X \alpha_1$

# using b1 which is prediction of β
pd.DataFrame({r"$Q y_1+b_1$": Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + b1, r"$y_1-X \alpha_1$": Y1 - X.dot(alpha1)}).plot()

つまり、$y_1$ から $X$ の成分を除去した系列 $y_1 - X \alpha_1$ は、$Q y_1' + b_1$ で近似できることを意味しています

おまけ

$\alpha_1$ の代わりに推定値の $a_1$ を用いると、理論通り完全に一致します

pd.DataFrame({r"$Q y_1'+b_1$": Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + b1, r"$y_1-X a_1$": Y1 - X.dot(a1)}).plot()

推定値 $b_1$ が示す意味の理解を深めるため $\beta_1$ と比較してみましょう

$y_1'$ の $X'$ の直交補空間上の像 $Q y_1'$ に、実際の定数値 $\beta_1$ と未知の変動$U_1$の平均値$\bar{U_1}$ の和 $\beta_1 + \bar{U_1}$ を加えると、$y_1$ から $X$ の成分 $X \alpha_1$ を除いた系列とほぼ一致します
- 実用上は、$\beta_1$ は未知のため利用できない点に注意が必要です

pd.DataFrame({r"$Qy_1'+\beta_1+\bar{U_1}$": Q(X - X.mean(axis=0)).dot(Y1 - Y1.mean(axis=0)) + beta1 + U1.mean(axis=0), r"$y1-X \alpha_1$": Y1 - X.dot(alpha1)}).plot()

このことから、推定値 $b_1$ は、実際の定数値 $\beta_1$ と未知変動$U_1$の平均値$\bar{U_1}$ の和 $\beta_1 + \bar{U_1}$ の推定値（近似）とみなせるでしょう

まとめ

重回帰モデルの一般的ではない一般解を求めることができた
$y_1, y_2$ の共変量 $X$ に対する偏相関係数は、それぞれに対して共変量 $X$ を説明変数とした重回帰モデルを構築した予測との残差 $e_{y_1}, e_{y_2}$ の相関係数と一致します
予測系列（ベクトル）と直交射影空間上の像系列（ベクトル）と元系列の平均値（ベクトル）との和は、一致します
- $\hat{y} = P y' + \bar{y}$
残差系列（ベクトル）と直交補空間上の像系列（ベクトル）は、一致します
- $e_{y} = Q y'$
おまけに、$y'$ の $X'$ を用いた予測残差 $e_{y'}$ は、 $y$ の $X$ を用いた予測残差 $e_{y}$ に一致します
- $e_{y'} = e_{y}$
系列（ベクトル） $y$ を $X'$ への直交射影と直交補空間上の射影を使って分解できます
- $y = P y' + \bar{y} + Q y'$
$y_1$ から $X$ の成分を抽出した系列（ベクトル） $X \alpha_1$ は、$P y_1' + \bar{X a_1}$ で近似できます
$y_1$ から $X$ の成分を除去した系列（ベクトル） $y_1 - X \alpha_1$ は、$Q y_1' + b_1$ で近似できます
- 特に、$P y' + \bar{X a} + b_1 + Q y' = P y' + \bar{\hat{y}} + Q y' \approx P y' + \bar{y} + Q y' = y$ であることに注意します
$b_1$ は、実際の定数値 $\beta_1$ と未知の変動 $U_1$ の平均値$\bar{U_1}$ の和 $\beta_1 + \bar{U_1}$ の推定値（近似）とみなせます

参考文献

理論全般
- 前回の記事公開！週末研究ノート１０　ー　偏相関行列の一般式
正規方程式
- Chainer Tutorial 7. 単回帰分析と重回帰分析