python

pythonのnumpyだけで重回帰をフルスクラッチ実装!!

こんにちは。

みなさん、重回帰してますか?

僕は研究でよくしてます。計算が速い上に安定してそこそこの精度を出すRidge回帰が好きです。

通知表で38/45をとるそこそこの優等生として見てます。

Lasso回帰は多くのパラメータを0にするので説明するのに役立つ、と巷では言われていますが、個人的には精度や、なぜその変数を0にしたのかを、万人が納得できるように説明するのが難しいので好きではありません。

今回のブログでは、これらの回帰の元となる重回帰をスクラッチで実装をし、その続編としてRidge回帰, Lasso回帰をスクラッチで実装していきます。

重回帰とは?

文字の説明

  • $y$ : 目的変数 ( 従属変数 )
  • $x$: 説明変数 ( 独立変数 )
  • $ϕ_j$ : 基底関数 ( $x_j$, $x_j^2$, $ln{(x_j)}$ など )
  • $ω_j$ : 回帰係数
  • $ε$ : 撹乱項、誤差項

どんな考え方をするの?

数学的な部分はものすごく軽く触れます。

重回帰は独立変数が複数存在するモデルを使った回帰のことです。

$y = f(x_1,x_2,…,x_n)$

超簡単に説明をすると、上記の式に入力の$x_n$を用いて、$y$をもとめるための$f(x)$という良い感じの式を求めることです。

$y =ω_0+ω_1 ϕ_1 (x)+ω_2 ϕ_2 (x)+⋯+ω_{K-1} ϕ_{K-1} (x)+ε$
   $=w^T ϕ+ε$

ただし、

$w= (ω_0,ω_1,⋯,ω_{K-1})^T,   ϕ=(ϕ_0 (x),ϕ_1 (x) , ⋯,ϕ_{K-1} (x))^T$

とします。

すなわち、&y&を求めるために必要なものは、

  1. 重みとなる$ω$
  2. 独立関数$ϕ$

この二つです。

$ϕ$については、$x$の何乗までを入力にしてあげるかを関数の中で決めれるようにします。

なので、ここでは解説を飛ばします。

どうやってωを求めるの?

良い感じの$f(x)$について考えるために誤差$E$を調べます。

$\displaystyle E = \frac{1}{2}\sum^N_{n=1}( f(x_n) – y_n )$

$f(x_n) = w_n^T ϕ_n$と表すことができるので、

$E = \frac{1}{2}||w^T ϕ – y || $

$E$が最小になってくれると嬉しいので、最小二乗法を用いて、$E$が最小となる$ω$を計算します。

ちなみに、$E$についている$\frac{1}{2}$ですが、微分をした後に良い感じの式になってほしいからついている係数です。本質的に意味は無いので無視して構いません。

最初値を求めたいので、$E$を$ω$で微分した結果が0になればよいです。微分を計算をして$ω$について整理すると、

$ω=(ϕ^T ϕ )^{-1} ϕ^T y$

となります。

これで、$w$を求めるコトができました。

$y$を求める

$y=w^T ϕ$

と適応をしてあげれば完成です。

重回帰の実装

以下、実装をしていきます。

使い勝手を考えてクラスとして定義をします。

fitの際指定をする、dimで入力の$x$を何乗したものまでを入力変数にするのかを決定します。

ちなみに、fitの部分では一番最初のインデックスには切片のための全て1の列が入っています。

predictの部分ではfitの際取得した切片、self.intercept_があるので、全て1の列は削除しています。

個人的な頑張ったポイントとしては、reshapeによる挙動を確かめながら、どの変数がどの列に対応するのかを対応づけたところです。皆さんも、reshapeを使用する際は細心の注意を払うことをオススメします。

import numpy as np

class MyMultipleRegressor:
    def __init__(self):
        self.dim = None
        self.coef_ = None
        self.intercept_ = None
        self.n_feature = None
        self.__name__ = "Multiple"

    def fit(self, X, y, dim=1):
        """
        Calculate Linear Regression
        Parameters
        X : array_like (n_sample, n_feature)
            The matrix of explanatory variables
        y : array_like (n_sample, 1)
            Objective variable
        dim : int, default 1
            Dimension of multiple regression
        """
        n_sample, n_feature = X.shape
        # X (n_sample, dim, n_feature)
        # make basis function φ
        X = np.array([[x**i for i in range(dim + 1)] for x in X])
        # delete duplicated ones array
        # ones array is Intercept
        X = X.reshape(n_sample, -1)[:, n_feature - 1:]
        # calculate coefficient
        weight = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), y)
        self.coef_ = weight
        self.intercept_ = weight[0]
        self.dim = dim
        self.n_feature = n_feature

    def predict(self, X):
        """
        Calculate predicted the result of multiple linear regression
        Parameters
        X : array_like (n_sample, n_feature)
            The matrix of explanatory variables in test dataset
        Returns
            Predicted objective variable
        """
        n_sample, n_feature = X.shape
        # X (n_sample, dim, n_feature)
        X = np.array([[x**i for i in range(self.dim + 1)] for x in X])
        # delete ones array
        X = X.reshape(n_sample, -1)[:, n_feature:]
        return (np.dot(X, self.coef_) + self.intercept_)

    def mse(self, X, y):
        """
        Calculate mean square error 
        Parameters
        X : array_like (n_sample, n_feature)
            The matrix of explanatory variables in test dataset
        y : Ground truth of objective variable
        Returns
            Mean square error
        """
        return np.mean((self.predict(X) - y)**2)

    def title(self):
        """
        Make title for plot
        Return
            string of equation
        """
        Xn = [f"x_{k}" for k in range(1, self.n_feature + 1)]
        title = "$f(x) = $" + f"{float(self.intercept_): .2f}"
        for d in range(1, self.dim + 1):
            for k in range(self.n_feature):
                if self.coef_[(d - 1) * self.n_feature + k] > 0:
                    sign = "+"
                elif self.coef_[(d - 1) * self.n_feature + k] < 0:
                    sign = "-"
                else:
                    continue
                if d == 1:
                    title += f"${sign}{abs(float(self.coef_[(d-1) * self.n_feature + k])):.2f}{Xn[k]}$"
                else:
                    title += f"${sign}{abs(float(self.coef_[(d-1) * self.n_feature + k])):.2f}{Xn[k]}^{d}$"
        return title

 

titleメソッドはplotをするときに式をきれいに表示させるために作成をしました。

ここも頑張ったので、みんなに使ってもらいたい!笑

データの確認

今回は以下の3種類のデータに対して重回帰をしていきます。

data3については3次元のデータなのでぱっと見ではどんな係数にすればよいのか見分けがつきにくいです。

重回帰を適応

コードを確認

それぞれのデータは、 data1.csv, data2.csv, data3.csvと保存されています。

import argparse
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from sklearn.model_selection import KFold

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("-d", "--dim", type=int, default=1, help="回帰の次元数.デフォルトは1")
    parser.add_argument("-f", "--fname", type=str, default="data1.csv", help="読み込むファイルの名前.デフォルトは'data1.csv'")
    args = parser.parse_args()

    # define regression model
    dim = args.dim
    model = MyMultipleRegressor()
    # load data
    fname = args.fname
    data = pd.read_csv('../' + fname)
    # split data
    n_sample, n_cols = data.shape
    feature_cols = list(data.columns[:-1])
    target_col = list(data.columns[-1:])
    kf = KFold(n_splits=2, shuffle=True)
    for train_idx, test_idx in kf.split(data):
        break
    # train data
    train_X = data.loc[train_idx, feature_cols].values.reshape(-1, n_cols - 1)
    train_y = data.loc[train_idx, target_col].values.reshape(-1, 1)
    model.fit(train_X, train_y, dim=dim)
    train_mse = model.mse(train_X, train_y)
    # test data
    test_X = data.loc[test_idx, feature_cols].values.reshape(-1, n_cols - 1)
    test_y = data.loc[test_idx, target_col].values.reshape(-1, 1)
    test_mse = model.mse(test_X, test_y)
    title = model.title()
    if n_cols == 2:
        # plot original data
        plt.figure(figsize=(10, 8))
        plt.xlabel(feature_cols[0], size=16)
        plt.ylabel(target_col[0], size=16)
        plt.title(f"{model.__name__} Regression\n" + title + f"\nTrain MSE : {train_mse:.4f}, Test MSE : {test_mse:.4f}", size=20)
        plt.scatter(train_X, train_y, label="Train Observed")
        plt.scatter(test_X, test_y, label="Test Observed")
        # plot predicted data
        X1 = np.linspace(data[feature_cols].min(), data[feature_cols].max(), 500)
        y1_hat = model.predict(X1)
        plt.plot(X1, y1_hat, c="r", label="Predicted")
        plt.legend(fontsize=16)
        plt.savefig(fname[:-4] + f"{model.__name__}_{dim}D.png", bbox_inches="tight")
        plt.show()
    elif n_cols == 3:
        # plot original data
        fig = plt.figure(figsize=(10, 8))
        ax = Axes3D(fig)
        ax.set_xlabel(feature_cols[0])
        ax.set_ylabel(feature_cols)
        ax.set_zlabel(target_col[0])
        plt.title(f"{model.__name__} Regression\n" + title + f"\nTrain MSE : {train_mse:.4f}, Test MSE : {test_mse:.4f}")
        ax.plot(data.loc[train_idx, feature_cols[0]].values, data.loc[train_idx, feature_cols].values,
                data.loc[train_idx, target_col].values.flatten(),
                marker="o", linestyle='None', label="Train Observed")
        ax.plot(data.loc[test_idx, feature_cols[0]].values, data.loc[test_idx, feature_cols].values,
                data.loc[test_idx, target_col].values.flatten(),
                marker="o", linestyle='None', label="Test Observed")
        # plot predicted data
        X1 = np.linspace(data[feature_cols[0]].min(), data[feature_cols[0]].max(), 30)
        X2 = np.linspace(data[feature_cols].min(), data[feature_cols].max(), 30)
        X1, X2 = np.meshgrid(X1, X2)
        X3 = np.concatenate([X1.reshape(-1, 1), X2.reshape(-1, 1)], axis=1)
        y3_hat = model.predict(X3).reshape(X1.shape)
        ax.plot_wireframe(X1, X2, y3_hat, color="red", label="Predictied")
        plt.legend(loc='center left', fontsize=16)
        plt.savefig(fname[:-4] + f"{model.__name__}_{dim}D.png", bbox_inches="tight")
        plt.show()


if __name__ == "__main__":
    main()

 

こんな感じで実行できます。

pythonでの3次元データのplotも参考になれば幸いです!

結果を確認

data1

見た目通り1次関数でのフィッティングが良い感じですね。

3次関数の方がTrain MSEは小さくなっていますが、Test MSEは大きくなっています。

これがいわゆる、過学習というやつです。

実世界でも、重要なのはTestデータで高い精度が出ることなので、正しい次数の選択が重要であることが分かります。

 

data2

今度は逆に3次関数でのフィッティングが良い感じですね。

1次関数の方ではTrain MSE、Test MSEともに大きくなっています。

これは表現力不足な状態です。

はやり、ほどよい次数選択が重要です。

data3

1~3次元でフィッティングをしてみました。

Test MSEに着目をすると、2次元でのフィッティングが一番良いです。

ここで、2次元でフィッティングした式の係数に着目して、

$f(x) = -0.12+1.95{x_1} – 0.61{x_2} + 0.10{x_1^2}+3.08{x_2^2}$

となっていることから、このデータにとって重要な要素は$x_1$, $x_2^2$であると当たりをつけることができます。

果たして、その解釈は正しいのかを確かめていきたいと思います。

 

係数を参考にしたデータの解釈

左上から順に、$(x1, x2)$, $(x2, x3)$, $(x3, x1)$の関係を表しています。

右下は、これらの相関係数を表したものです。

data3は、$(x2, x3)$が強く関係をしていて、それ以外はほとんど関連がないことが分かります。

 

以上より、$x_1$がなんか重要そうに見えたのは、偶然だということが分かりました。

このような結果になった理由は、恐らく、訓練データとテストデータで分割をする際に、訓練データの方でたまたま$x_1$に相関がある別れ方をしたからだと考えられます。

 

特徴量が1つ増えるだけでも人間にとって理解をすることが難しくなることがわかりました。

最後に

今回は重回帰について実装をしました。

また、データの解釈についても軽く触れてみました。

モデルに任せればオールオッケーという訳ではないと言うことも理解できました。

やはり、人間の解釈が重要となります。

優れたモデルを扱うときほどよくよく胸に刻んでおきたいです。

次回予告

実は、重回帰ではオーバーフィッティングしてしまうことが多いので、それを押えるために生まれたのが、Ridge, Lasso回帰です。

次回はRidge回帰を実装していきます。

こんな感じのブログを書いていきます。

もし興味を持っていただけたら、フォローお願いします。

pythonのnumpyだけでリッジ(Ridge)回帰をフルスクラッチ実装!!今回はRidge回帰を実装していきます! 実はリッジ回帰は重回帰とほとんど同じコードで作成することができるのです! なので、...
pythonのnumpyだけでラッソ(Lasso)回帰をフルスクラッチ実装!今回は、 Lasso回帰は多くのパラメータを0にするので説明するのに役立つ と巷では言われるLasso回帰を実装して...

 

オススメのプログラミングスクールをご紹介

タイピングもままならない完全にプログラミング初心者から

アホいぶきんぐ
アホいぶきんぐ
プログラミングってどこの国の言語なの~?

たった二ヶ月で

いぶきんぐ
いぶきんぐ
え!?人工知能めっちゃ簡単にできるじゃん!

応用も簡単にできる…!!

という状態になるまで、一気に成長させてくれたオススメのプログラミングスクールをご紹介します!

テックアカデミーのPython+AIコースを受講した僕が本音のレビュー・割引あり! というプログラミング完全初心者だった僕が Tech Academy(テックアカデミー)のPython×AIコース を二ヶ月間...

COMMENT

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です