回帰直線と射影
説明が明確で、分かりやすいストラング：線形代数イントロダクションは、線形代数の入門書として多くの人に支持されています。しかし、ハードルもあります。理解に多くに時間を費やしてしまうこともあります。その一つが、4.3の最小二乗法です。最小二乗法そのものは難しいものではありません。回帰直線を求めるのに最初に習うのが最小二乗法です。誤差を最小にするように回帰係数を決めましょうということなので、なんら難しさはありません。しかし、これを射影と合わせて理解しようとすると、それはやさしいことではありません。それは図4.6を見ればわかります。
左の図は回帰直線を幾何的に説明した図で、右側は線形代数的に表現したものです。ここでは時間の経過$t$を0,1,2としています。また、その観測値$\mathbf{b}$を6,0,0としています。これを
$\mathbf{b}=C+Dt$
という式で当てはめます。
そうすると
$C+D\cdot 0=6$
$C+D\cdot 1=0$
$C+D\cdot 2=0$
となります。しかし、この式は一意の解をもちません。なぜなら、この連立方程式は、回帰係数が2つなのに式は3つあり、式の数が多すぎます。これを
A=\left[\begin{array}{c,c}1&0\\1&1\\ 1&2 \end{array}\right]
\mathbf{x}=\left[\begin{array}{l} C\\ D\end{array} \right]
\mathbf{b}=\left[\begin{array}{l} 6\\ 0 \\ 0\end{array} \right]
と書きなおします。$A$は$m \times n$行列でこの場合は$m=3,n=2$です。$A\mathbf{x}=\mathbf{b}$は一意の解を持ちません。
そこで、書き直し
$\mathbf{y}=\hat{C}+\hat{D}t+e$
とします。$\hat{C}$と$\hat{D}$は推定値です。$\mathbf{b}$の予測値を求めるのですが、それでも説明できない部分があるのでそれを$\mathbf{e}$で表しています。これを誤差(残差)と呼びます。
$\mathbf{x}$を$\hat{\mathbf{x}}$として、条件を加えて解きます。$\mathbf{p}=A\hat{\mathbf{x}}$とします。$\mathbf{p}$は射影と呼ばれるものです。図4.6の左側の図では横軸が時間で縦軸は観測値、または予測値、射影です。2次元の平面上に描かれています。$\mathbf{b}$と$\mathbf{p}$が一致していないのが分かります。その差$\mathbf{e}$を誤差(残差)といいます。幾何学的には、この誤差(残差)の二乗和を最小にするように$\hat{\mathbf{x}}$を求めます。
右の図も同じことを説明していますが、表現の仕方が異なります。$a_1,a_2$はベクトルです。$a_1$は$C$を推定するために設けられたベクトルです。$a_2$は時間です。この2つのベクトルは行列$A$の中の縦ベクトルです。この2つのベクトルの線形結合は列空間を作ります。$\mathbf{p}$はこの列空間上の特定の点となります。$\mathbf{e}=\mathbf{b}-\mathbf{p}$ですが、これを誤差(残差)ベクトルといいます。この誤差が最小になるように$\hat{\mathbf{x}}$を決めます。つまり、この誤差(残差)ベクトルと列空間が直交するように$\hat{\mathbf{x}}$を求めます。
幾何的な表現では2次元の平面の上にいくつでも時間を発展させていくことができます。代数的な表現でも同じように時間を発展させていくことができますが、可視化しようとするとそこには限界があり、時間にすると0,1,2までです。これは縦ベクトルの要素の数は次元の数になるので、可視化では3次元が限界になります。このことをもう少し深く理解するために列空間を考えてみましょう。
列空間の理解
まず、$C,D$を任意にしてみましょう。そうするとそれは列ベクトルが張る列空間になります。実際には$C,D$は無限になるのですが、ここではPythonを用いて乱数を発生され無限の$C,D$を疑似的に再現します。乱数で生成した$C,D$を$(C^*,D^*)=\mathbf{x}^*$として
$\mathbf{q}=A\mathbf{x}^*$
とします。
時間が0のみの場合は$A=[1, 0]$になり、$C,D$を任意の数とすると、直線上に広がります。これは自明なのでPythonでの描画は省略します。
時間が0と1の場合は、
A=\left[\begin{array}{c,c}1&0\\1&1 \end{array}\right]
となり、やはり$C,D$を任意の数とすると、これは平面上に広がります。Pythonを使って、$C,D$を-10から10の間の任意の数として平面を描いてみましょう。右肩上がりの平面に見えますが、これは$A$の列空間ですのでz軸の値は$z=0$です。列空間は３次元上にあり、要因の数が２つなので、３次元空間上のAの線形結合の作る平面となります。
\mathbf{q}=\left[\begin{array}{c,c}1&0\\1&1 \end{array}\right]\left[\begin{array}{c}C^*\\D^* \end{array}\right]=\left[\begin{array}{c}q_1\\q_2 \end{array}\right]
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (10, 3))
ax = fig.add_subplot(131, projection='3d')
ax.set_xlabel("q1", size = 10, color = "r")
ax.set_ylabel("q2", size = 10, color = "r")
ax.set_xticks([-10.0, -5, 0.0, 5, 10.0])
ax.set_yticks([-10.0, -5, 0.0, 5, 10.0])
xx=np.random.rand(10000,2)*20-10
A=np.array([[1,0],
[1,1]])
q=[A@xxx for xxx in xx]
q1=(np.transpose(q))[0]
q2=(np.transpose(q))[1]
ax.scatter(q1, q2, s = 1, c = "blue")
ax2 = fig.add_subplot(132, projection='3d')
ax2.set_zlabel("z", size = 10, color = "r")
ax2.set_ylabel("q2", size = 10, color = "r")
ax2.scatter(q1, q2, s = 1, c = "blue")
ax2.view_init(elev=0, azim=0)
ax3 = fig.add_subplot(133, projection='3d')
ax3.set_xlabel("q1", size = 10, color = "r")
ax3.set_ylabel("q2", size = 10, color = "r")
ax3.scatter(q1, q2, s = 1, c = "blue")
ax3.view_init(elev=90, azim=0)
plt.tight_layout()
plt.show()
つぎに同じことを時間が0,1,2の場合について行ってみましょう。列空間は3次元空間上の平面を作っているのが分かります。この場合には少し斜めになっているのが分かります。これは$A$の線形結合の集合を示しています。
\mathbf{q}=\left[\begin{array}{c,c}1&0\\1&1\\1&2 \end{array}\right]\left[\begin{array}{c}C^*\\D^* \end{array}\right]=\left[\begin{array}{c}q_1\\q_2\\q_3 \end{array}\right]
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (3, 3))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("q1", size = 10, color = "r")
ax.set_ylabel("q2", size = 10, color = "r")
ax.set_zlabel("q3", size = 10, color = "r")
ax.set_xticks([-10.0, -5, 0.0, 5, 10.0])
ax.set_yticks([-10.0, -5, 0.0, 5, 10.0])
xx=np.random.rand(10000,2)*20-10
A=np.array([[1,0],
[1,1],
[1,2]])
q=[A@xxx for xxx in xx]
q1=(np.transpose(q))[0]
q2=(np.transpose(q))[1]
q3=(np.transpose(q))[2]
ax.scatter(q1, q2, q3, s = 1, c = "blue")
plt.show()
つぎに時間を0,1,2,3にしてみましょう。
\mathbf{q}=\left[\begin{array}{c,c}1&0\\1&1\\1&2 \\ 1&3 \end{array}\right]\left[\begin{array}{c}C^*\\D^* \end{array}\right]=\left[\begin{array}{c}q_1\\q_2\\q_3\\q_4 \end{array}\right]
列空間は4次元空間上にあるのですが、それは表現できません。下の図は3次元の図として描いていますが、これは正確ではありません。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (3, 3))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("q1", size = 10, color = "r")
ax.set_ylabel("q2", size = 10, color = "r")
ax.set_zlabel("q3", size = 10, color = "r")
ax.set_xticks([-10.0, -5, 0.0, 5, 10.0])
ax.set_yticks([-10.0, -5, 0.0, 5, 10.0])
xx=np.random.rand(10000,2)*20-10
A=np.array([[1,0],
[1,1],
[1,2],
[1,3]])
q=[A@xxx for xxx in xx]
q1=(np.transpose(q))[0]
q2=(np.transpose(q))[1]
q3=(np.transpose(q))[2]
q4=(np.transpose(q))[3]
ax.scatter(q1, q2,q3, s = 1, c = "blue")
plt.show()
そこで4次元目をおなじｚ軸上に描いてみました。赤の面が4次元目のw軸ですが面が青い3次元の面と一致しないのが分かります。もし4次元の可視化が可能なら、列空間は平面になることが期待されます。
fig = plt.figure(figsize = (3, 3))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("q1", size = 10, color = "r")
ax.set_ylabel("q2", size = 10, color = "r")
ax.set_zlabel("q4 or q3", size = 10, color = "r")
ax.set_xticks([-10.0, -5, 0.0, 5, 10.0])
ax.set_yticks([-10.0, -5, 0.0, 5, 10.0])
ax.scatter(q1, q2,q3, s = 1, c = "blue")
ax.scatter(q1, q2,q4, s = 1, c = "red")
plt.show()
いよいよ、条件を付けて$\hat{\mathbf{x}}$を求めてみます。
A=\left[\begin{array}{c,c}1&0\\1&1\\ 1&2 \end{array}\right]
\mathbf{x}=\left[\begin{array}{l} C\\ D\end{array} \right]
b=\left[\begin{array}{l} 6\\ 0 \\ 0\end{array} \right]
$A\mathbf{x}=\mathbf{b}$を解きますが、まず$A\hat{\mathbf{x}}=\mathbf{p}$と置きます。$A$の部分空間である$\mathbf{p}$と$\mathbf{e}=\mathbf{b}-A\hat{\mathbf{x}}$のベクトルが直交するようにします。
$A^T(\mathbf{b}-A\hat{\mathbf{x}})=0$ または $A^T\mathbf{b}=A^TA\hat{\mathbf{x}}$
これを$\hat{\mathbf{x}}$について解くと
$\hat{\mathbf{x}}=(A^TA)^{-1}A^T\mathbf{b}$
が得られます。よって
\hat{\mathbf{x}}=\left[\begin{array}{c}5\\-3 \end{array}\right]
が得られます。したがって、$\mathbf{p}=A\hat{\mathbf{x}}$より
\mathbf{p}=\left[\begin{array}{c}5\\2\\ -1 \end{array}\right]
となります。つまり誤差を最小にする$\hat{\mathbf{x}}$は一意に定まり、その誤差ベクトルは
\mathbf{e}=\left[\begin{array}{c}1\\-2\\ 1 \end{array}\right]
となります。
b=np.array([[6],[0],[0]])
A=np.array([[1,0],
[1,1],
[1,2]])
print(A.T@A)
print(A.T@b)
x_hat=np.linalg.inv((A.T@A))@(A.T@b)
p=A@x_hat
x_hat,p,b-p
これを赤い点で列空間上に描いてみましょう。
fig = plt.figure(figsize = (3, 3))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("p1", size = 10, color = "r")
ax.set_ylabel("p2", size = 10, color = "r")
ax.set_zlabel("p3", size = 10, color = "r")
ax.set_xticks([-10.0, -5, 0.0, 5, 10.0])
xx=np.random.rand(1000,2)*20-10
A=np.array([[1,0],
[1,1],
[1,2]])
p=[A@xxx for xxx in xx]
p1=(np.transpose(p))[0]
p2=(np.transpose(p))[1]
p3=(np.transpose(p))[2]
ax.scatter(p1, p2, p3, s = 1, c = "blue")
ax.scatter(5,2,-1, s = 10, c = "red")
plt.show()
この赤い点は$\mathbf{b}$と、$\mathbf{e}$が列空間と直交するという条件を与えたときの列空間の1点です。直交という言葉から列空間の赤い点の上に$\mathbf{b}$がありそうですが、$\mathbf{b}$はベクトルで図4.6の右側の図にあるように空間的には列空間とは別の場所にあると考えた方がよさそうです。図4.7では$\mathbf{b}$は列空間内でもなく、$A^T$の零空間でもないところにあります。$\mathbf{e}$は$A^T$の零空間のへりにあります。
つぎにこれを横軸を時間、縦軸を予測値・観測値とした図を描いてみると直線であるのが分かります。これは射影軸です。
fig = plt.figure(figsize = (3, 3))
A0=np.linspace(-1,3,100)
O=np.ones(100)
A=np.stack([O,A0],1)
plt.scatter(A0,A@x_hat,s=1,c='blue')
plt.scatter(0,6,s=10,c='red')
plt.scatter(1,0,s=10,c='red')
plt.scatter(2,0,s=10,c='red')
plt.scatter(0,5,s=20,c='green')
plt.scatter(1,2,s=20,c='green')
plt.scatter(2,-1,s=20,c='green')
fig = plt.figure(figsize = (3, 3))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("p1", size = 10, color = "r")
ax.set_ylabel("p2", size = 10, color = "r")
ax.set_zlabel("p3", size = 10, color = "r")
ax.set_xticks([-10.0, -5, 0.0, 5, 10.0])
xx=np.array([[5,-3]])
A=np.array([[1,0],
[1,1],
[1,2]])
p=[A@xxx for xxx in xx]
p1=(np.transpose(p))[0]
p2=(np.transpose(p))[1]
p3=(np.transpose(p))[2]
ax.scatter(p1, p2, p3, s = 100, c = "blue")
#ax.scatter(5,2,-1, s = 10, c = "red")
plt.show()
p.232の図4.6の右と左の図は賢明に読み取れば、右の図と左の図を矢印で結び付けて左図の回帰直線を右図のどこかに見出せて、また、左図のどこかに直交を見出せるのではないかと思いましたが、それは考えるだけ無駄でした。
おまけ
観測値$\mathbf{b}$を乱数により生成し、それを乱数により構成した$A$から最小二乗法により予測値を計算し、この予測値$\mathbf{p}$と$\mathbf{b}$の線形性をグラフによりチェックします。$\hat{x}$はモデルパラメータの推定値、$\mathbf{b}$は残差です。
m=3
b=np.random.rand(1,m)
b=b.T
print('b',b)
A=np.random.rand(m,2)
print('A',A)
x_hat=np.linalg.inv((A.T@A))@(A.T@b)
p=A@x_hat
print('hat_x',x_hat,'p',p,'e',b-p)
fig = plt.figure(figsize = (3, 3))
plt.scatter(b,p)
plt.xlabel("b", size = 10, color = "r")
plt.ylabel("p", size = 10, color = "b")
plt.show()
もともと$A$の線形結合でない$\mathbf{b}$は射影軸とは直線の関係になりません。
つぎにモデルパラメータ$\mathbf{x}$を乱数で生成し、その説明変数$A$を乱数で生成し、線形結合$A\mathbf{x}$により$\mathbf{p}$を生成します。それにノイズを載せます。ノイズは小さなノイズと大きなノイズを作ります。ノイズが小さければ線形結合の関係が維持され、射影軸は直線になります。ノイズが大きければ線形性は壊れ、射影軸は直線になりません。
m=3
n=2
x=np.random.rand(1,n)
x=x.T
A=np.random.rand(m,2)
p=A@x
e=np.random.rand(m,1)
b=p-e*0.01
x_hat=np.linalg.inv((A.T@A))@(A.T@b)
p=A@x_hat
fig = plt.figure(figsize = (8, 3))
ax = fig.add_subplot(121)
ax.scatter(b,p)
ax.set_title('small noise:linear combination ')
ax.set_xlabel("b", size = 10, color = "r")
ax.set_ylabel("p", size = 10, color = "r")
b=p-e
x_hat=np.linalg.inv((A.T@A))@(A.T@b)
p=A@x_hat
ax2 = fig.add_subplot(122)
ax2.scatter(b,p)
ax2.set_title('large noise')
ax2.set_xlabel("b", size = 10, color = "r")
ax2.set_ylabel("p", size = 10, color = "r")
plt.show()
Python3ではじめるシステムトレード【第2版】環境構築と売買戦略
「画像をクリックしていただくとpanrollingのホームページから書籍を購入していただけます。