ほぼPython

Not技術ブログBut勉強ブログ 内容には誤りがあることが多いです

目で見て楽しむ重回帰分析

前回のブログで、説明変数が1つの場合の回帰分析を行いました。
最適な解を求めるのに、わざわざ勾配法を使いましたが、普通に偏微分方程式を解けば、解析解を得ることが出来ます。
つまり、勾配ベクトルの成分=0として解いていけば、求められます。説明変数が1つの場合はとても簡単に行うことが出来るので飛ばして、
今回は、説明変数が2つの場合で行いたいと思います。

まずはデータを用意して、それがどんなグラフになっているか確認します。今回はy = 3*x0 + 2*x1 + 5 という方程式を思い浮かべながら適当に値に誤差を持たせるという作業を手動で行い、データを用意しました。

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

#データの準備
X0=np.array([0,1,1,2,-1,3,-2,0,2,-4]) #説明変数1
X1=np.array([0,1,3,-1,-5,1,1,3,-5,1]) #説明変数2
T=np.array([5.2,10.4,13.5,8.5,-7.6,15.1,0.8,8.9,2.4,-4.2]) #被説明変数

#グラフ描画
plt.figure(figsize=(5,5))
ax =plt.subplot(1,1,1,projection='3d')
for i in range(len(X0)):
    ax.plot([X0[i],X0[i]],[X1[i],X1[i]],[-7.6,T[i]],color='gray')
ax.plot(X0,X1,T,marker='o',linestyle='None')
ax.set_xlabel("X0")
ax.set_ylabel("X1")
ax.set_zlabel("T")
plt.show()

すると、以下のような3次元のグラフを得ることが出来ます。あるX0とX1により、Tが定まるというグラフです。

f:id:short_4010:20180301184733p:plain

それでは、これらの点すべてにもっとも当てはまりのよい平面の方程式を求めていきたいとおもいます。
平面の方程式を y = w0*x0 + w1*x1 + w2 とします。

#解析解
def fit_plane(x0,x1,t):
    c_tx0 = np.mean(t*x0) - np.mean(t)*np.mean(x0)
    c_tx1 = np.mean(t*x1) - np.mean(t)*np.mean(x1)
    c_x0x1 = np.mean(x0*x1) - np.mean(x0)*np.mean(x1)
    v_x0 = np.var(x0)
    v_x1 = np.var(x1)
    w0 = (c_tx1*c_x0x1 - v_x1 * c_tx0)/(c_x0x1**2 - v_x0*v_x1) 
    w1 = (c_tx0*c_x0x1 - v_x0*c_tx1)/(c_x0x1**2 - v_x0 * v_x1)
    w2 = -w0*np.mean(x0) - w1*np.mean(x1)+np.mean(t)
    return np.array([w0,w1,w2])

たったこれだけで方程式の係数を求められました。結果はw0=2.86, w1=1.76, w2=4.90となっており、まあ僕のイメージしていた方程式に適当に誤差を含めたらこんな感じになるだろうという結果になりました。

それでは、実際に、この方程式が元のデータにどれくらい当てはまっているか可視化してみましょう。

def show_plane(w):
    x0 = np.linspace(-5,5,5)
    x1 = np.linspace(-5,5,5)
    xx0,xx1 = np.meshgrid(x0,x1)
    w0,w1,w2 = w
    y = w0*xx0 + w1*xx1 +w2
    ax.plot_surface(xx0,xx1,y,rstride=1,cstride=1,alpha=0.3,color='blue',edgecolor='black')

という関数を追加して、一番最初のグラフに重ねてみます。
f:id:short_4010:20180301224335g:plain

いい感じにフィッティングしていますね。(見ているのが楽しくて10分くらい何も考えずにぐるぐる動かしてました。)

ということで、与えられたデータから平面の方程式
y = 2.86* x0 + 1.76 * x1 + 4.90 を求められたので与えられていないデータについても予測が行えるようになりました。
病気の治癒率や不動産価格、その他もろもろいくつかの変数が結果に影響を与えていて、それが線形結合の方程式に沿いそうだというとき、使えそうですね。(本当は最小二乗法を使うためには結構強い仮定が必要なのでその仮定に沿っている必要はあるようです、勉強不足なのでまだ詳しくわかりません。)

最後に、解析解のことについて疑問を投げて終わりにしようと思います。解析解は、まず、二乗誤差を求め、その二乗誤差を各係数(w0,w1,w2)で偏微分してそれを=0として解いていくことで、fit_planeの中に書いたような解としてw0,w1,w2が求まります。

と、いろんなところに書いてありますが、偏微分して=0というのは極値としての必要条件であって、極値をとるとは言い切れない気がします。そしてそもそも極値の中でも極小値をとると言い切れないと思います。

しかし、なぜそのようなことがいえるのか謎です。数学に詳しくないので詳しい方教えてください・・・。

追記

なぜ偏微分=0で極小値と言えるのかわかりました。回帰モデルが1次式の場合、残差の二乗和は、かならず、二次式になりますね。
しかも残差二乗和は全体として正なので、下に凸なんですね。だから偏微分=0だけで極小(最小)と言えるんですね。(たぶん)