どこから見てもメンダコ

軟体動物門頭足綱八腕類メンダコ科

遺伝的アルゴリズムで特徴量選択

はじめに

f:id:horomary:20190310190326p:plain:w300

pythonによる進化計算アルゴリズムフレームワークであるdeapを使用して、機械学習モデリングのための特徴選択を遺伝的アルゴリズムで行ってみます。

github.com

逐次的な特徴量選択アルゴリズムと比較した場合、遺伝的アルゴリズムを使って嬉しいのは ”特徴量の数は30個以内で、最大の精度がでる特徴量セット” を探索する、というように使用する特徴量の数を明示的に制限することができる点です。

また、使用する特徴量の数の制限を厳しくした場合には逐次的な特徴量選択手法よりも精度のよいモデルが得られることが(経験則ですが)多いです。

嬉しくないのはGAの計算コストがかなり重いことです。


ナップザック問題としての特徴量選択

N個の特徴量のそれぞれについて使用するか使用しないかを選択する。 使用する特徴量の数は設定値以下に制限する。

というナップサック問題と捉えることができるのでGAで容易に解決できます。


データセットの準備

みんな大好きボストン住宅データセットを使いますが、 そのままでは特徴量の数が13と寂しいのでPolynomialFeaturesで特徴量を増やします。

import pandas as pd
import numpy as np
from sklearn.datasets import load_boston
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler


# データセットのロード
X = pd.DataFrame(load_boston().data, columns=load_boston().feature_names)
y = load_boston().target

# 多項式特徴量を追加
poly = PolynomialFeatures(2)
poly.fit(X)
X_poly = pd.DataFrame(poly.transform(X), columns=poly.get_feature_names(input_features=X.columns))

# 標準化
StSc = StandardScaler()
X_sc = pd.DataFrame(StSc.fit_transform(X), columns=X.columns)
X_poly_sc = pd.DataFrame(StSc.fit_transform(X_poly), columns=X_poly.columns)

これによって特徴量数105, サンプル数506のデータセットが出来上がりました。

ベンチマーク

特徴量選択をする前に、そのままのデータセット (特徴数13) と 多項式特徴量を追加したデータセット(特徴数105) でどの程度の予測精度が出るのか見ておきましょう。

from sklearn.linear_model import RidgeCV
from sklearn.model_selection import train_test_split

# そのままのデータセット
scores = []
for _ in range(30):
    X_train, X_test, y_train, y_test = train_test_split(X_sc, y, test_size=0.4)
    model = RidgeCV()
    model.fit(X_train, y_train)
    scores.append(model.score(X_test, y_test))

print(np.array(scores).mean()) # 0.70
# 多項式特徴量を追加した場合
scores = []
for _ in range(300):
    X_train, X_test, y_train, y_test = train_test_split(X_poly_sc, y, test_size=0.4)
    model = RidgeCV()
    model.fit(X_train, y_train)
    scores.append(model.score(X_test, y_test))

print(np.array(scores).mean())    # 0.82

[結果]

特徴量の数 13 105
精度 0.70 0.82

多項式特徴量を追加しただけでそこそこ精度が上がりました。
ちなみに多項式特徴量で回帰係数が最大だったのはRM2 (部屋の数の二乗) でした。尤もらしいですね。


GAによる変数選択

いよいよ遺伝的アルゴリズムで変数選択をします。
deapのナップザック問題のサンプルを基にコーディングしました。

deap/knapsack.py at master · DEAP/deap · GitHub

この例では次世代に引き継ぐ個体選択にはNSGA-Ⅱを採用しています。
GAでは思考停止でこれを使っとけばだいたいうまくいきます(※個人の印象です)

from deap import algorithms
from deap import base
from deap import creator
from deap import tools


def evalIndividual(individual):
    """ 個体評価のための関数
        Prameter:
        individual: 要素が0or1のリスト  長さは入力変数の数と同じ
        Return:
        n_features: 使った特徴量の数
        score: モデルの精度
    """
    n_features = sum(individual)

    if n_features == 0:
        return 9999, -9999

    X_temp = X_poly_sc.iloc[:, [bool(val) for val in individual]]
    
    scores = []
    for _ in range(30):
        X_train, X_test, y_train, y_test = train_test_split(X_temp, y, test_size=0.4) 
        model = RidgeCV()
        model.fit(X_train, y_train)
        scores.append(model.score(X_test, y_test))

    score = np.array(scores).mean()

    return n_features, score


creator.create("Fitness", base.Fitness, weights=(-1.0, 1.0))   
creator.create("Individual", list, fitness=creator.Fitness)


toolbox = base.Toolbox()
toolbox.register("attr_bool", random.randint, 0, 1)
toolbox.register("individual", tools.initRepeat, creator.Individual, 
                 toolbox.attr_bool, X_poly_sc.shape[1])
toolbox.register("population", tools.initRepeat, list, toolbox.individual)


toolbox.register("evaluate", evalIndividual)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutFlipBit, indpb=0.05)
toolbox.register("select", tools.selNSGA2)

"""遺伝的アルゴリズム設定
50世代まで
1世代の個体数300
次世代に引き継ぐ個体数100
交叉確率0.7
突然変異確率0.1
"""
NGEN = 50
MU = 100
LAMBDA = 300
CXPB = 0.7
MUTPB = 0.1

pop = toolbox.population(n=MU)
hof = tools.ParetoFront()
stats = tools.Statistics(lambda ind: ind.fitness.values)
stats.register("avg", np.mean, axis=0)
stats.register("std", np.std, axis=0)
stats.register("min", np.min, axis=0)
stats.register("max", np.max, axis=0)

pop, log = algorithms.eaMuPlusLambda(pop, toolbox, MU, LAMBDA, CXPB,
                                     MUTPB, NGEN, stats, halloffame=hof)

30分くらいで終わりました。
パレート解はhofに格納されています。popには最終世代の100個体が格納されています。

※パレート解とは
CAE用語 パレート解:Ansysの導入ならCAE30年のサポート実績:サイバネット

結果

パレートフロンティアは以下のような感じ。
N_featuresが特徴量の数、Scoreが精度です。

f:id:horomary:20190310183819p:plain

特徴量を105から16まで削っても、削減前と同等の精度が出ています。

※選ばれた16の特徴量 f:id:horomary:20190310185730p:plain

パレートフロンティアを可視化。

f:id:horomary:20190310184859p:plain

赤点がパレートフロンティア、青点が最終世代の100個体です。

まとめ

deapは優秀ですね!