ヤナギの木の下に人影が…！

f:id:tatamiyatamiyatatatamiya:20190323175407p:plain

ということで，怪奇…もとい，回帰アルゴリズムを実装していきたいと思います。*1

前回までは決定木・ランダムフォレストの分類アルゴリズムを実装していきました。今回は連続的な値を予測する回帰を扱えるように拡張を行います。

指針

前回までのあらすじと今回やること

今まで扱ってきた分類問題では，不純度が小さくなるようなデータ分割方法を探索し，学習モデルを構築していました。

一方で，株価の値動きなどのように連続的な値をとるような場合でも，不純度の代わりに誤差を用いることで，今までと同様に学習モデルを作ることができます。

最小化する指標

分類問題では，データセット内に含まれる各クラスの構成率 $p_k (k \in \{1,2,3,...,K\})$ をもとに計算したGini係数

$\sum _ {k=1} ^{K} (1-p_k) p_k = 1 - \sum _{k=1} ^{K} p_k ^2$

の最小化を行いました。*2

一方回帰では，分割を行なった左右のデータ断片について，それぞれ平均値からの逸脱度を計算します。そして左右の逸脱度の和が最小になるように分割を行なっていきます。

逸脱度としては，ここでは二乗和を用いることにします。*3

分割の数式表現

N個の学習データセット $(X_i, y_i)_{i \in \{1,2,...,N\}}$ があるとします。 $X_i = (x_i^{(1)}, x_i^{(2)}, ..., x_i^{(D)})$ はD次元の特徴量ベクトルです。

この時，二乗和による評価指標は以下のようにかけます：

$E(j,s) = \sum_{i \in R_{\rm left}(j,s)} (y_{i} - c_{\rm left})^{2} + \sum_{i \in R_{\rm right}(j,s)} (y_{i} - c_{\rm right})^{2}$

ここで， $R_{\rm left/right}(j,s)$ は分割後左or右断片に入るデータインデックスの集合を， $c_{\rm left/right}$ は各断片内での平均値を表します：

$R_{\rm left}(j,s) = \{ i | x_i^{(j)} \le s \}, \quad R_{\rm right}(j,s) = \{ i | x_i^{(j)} > s \}$

$c_{\rm left/right} = \frac{1}{N_{\rm left/right}}\sum_{i \in R_{\rm left/right}(j,s)} y_{I}, \quad N_{\rm left/right} = \sum_{i \in R_{\rm left/right}(j,s)}1$

$j$ は分割に用いる特徴量のインデックス， $s$ は分割の閾値を表します。

最適な分割は，上記の評価指標 $E(j,s)$ が最小となる $(j,s)$ の組み合わせとして表されます。

実装

サンプルデータ

分類問題ではサンプルデータとしてirisを用いて，葉・萼片の長さと幅から品種を予測していました。

回帰問題を扱うにあたっては，同じくよく使われるデータセットとして， boston housingを用います。

データの説明については，こちらを参照願います。

scikit-learnから，以下のように読み込みます。

from sklearn import datasets

boston = datasets.load_boston()
X = boston.data
y = boston.target

13種類の特徴量からなるXと予測対象となる連続値yの組が，全部で506個含まれています。

最適な分割を探る

探索する特徴量を一つに固定した場合，二乗誤差和が最小になるような分割点は以下のような関数で求められます：

def find_optimal_division(x, y):
    list_criterion = []
    x_unique = np.unique(x)

    if len(x_unique) == 1:
        return None, func_criterion(y)

    for threshold in x_unique[:-1]:

        mask_divide = x > threshold
        y_left = y[mask_divide]
        y_right = y[~mask_divide]

        criterion_divide = criterion_lr(y_left, y_right, func_criterion)
        list_criterion.append(criterion_divide)
        
    array_criterion = np.array(list_criterion)
    i_div_opt = np.argmin(array_criterion)
    
    return x_unique[i_div_opt], array_criterion[i_div_opt]

基本的には分類の時と同じですが，指標の計算部分だけ異なります。

criterion_divide = criterion_lr(y_left, y_right, func_criterion)

ここの部分は，あらかじめ下記のように定義しておきます。

func_criterion = np.var

def criterion_lr(y_left, y_right, func_criterion):
    
    return (func_criterion(y_left) * len(y_left) + func_criterion(y_right) * len(y_right) ) / (len(y_left) + len(y_right))

誤差の二乗和を求めるにあたり，numpyのnp.varを使い分散を求めたうえでデータ数をかけるという方法をとりました。

無駄があるように見えますが，

numpyの最適化されている関数を有効活用できる
別の指標を扱う際にもコードの改変が少ない（後述）

というメリットがあります。

なお，求めた誤差二乗和をデータ数で割ることで，平均二乗誤差(Mean Squared Error, MSE)にしてあります。平均を取らなくても計算上問題ありませんが，こちらも別指標を使う際に共通化できるようにこのようなにしてあります。

これを例によってnumpyのapply_alogn_axis関数に適用することで，全特徴量から最適な分割を探索できます。

def divide(X, y):

    results = np.apply_along_axis(find_optimal_division, 0, X, y)

    arg_div = np.argmin(results[1])
    x_div = results[0, arg_div]
    criterion_opt = results[1, arg_div]

    return arg_div, x_div, criterion_opt

あとは今までと同一です。

使ってみる

以上をクラスにしてまとめ，実行してみました。コードはGitHubにあげてあります。

github.com

boston = datasets.load_boston()
X = boston.data
y = boston.target

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y)

tree = MyTree(criterion='mse', max_depth=5, verbose=True)
tree.fit(X_train, y_train)

y_pred = tree.predict(X_test)
1 - np.sum((y_test - y_pred) ** 2) / len(y_test) / y_test.var()
# 0.7113242620957975

モデルの評価には，決定係数を用いました：

$R^{2} = 1 - \frac{\sum _ {i=1} ^{N} (y_i - y_i^{\rm (pred)})^{2}}{\sum _ {i=1} ^{N} (y_i - y_i^{\rm (mean)})^{2}}$

データ自身のばらつきと比べて，予測値と実際の値のズレが小さいほど，この値は1に近い値となります。

今回作成したモデルでは決定係数は0.71でした。モデルの誤差がデータ自身のばらつきの3割程度に抑えられているので，パラメータチューニングを行なっていないにも関わらずそれなりによい値のように思えます。

まとめ

回帰についても，評価指標の部分を変えるだけで実装することができました。

なお，評価指標の部分で，criterion_lr(y_left, y_right, func_criterion)という関数を導入しました。ここの部分の引数func_criterionに別の関数を入れてやることで，別の指標を使うことができます。

例えば，Gini係数を下記のように定義して代入すると，左右の分割断片断片でのGini係数の平均値（データ数で重み付き）を出すことができます。

def gini(x):
    _, counts = np.unique(x, return_counts=True)
    prob = counts / len(x)
    
    return 1 - (prob * prob).sum()

func_criterion = gini
criterion_divide = criterion_lr(y_left, y_right, func_criterion)

このようにして，異なる指標でも引数に入れる関数を変えるだけで共通のコードを使い回せるようにしてあります。

なお，MyTreeクラスにまとめ直す際に，引数にcriterionを加えました。これをginiもしくはmseにすることで，分類と回帰を切り替えられるようにしてあります。

今回は，決定木による回帰を，以前作成した分類との共通部分を意識しながら実装していきました。

しかしscikit-learnの実装では，分類はDecisionTreeClassifier，回帰はDecisionTreeRegressorと，異なるインターフェイスに分けています。とはいっても，共通する部分は継承をすることで，重複したコードが少なくなるように実装されているはずです。

自己流での実装もひと段落しましたので，次回からはいよいよscikit-learnの中身を覗いていきたいと思います。

*1:だんだん冒頭ネタもいらすとも苦しくなってきた…

*2:Gini係数の他に，情報エントロピー $-\sum_{k=1}^{K}p_k \log p_k$ を使う場合もあります。