最小二乗法による単回帰分析

最小二乗法による線形解析は、多変数の場合は重回帰分析、1変数の場合単回帰分析と呼ばれます。単回帰分析は、つまり1次関数の式$y=ax+b$です。

回帰分析なのでフィットする直線を求めたい、つまり$a$と$b$の最適値を求めることになります。xとyはデータとして与えられています。

やり方を調べていて、scikit-learnとかのライブラリを使ったり、共分散などを計算してガチに解いたりといろいろありますが、最小二乗法なので簡単にできるのでは？と思いました。それを試します。

知りたい$a,b$をベクトルとして

  A = \left(
    \begin{array}{c}
      a \\
      b \\
    \end{array}
  \right)

と書きます。データ群$X,Y$もベクトルで表せますが、ここで行列（連立方程式）の形にしたいので、$X$は以下のように書きます。

  X = \left(
    \begin{array}{cc}
      x_1 & 1\\
      x_2 & 1\\
      \vdots \\
      x_n & 1
    \end{array}
  \right)

とするのがポイントです。つまり連立方程式は

XA = Y \\  
\left(
    \begin{array}{cc}
      x_1 & 1\\
      x_2 & 1\\
      \vdots \\
      x_n & 1
    \end{array}
  \right)\left(
    \begin{array}{c}
      a \\
      b \\
    \end{array}
  \right)
=
\left(\begin{array}{c}
      y_1\\
      y_2\\
      \vdots \\
      y_n
    \end{array}\right)

ですね。次元解析するとN×2・2×1=N×1 となっていることがわかります。あとは一般化逆行列$X^\dagger$を使えば

  A = X^\dagger Y

と一発で解くことができます。

このようにすると実装も簡単です。以下にテスト用コードを張ります。

import numpy as np
import random

def linear_regression():

    # 答えを最初に作っておく
    a = 2
    b = 10
    x_true = np.arange(0.0,50.0,1.0)
    y_true = a * x_true + b


    # 正解値からランダムにずらしてデータを作る
    xd = np.zeros(len(x_true))
    yd = np.zeros(len(y_true))
    for i in range(len(xd)):
        yd[i] = y_true[i] + 100*(random.randint(-10,10)/100)
    for i in range(len(xd)):
        xd[i] = x_true[i] + 10*(random.randint(-10,10)/100)
    print(xd)
    print(yd)

    # データ群行列
    X = np.c_[xd, np.linspace(1,1,len(xd))]
    print(X)

    # 最小二乗法　一般化逆行列を左から掛けるだけ
    A = np.linalg.pinv(X) @ yd
    y_est = x_true * A[0] + A[1]

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots()
    ax.scatter(xd, yd, label='data')
    ax.plot(x_true,y_true, label='true')
    ax.plot(x_true,y_est, label='linear regression')
    ax.legend()
    ax.grid()
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    fig.tight_layout()
    plt.show()

    return

if __name__ == '__main__' :
    linear_regression()

結果は次のようになります。

trueは元になる直線です。そこからランダムに値を変えているので一致はしませんが、合ってそうです。これならライブラリなどに頼らなくても簡単に回帰分析ができそうです。