[Python] 最小自乗法で重回帰と単回帰を求めてみた
2018-04-25

最小自乗法による回帰分析の練習をしました。

最小自乗法(最小二乗法)は分布した点のばらつきが最も小さくなるような関数を求める手法です。 残差の自乗和が最小になる関数を求めます。

単回帰

単回帰とは 1つの説明変数から 目的変数を求める回帰分析です。 ここでは 以下の方法で求める

  • 偏導関数の連立方程式を解く
  • 共分散と分散から回帰係数を算出して解く

データは適当に考えた身長と体重のデータを使います。

In [1]:
%autosave 0
%matplotlib inline
Autosave disabled
In [2]:
import numpy as np
from IPython.display import display
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy import symbols, Eq, init_printing, diff, solve
init_printing()
In [3]:
a, b, X, Y, Cov, σ, σx, σy, μx, μy = symbols('a b X Y Cov σ σ_x σ_y μ_x μ_y')
In [4]:
data = np.array(sorted([  # ソートしたほうが使いやすいのでさきにやっておく
    [172, 58],
    [180, 80],
    [175, 60],
    [162, 50],
    [167, 63],
    [151, 45],
    [188, 70],
    [170, 70],
    [166, 61],
    [100, 20],
    [110, 30],
    [90, 25],
]))
In [5]:
plt.gca().set_aspect('equal', adjustable='box')
plt.plot(data[:,0], data[:,1], 'ro')
plt.show()
In [6]:
d = 0
for x, y in data:
    d += (y - (a * x + b)) ** 2

d.expand()
Out[6]:
$$291543 a^{2} + 3662 a b - 206056 a + 12 b^{2} - 1264 b + 37304$$
In [7]:
# a, b それぞれについて 偏微分 する
expr_a = diff(d, a)
expr_b = diff(d, b)
expr_a, expr_b
Out[7]:
$$\left ( 583086 a + 3662 b - 206056, \quad 3662 a + 24 b - 1264\right )$$
In [8]:
# 偏微分したものを連立して解く
solved = solve([expr_a, expr_b])
solved
Out[8]:
$$\left \{ a : \frac{79144}{145955}, \quad b : - \frac{4389092}{145955}\right \}$$
In [9]:
xs = np.linspace(data[0,0], data[-1,0], 2)

plt.plot(xs, solved[a] * xs + solved[b])
plt.plot(data[:,0], data[:,1], 'ro')
plt.show()
In [10]:
# 今度は
# 共分散 (偏差の積和の平均) と 分散 (偏差の2乗和の平均) から 回帰係数 を求める
# 相関係数ではないので注意 (近い値になることはあるけど)
Eq(a, Cov(X, Y) / (σx ** 2)), Eq(b, μy - a * μx)
Out[10]:
$$\left ( a = \frac{1}{σ_{x}^{2}} \operatorname{Cov}{\left (X,Y \right )}, \quad b = - a μ_{x} + μ_{y}\right )$$
In [11]:
# 共分散行列を求める
# Sxx Sxy 
# Sxy Syy
cov = np.cov(data[:,0], data[:,1], bias=True)
cov
Out[11]:
array([[1013.57638889,  549.61111111],
       [ 549.61111111,  334.88888889]])
In [12]:
# 共分散 / Xの分散
gradient = cov[0, 1] / cov[0, 0]
gradient
Out[12]:
$$0.5422493234216025$$
In [13]:
# 傾きから Y切片 を求める
# Yの平均 - 傾き * Xの平均
intercept = data[:, 1].mean() - gradient * data[:, 0].mean()
intercept
Out[13]:
$$-30.071542598746184$$
In [14]:
xs = np.linspace(data[0,0], data[-1,0], 2)

plt.plot(xs, gradient * xs + intercept)
plt.plot(data[:,0], data[:,1], 'ro')
plt.show()

重回帰

単回帰とは 複数の説明変数から 目的変数を求める回帰分析です。 ここでは 以下の方法で求める

  • 偏導関数の連立方程式を解く (先程と同じ)
  • 擬似逆行列を使って解く

先程の身長、体重に腹囲を追加したデータを使います。

In [1]:
%autosave 0
%matplotlib inline
Autosave disabled
In [2]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy import (
    symbols, Eq, init_printing, diff, solve, Symbol, 
    Matrix, MatrixSymbol
)
init_printing()
In [3]:
# さっきのデータに 腹囲 を足してみる。値は適当。
data = np.array([
    [172, 58, 63],
    [180, 80, 85],
    [175, 60, 70],
    [162, 50, 65],
    [167, 63, 70],
    [151, 45, 59],
    [188, 70, 90],
    [170, 70, 70],
    [166, 61, 66],
    [100, 20, 40],
    [110, 30, 50],
    [90, 25, 30],
])
In [4]:
a1, a2, b = symbols('a_1 a_2 b')
In [5]:
d = 0
for x, y, z in data:
    d += (z - (a1 * x + a2 * y + b)) ** 2

d.expand()
Out[5]:
$$291543 a_{1}^{2} + 206056 a_{1} a_{2} + 3662 a_{1} b - 242982 a_{1} + 37304 a_{2}^{2} + 1264 a_{2} b - 86490 a_{2} + 12 b^{2} - 1516 b + 51056$$
In [6]:
expr_a1 = diff(d, a1)
expr_a2 = diff(d, a2)
expr_b = diff(d, b)
expr_a1, expr_a2, expr_b
Out[6]:
$$\left ( 583086 a_{1} + 206056 a_{2} + 3662 b - 242982, \quad 206056 a_{1} + 74608 a_{2} + 1264 b - 86490, \quad 3662 a_{1} + 1264 a_{2} + 24 b - 1516\right )$$
In [7]:
solved = solve([expr_a1, expr_a2, expr_b])
solved
Out[7]:
$$\left \{ a_{1} : \frac{2279285}{8070429}, \quad a_{2} : \frac{23472007}{64563432}, \quad b : \frac{14954299}{16140858}\right \}$$
In [8]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('Height')
ax.set_ylabel('Weight')
ax.set_zlabel('Waist circumference')

xs = np.linspace(data[:,0].min(), data[:,0].max(), 2)
ys = np.linspace(data[:,1].min(), data[:,1].max(), 2)
zs = solved[a1] * xs + solved[a2] * ys + solved[b]
ax.plot(xs, ys, zs)
ax.scatter(data[:,0], data[:,1], data[:,2])
Out[8]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f12f4444dd8>
In [9]:
# 続いて 擬似逆行列を用いて算出する
num = len(data)
Y = MatrixSymbol('Y', num, 1)
X = MatrixSymbol('X', num, 3)
A = MatrixSymbol('A', 3, 1)
In [10]:
xs = data[0:, 0]
ys = data[0:, 1]
zs = data[0:, 2]
# pd.DataFrame([(x, y, z) for x, y, z in zip(xs, ys, zs)], columns=['x', 'y', 'z'])
In [11]:
# X を [xs, ys, 1]
xm = Matrix(np.array([xs, ys, [1] * num]).T)
xm
Out[11]:
$$\left[\begin{matrix}172 & 58 & 1\\180 & 80 & 1\\175 & 60 & 1\\162 & 50 & 1\\167 & 63 & 1\\151 & 45 & 1\\188 & 70 & 1\\170 & 70 & 1\\166 & 61 & 1\\100 & 20 & 1\\110 & 30 & 1\\90 & 25 & 1\end{matrix}\right]$$
In [12]:
# A を [a1, a2, b]
am = Matrix([a1, a2, b])
am
Out[12]:
$$\left[\begin{matrix}a_{1}\\a_{2}\\b\end{matrix}\right]$$
In [13]:
# Y を zs
ym = Matrix(zs)
ym
Out[13]:
$$\left[\begin{matrix}63\\85\\70\\65\\70\\59\\90\\70\\66\\40\\50\\30\end{matrix}\right]$$
In [14]:
# とすると 以下のように表せ、A について求めればそれぞれの傾きと切片が求まる
Eq(Y, X * A)
Out[14]:
$$Y = X A$$
In [15]:
# こういうこと
Eq(ym, xm * am)
Out[15]:
$$\left[\begin{matrix}63\\85\\70\\65\\70\\59\\90\\70\\66\\40\\50\\30\end{matrix}\right] = \left[\begin{matrix}172 a_{1} + 58 a_{2} + b\\180 a_{1} + 80 a_{2} + b\\175 a_{1} + 60 a_{2} + b\\162 a_{1} + 50 a_{2} + b\\167 a_{1} + 63 a_{2} + b\\151 a_{1} + 45 a_{2} + b\\188 a_{1} + 70 a_{2} + b\\170 a_{1} + 70 a_{2} + b\\166 a_{1} + 61 a_{2} + b\\100 a_{1} + 20 a_{2} + b\\110 a_{1} + 30 a_{2} + b\\90 a_{1} + 25 a_{2} + b\end{matrix}\right]$$
In [16]:
# ただし、Xは正方行列でないため 転置行列を両辺にかけて正方行列にする
Eq(X.T * Y, X.T * X * A)
Out[16]:
$$X^T Y = X^T X A$$
In [17]:
# A を括りだすために XTX の逆行列をかける
Eq((X.T * X) ** -1 * (X.T * Y), A)
Out[17]:
$$\left(X^T X\right)^{-1} X^T Y = A$$
In [18]:
# (A以外を)代入すると..
a = (xm.T * xm).inv().dot(xm.T * ym)
a
Out[18]:
$$\left [ \frac{2279285}{8070429}, \quad \frac{23472007}{64563432}, \quad \frac{14954299}{16140858}\right ]$$
In [19]:
# aを用いて線を描画する
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('Height')
ax.set_ylabel('Weight')
ax.set_zlabel('Waist circumference')

xs = np.linspace(data[:,0].min(), data[:,0].max(), 2)
ys = np.linspace(data[:,1].min(), data[:,1].max(), 2)
zs = a[0] * xs + a[1] * ys + a[2]
ax.plot(xs, ys, zs)
ax.scatter(data[:,0], data[:,1], data[:,2])
Out[19]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f12f6559eb8>

二次関数の回帰

上記はいずれも一次関数をモデルとしていましたが、それ以外の関数にフィッティングさせることもできます。

今回は近似する二次関数を求めてみます。先程使った 擬似逆行列 で求めてみます。データは適当。

In [1]:
%autosave 0
%matplotlib inline
Autosave disabled
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from sympy import (
    symbols, Eq, init_printing,
    Matrix, MatrixSymbol
)
init_printing()
In [3]:
y, x, a0, a1, a2 = symbols('y x a0, a1 a2')
In [4]:
# 適当に座標を作る
data = np.array(sorted([
    [-5, 1000],
    [-4, 750],
    [-3, 600],
    [-2, 300],
    [-1, 70],
    [0, 0],
    [1, 50],
    [2, 200],
    [3, 500],
    [4, 800],
    [5, 1000],
    [6, 1500],
]))
In [5]:
xs = data[0:, 0]
ys = data[0:, 1]
In [6]:
num = len(xs)
In [7]:
Y = MatrixSymbol('Y', num, 1)
X = MatrixSymbol('X', num, 3)
A = MatrixSymbol('A', 3, 1)
In [8]:
# 二次関数の多項式の係数をそれぞれ a0, a1, a2 とする
Eq(y, a0 + a1 * x + a2 * x ** 2)
Out[8]:
$$y = a_{0} + a_{1} x + a_{2} x^{2}$$
In [9]:
# X を [1, x, x^2] の行列
xm = Matrix(np.array([[1] * num, xs, xs ** 2]).T)
xm
Out[9]:
$$\left[\begin{matrix}1 & -5 & 25\\1 & -4 & 16\\1 & -3 & 9\\1 & -2 & 4\\1 & -1 & 1\\1 & 0 & 0\\1 & 1 & 1\\1 & 2 & 4\\1 & 3 & 9\\1 & 4 & 16\\1 & 5 & 25\\1 & 6 & 36\end{matrix}\right]$$
In [10]:
# A を [a0, a1, a2]
am = Matrix([a0, a1, a2])
am
Out[10]:
$$\left[\begin{matrix}a_{0}\\a_{1}\\a_{2}\end{matrix}\right]$$
In [11]:
# Y を ys
ym = Matrix(ys)
ym
Out[11]:
$$\left[\begin{matrix}1000\\750\\600\\300\\70\\0\\50\\200\\500\\800\\1000\\1500\end{matrix}\right]$$
In [12]:
# とすると 先程と同様に それぞれの係数を求められる
Eq(Y, X * A)
Out[12]:
$$Y = X A$$
In [13]:
# こういうこと
Eq(ym, xm * am)
Out[13]:
$$\left[\begin{matrix}1000\\750\\600\\300\\70\\0\\50\\200\\500\\800\\1000\\1500\end{matrix}\right] = \left[\begin{matrix}a_{0} - 5 a_{1} + 25 a_{2}\\a_{0} - 4 a_{1} + 16 a_{2}\\a_{0} - 3 a_{1} + 9 a_{2}\\a_{0} - 2 a_{1} + 4 a_{2}\\a_{0} - a_{1} + a_{2}\\a_{0}\\a_{0} + a_{1} + a_{2}\\a_{0} + 2 a_{1} + 4 a_{2}\\a_{0} + 3 a_{1} + 9 a_{2}\\a_{0} + 4 a_{1} + 16 a_{2}\\a_{0} + 5 a_{1} + 25 a_{2}\\a_{0} + 6 a_{1} + 36 a_{2}\end{matrix}\right]$$
In [14]:
# 解き方の復習
Eq(A, (X.T * X) ** -1 * (X.T * Y))
Out[14]:
$$A = \left(X^T X\right)^{-1} X^T Y$$
In [15]:
# それぞれ x^0, x^1, x^2 の係数
a = (xm.T * xm).inv().dot(xm.T * ym)
a
Out[15]:
$$\left [ \frac{23165}{286}, \quad - \frac{5605}{2002}, \quad \frac{79735}{2002}\right ]$$
In [16]:
# aを用いて線を描画する
plt.plot(xs, a[0] + a[1] * xs + a[2] * xs ** 2)
plt.plot(data[:,0], data[:,1], 'ro')
plt.show()

備考

今回は理解のために擬似逆行列を用いましたが、最小二乗法の解法としてはあまり用いられないようです。

参考