LoginSignup
4
4

More than 3 years have passed since last update.

3Dアニメーションで見る最小二乗法のこころ

Last updated at Posted at 2019-10-20

はじめに

最小二乗法はなぜ、差分の二乗和を考えるのか?という疑問について答えるのが本稿の目的である。
ここでは、最小二乗法の一般的な説明
1.データと近似関数値の差分の二乗和を作る。
2.上の二乗和をパラメータごとに偏微分して、その共通零点を求める。
には殆どふれずに、上の1.に至るまでの思想について自分なりに記述してみたい。もちろん、最小二乗法はとても広い応用を持つが、本稿では直線近似に限定して、「最小二乗法のこころ」に迫りたい。

3点を直線近似する

平面にいくつかの点を与えたときに、その点を通る直線はどれだけあるか?ということから考え始めてみる。

1点を通る直線は無限個存在する。

image.png

異なる2点を通る直線は1つだけ存在する。

image.png

互いに異なる3点を通る直線は一般には存在しない。

(もちろん、最初から直線上にある3点を選んだとすれば直線が存在する。)
image.png
純粋な幾何学ならば「ここから先は考えない!」と言い切っても良いだろう。ただ、応用的な立場に立つならば「3点をちょうど通る直線は存在しなくても、例えば上図の赤い直線のようにできるだけ3点に近い直線を求めたい。」ということになる。そこで、直線近似の出番である。

素朴に考えてみる

上図の$P_1$と$P_2$に直線を近づけると$P_3$からは遠ざかる。
また、$P_2$と$P_3$に直線を近づけると$P_1$からは遠ざかり、
同様に$P_3$と$P_1$に直線を近づけると$P_2$からは遠ざかる。
どうも$P_1,P_2,P_3$を別々に考えていると、あちらを立てればこちらが立たずとなってしまい、考えが発散しそうになる。また「できるだけ3点に近い直線」とは何か?...等々思いめぐらしてみた。

最小二乗法のこころに迫る

そして、次のことが大事であると思うようになった。
平面上の3点 $P_1,P_2,P_3$を一組の対象(一つの点)として考えること
これは、直線をいい具合に$P_1,P_2,P_3$に近づけるために必要と考えた。そして、これができれば「できるだけ3点に近い直線」を定義できるのではないか?
そこで、座標系を導入して実際に平面上の3点を1点としてみたい。下の左図ように通常の座標系 $(x,y)$ を導入して、$P_1$、$P_2$、$P_3$の座標がそれぞれ $(x_1,y_1)$、$(x_2,y_2)$、$(x_3,y_3)$であるとする。また、今は直線近似したいので、直線の式を $y=f(x)=ax+b$ と表しておく。さらに、$x=x_1$ における直線上の点を $Q_1$、$x=x_2$ における直線上の点を $Q_2$、$x=x_3$ における直線上の点を $Q_3$ とすれば、座標はそれぞれ$(x_1,f(x_1))$、$(x_2,f(x_2))$、$(x_3,f(x_3))$ となる。

grf1.png

$Q_1$、$Q_2$、$Q_3$ が $P_1$、$P_2$、$P_3$ に近づけば、直線が $P_1$、$P_2$、$P_3$ に近づく。 $Q_i (i=1,2,3)$ と $P_i (i=1,2,3)$ を比較すると、$x$ 座標は共通なので、$f(x_i)(i=1,2,3)$ を $y_i(i=1,2,3)$ に近づければよい。そこで、上の右図のように新たに3次元空間(ここではデータベクトル化空間と呼んでおく)に座標系 $(Y_1,Y_2,Y_3)$ を導入する。そして、 $(y_1,y_2,y_3)$ なるベクトルを考えて$\tilde P$ とする。これで平面上の3点を1点にできた。同様に $(f(x_1),f(x_2),f(x_3))$ なるベクトルを考えて $\tilde Q$ とする。さらに、$\tilde P$ と $\tilde Q$ 2点間の距離を $d$ とする。いよいよ、「最小二乗法のこころ」について述べたい。上で定義してきた記号を用いると、それは、「$\tilde P$ と $\tilde Q$ 2点間の距離 $d$ を最小にするような$a,b$を求めること。」となる。
そして、そのような$a,b$に対応する直線 $y=ax+b$を「できるだけ3点に近い直線」と定義する。
ところで、$d$ は三平方の定理を応用すれば下の式で求められる。


d = \sqrt{(f(x_1)-y_1)^2 + (f(x_2)-y_2)^2 + (f(x_3)-y_3)^2} 

そして、$d$を最小にするのは、$d^2$を最小にするのと同じである。すなわち、


d^2 = (f(x_1)-y_1)^2 + (f(x_2)-y_2)^2 + (f(x_3)-y_3)^2 (差分の二乗和!)を最小にするべき。

かくして、冒頭の「1.データと近似関数値の差分の二乗和を作る。」までたどり着いたことになる。

3Dアニメーションで見てみる。

以上で述べた「最小二乗法のこころ」を可視化するために、Pythonで作成した3Dアニメーションを下に添付した。グラフで用いた記号などは上述してあるので、必要に応じて今までの説明を参照願いたい。
左上:$a,b$ に対する $d^2$ の値をプロットした3Dグラフ
右上:データベクトル化空間に、$a,b$ に対する $\tilde Q$ (紫)をプロットした3Dグラフ。赤点は $\tilde P$ を示す。
左下:$a,b$ に対する直線をプロットしたグラフ
ここでは、最初に与える3点の座標は $P_1:(-7, -8), P_2:(4, 3), P_3:(8, 9)$ としてある。
最小二乗法の解は $a=200/181,b=-92/181$ である。また、アニメーションの1コマごとに$a$の値を変化させている。$b$ については最小二乗法の解$b=-92/181$ に固定してある。
そして、最後の1コマは最小二乗法の解 $a=200/181,b=-92/181$ におけるグラフになっている。resulta.gif(上をクリックすると、アニメーションを再生)

Pythonプログラム

上記の3Dアニメーションを描くPythonプログラムのコードは以下のとおり。


# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

f = lambda a,b,x : a*x+b                    #直線の式
sqfy = lambda a,b,x,y:(f(a,b,x)-y)**2       #近似直線とデータの差分の2乗
d2fy = lambda a,b,x,y:sqfy(a,b,x,y).sum()   #近似直線とデータの差分の2乗和

srcdata = [[-7, -8], [4, 3], [8, 9]]        #データ点
ars = np.array(srcdata)
src_x = np.array([ars[0][0],ars[1][0],ars[2][0]]) #データ点のx座標
src_y = np.array([ars[0][1],ars[1][1],ars[2][1]]) #データ点のy座標

#グラフ準備
fig = plt.figure(figsize=4*plt.figaspect(1.0))
plt.subplots_adjust(hspace=0.2)
ax1 = fig.add_subplot(2, 2, 1, projection='3d') #a,b に対する d2 の値をプロット(3Dグラフ)
ax2 = fig.add_subplot(2, 2, 2, projection='3d') #a,b に対する Q~ (紫)をプロット(3Dグラフ)。固定点としては P~ も描画する。
ax3 = fig.add_subplot(2, 2, 3)                  #a,b に対する直線をプロットしたグラフ

imgs = []   #アニメーションのコマの配列
#a,bの値を設定
#最小二乗法の解は、a=200/181,b=-92/181
#最後のデータは、最小二乗法の解に対応するデータとする
#今回は理解しやすいようにbの値は1つに固定しておく
ao=np.arange(0.1,2.1,0.1)
a=np.insert(ao,20,[200/181])
bv=-92/181
b=bv*np.ones(21)

#アニメーションのコマで使うデータ配列
ys0=[]
ys1=[]
ys2=[]
d2=[]
fy0=[]
fy1=[]
fy2=[]

for ma in a:
    ys0.append(src_y[0])
    ys1.append(src_y[1])
    ys2.append(src_y[2])
    d2.append(d2fy(ma,bv,src_x,src_y))
    fy=f(ma,bv,src_x)
    fy0.append(fy[0])
    fy1.append(fy[1])
    fy2.append(fy[2])

# 各グラフのタイトルや軸ラベル
ax1.set_title("Square of distance depends on a&b",fontsize=18)
ax1.set_xlabel("a",fontsize=15)
ax1.set_ylabel("b",fontsize=15)
ax1.set_zlabel("d^2",fontsize=12)

ax2.set_title("(f(x1),f(x2),f(x3)) depends on a&b",fontsize=18)
ax2.set_xlabel("Y1",fontsize=15)
ax2.set_ylabel("Y2",fontsize=15)
ax2.set_zlabel("Y3",fontsize=15)

ax3.set_title("Approximate line depends on a&b",fontsize=18)
ax3.set_xlabel("x",fontsize=15)
ax3.set_ylabel("y",fontsize=15)
ax3.grid(which = "major", axis = "x", color = "green", alpha = 0.8,linestyle = "--", linewidth = 1)
ax3.grid(which = "major", axis = "y", color = "green", alpha = 0.8,linestyle = "--", linewidth = 1)

for i in range(len(a)):
    if (i < len(a)-1 ):
        #軌跡を描画
        img1 = ax1.plot(a[:i+1], b[:i+1], d2[:i+1],marker=".",color="blue")    
        img2 = ax2.plot(ys0[:i+1],ys1[:i+1],ys2[:i+1],marker="o",color="red")
        img2 = ax2.plot(fy0[:i+1], fy1[:i+1],fy2[:i+1], marker=".",color="violet")
    else:
        #最後の1コマは最小二乗法の解について描画
        img1 = ax1.plot([a[i]], [b[i]], [d2[i]],marker="o",color="blue")    
        img2 = ax2.plot([ys0[i]],[ys1[i]],[ys2[i]],marker="o",color="red")
        img2 = ax2.plot([fy0[i]], [fy1[i]],[fy2[i]], marker="o",color="purple")

    #a,bの値に対応する直線を描画する
    lptX=[]
    lptY=[]
    lptX.append(0)  
    lptX.append(src_x[0])
    lptX.append(src_x[2])
    lptY.append(bv)  
    lptY.append(fy0[i])
    lptY.append(fy2[i])

    img3 = ax3.plot(src_x[0],src_y[0],marker="o",color="red")
    img3 = ax3.plot(src_x[1],src_y[1],marker="o",color="red")
    img3 = ax3.plot(src_x[2],src_y[2],marker="o",color="red")
    img3 = ax3.plot(lptX,lptY,color="violet")

    imgs.append(img1+img2+img3) #コマの情報を追加

anim = animation.ArtistAnimation(fig, imgs, interval=500) #アニメーションの作成
anim.save('resulta.gif', 'imagemagick') #gifファイルに保存

plt.show()

(付録)データベクトル化空間における幾何学と相関係数

上のアニメーションでは $b$ を固定していたが、$a,b$ ともに独立に変化させると、$ \tilde Q(a,b):=(-7a+b,4a+b,8a+b)$ 全体が作る軌跡は平面($H$ と表記)になる。$d$ が最小値 $d_{min}$になるとき、$ \tilde Q$ は $ \tilde P$ から平面$H$ におろした垂線の足になる。これは数式によって確かめることができる。$O$ をデータベクトル化空間の原点、$<,>$ を標準的な内積として、実際に最小二乗法を実行してみると、


\begin{eqnarray}
\vec{q}(a,b)&:=& \vec{O\tilde Q}\\\vec{p} &:=& \vec{O\tilde P}\\
d^2&=&<\vec{q}(a,b)-\vec{p},\vec{q}(a,b)-\vec{p}>\\
&=& <\vec{q}(a,b),\vec{q}(a,b)>-2<\vec{q}(a,b),\vec{p}>+<\vec{p},\vec{p}>\\
\frac{\partial (d^2)}{\partial a} &=& 2<\frac{\partial \vec{q}(a,b)}{\partial a},\vec{q}(a,b)> - 2<\frac{\partial \vec{q}(a,b)}{\partial a},\vec{p}>\\
&=& 2<\frac{\partial \vec{q}(a,b)}{\partial a},\vec{q}(a,b)-\vec{p}>\\
\frac{\partial (d^2)}{\partial b} &=& 2<\frac{\partial \vec{q}(a,b)}{\partial b},\vec{q}(a,b)-\vec{p}>\\
\frac{\partial (d^2)}{\partial a} &=& \frac{\partial (d^2)}{\partial b} = 0\\
&\Leftrightarrow& <\frac{\partial \vec{q}(a,b)}{\partial a},\vec{q}(a,b)-\vec{p}> = 0 ,<\frac{\partial \vec{q}(a,b)}{\partial b},\vec{q}(a,b)-\vec{p}> = 0\\
\therefore \space\space \frac{\partial \vec{q}(a,b)}{\partial a} &\perp& (\vec{q}(a,b)-\vec{p}),\frac{\partial \vec{q}(a,b)}{\partial b} \perp (\vec{q}(a,b)-\vec{p})
\end{eqnarray}

が一般に成り立つ。もし、$\phi:(a,b) \rightarrow \tilde {Q}(a,b)$ が埋め込み写像ならば、$\frac{\partial \vec{q}(a,b)}{\partial a}$ 、$\frac{\partial \vec{q}(a,b)}{\partial b}$ は多様体の$\tilde {Q}(a,b)$ における接平面を張るので、$\vec{q}(a,b)-\vec{p}$ は接平面に直交する。今の場合、$\phi$ は3次元空間の中への埋め込みで、平面 $H$ に写るので、$\vec{q}(a,b)-\vec{p}$ は平面 $H$ に直交する。下図参照。vecspc.png

また、$d$ が最小値 $d_{min}$になるとき、$\bar x$ を$x_i (i=1,2,3)$ の平均値、$\bar y$ を$y_i (i=1,2,3)$ の平均値とおくと、$a,b$ は下のように表せる。

\begin{eqnarray}
a &=& \frac{\sum_{i=1}^{3}(x_i - \bar x)(y_i - \bar y)}{\sum_{i=1}^{3}(x_i - \bar x)^2}\\
b &=& \bar y - a\bar x
\end{eqnarray}

上の式を用いて、さらに${d_{min}}^2$ を計算すると、下のようになり、相関係数と関係があることが分かる。

\begin{eqnarray}
{d_{min}}^2 &=& (\sum_{i=1}^{3}(y_i - \bar y)^2)(1-r^2)\\
ただし、r &=& \frac{\sum_{i=1}^{3}(x_i - \bar x)(y_i - \bar y)}{\sqrt{(\sum_{i=1}^{3}(x_i - \bar x)^2)(\sum_{i=1}^{3}(y_i - \bar y)^2)}} ({x_i}と{y_i}の相関係数)
\end{eqnarray}

このように、最小二乗法のこころには高次元の幾何学が隠れていて、統計量などとも関連を持つ事が分かった。
ここでは、3点を与えたのでデータベクトル化空間は3次元だったが、一般に$N$ 個の点を与えらえた場合にはデータベクトル化空間は$N$次元になる。そして、直線近似を他の関数による近似に置き換えることで、より豊かな高次元の幾何学が現れるのではないか?と期待して、本稿を終わりにする。

参考

1では、アニメーションで軌跡を描く方法を参考にさせていただいた。
2では、複数のグラフを一つのアニメーションに含める方法を参考にさせていただいた。

1.matplotlib で線分をアニメーションさせる
2.matplotlib アニメーション (Jupyter)

4
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4
4