GPでベットをするかしないかを予測する関数を作ります。
前回は過学習してしまい、テストデータでの利益を上げることができませんでした。
そこで、今回はすこし解候補の選択方法を変更して、過学習を抑えようと思います。
変更する選択方法をε-Lexicase Selectionというのですが、すごく単純に言うと、精度がめちゃめちゃ上がる手法です。neural networkにおけるbatch normalizationみたいな立ち位置だと思ってください。
この話を解説するとめちゃめちゃ長いので解説はしません。論文を読んでください。
Epsilon-Lexicase Selection for Regression | Proceedings of the Genetic and Evolutionary Computation Conference 2016
dl.acm.org
学習過程はipynbファイルで行いました。以下はその記録です。
前処理
まず、使いそうなライブラリをインポートします。
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display, Math学習データの読み込み
df = pd.read_pickle('/work/formatted_source_data/analysis_data02.pkl')次に学習データとテストデータ、バリデーションデータの切り分け、および入力データと教師データの切り分けをします。データの正規化もします。
教師データは単勝が当たったかどうか(tansyo_hit)と単勝が当たった時の支払い金額(tansyo_payout)を使います。
df = df.sort_values(['race_date']).reset_index(drop=True)
not_use_columns = [
'race_date', 'race_id', 'race_grade', 'name', 'jocky_name', 'odds', 'popular', 'rank', 'time', 'prize', 'tansyo_hit', 'tansyo_payout', 'hukusyo_hit', 'hukusyo_payout',
]
not_use_df = df[not_use_columns]
target_df = df.drop(columns=not_use_columns)
train_df = target_df.loc[not_use_df['race_date'] < '2018-01-01', :].reset_index(drop=True)
valid_df = target_df.loc[('2018-01-01' <= not_use_df['race_date']) & (not_use_df['race_date'] < '2022-01-01'), :].reset_index(drop=True)
test_df = target_df.loc['2022-01-01' <= not_use_df['race_date'], :].reset_index(drop=True)
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
train_X = pd.DataFrame(data=scaler.fit_transform(train_df), columns=train_df.columns)
valid_X = pd.DataFrame(data=scaler.transform(valid_df), columns=train_df.columns)
test_X = pd.DataFrame(data=scaler.transform(test_df), columns=train_df.columns)
train_X = train_X.fillna(0)
valid_X = valid_X.fillna(0)
test_X = test_X.fillna(0)
Y_columns = ['tansyo_hit', 'tansyo_payout']
train_Y = not_use_df.loc[not_use_df['race_date'] < '2018-01-01', Y_columns].reset_index(drop=True)
valid_Y = not_use_df.loc[('2018-01-01' <= not_use_df['race_date']) & (not_use_df['race_date'] < '2022-01-01'), Y_columns].reset_index(drop=True)
test_Y = not_use_df.loc['2022-01-01' <= not_use_df['race_date'], Y_columns].reset_index(drop=True)
train_Y = train_Y.rename(columns={'tansyo_hit':'hit', 'tansyo_payout':'payout'})
valid_Y = valid_Y.rename(columns={'tansyo_hit':'hit', 'tansyo_payout':'payout'})
test_Y = test_Y.rename(columns={'tansyo_hit':'hit', 'tansyo_payout':'payout'})学習
さて、ここからDeapを利用するコードを書いていきます。GPのサンプルコードは公式のGithub内にあるので、それらを参考にします。
まず、使うライブラリをインポートします。
import operator
from deap import base
from deap import creator
from deap import tools
from deap import gp
import random
from functools import partial割り算の演算関数を再定義します。pythonの素の割り算やnumpyの割り算は、0割するとエラーやwarnningがでてきますし、nanで埋められてしまいます。その値が他の関数の入力に来ると、nanがどんどん増殖してしまうため、nanを1として再定義します。
def protectedDiv(left, right):
with np.errstate(divide='ignore',invalid='ignore'):
x = np.divide(left, right)
if isinstance(x, np.ndarray):
x[np.isinf(x)] = 1
x[np.isnan(x)] = 1
elif np.isinf(x) or np.isnan(x):
x = 1
return x次にGPでつかう関数群を定義します。これはdeap.gp.PrimitiveSetで定義されます。
PrimitiveSetの第2引数は、この関数の入力数を表すので、特徴量数を指定します。
pset = gp.PrimitiveSet("MAIN", len(train_X.columns))
pset.addPrimitive(np.add, 2, name="add")
pset.addPrimitive(np.subtract, 2, name="sub")
pset.addPrimitive(np.multiply, 2, name="mul")
pset.addPrimitive(protectedDiv, 2, name='div')
pset.addPrimitive(np.negative, 1, name="neg")
pset.addPrimitive(np.cos, 1, name="cos")
pset.addPrimitive(np.sin, 1, name="sin")
pset.addEphemeralConstant("rand101", partial(random.randint, -1, 1))
pset.renameArguments(**{f'ARG{i}':f'x_{col}' for i, col in enumerate(train_X.columns)})オリジナルのlexicase selectionではFitness case num = 学習データ数ですが、このアルゴリズムは学習データ数が多くなると計算時間がとても増えます。
今回のデータ数は50万件以上なため、普通のPCだと計算が終わりません。
そこで、今回は学習データをいくつかのバッチで切り分け、その分けたバッチ数をfittness case numとします。
FITNESS_CASE_NUM = 100
creator.create("FitnessMax", base.Fitness, weights=[1.0]*FITNESS_CASE_NUM)
creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMax)
toolbox = base.Toolbox()
toolbox.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=2)
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("compile", gp.compile, pset=pset)評価関数を定義します。個体(解候補)が与えられたときに、その評価値を返せばよいです。出力値が0以上のときに賭け、そうでないときは賭けないようにして、利益を算出します。
その総利益を評価値として返します。
また、FITNESS_CASE_NUM個に学習データを切り分け、それぞれの評価値を計算し返します。
train_inputs = {f'x_{col}':train_X.loc[:, col].values for col in train_X.columns}
def eval_individual(individual):
weights = np.array(individual)
ret = neural_network(weights, train_X)
ret = ret[train_Y_mask]
bet = np.zeros_like(ret)
bet[ret > 0.5] = 1
win_money = bet*(train_hit_Y_flat * train_payout_Y_flat - 1.0)
fitnesses = []
for sp_win_money in np.array_split(win_money, FITNESS_CASE_NUM):
fitnesses.append(np.sum(sp_win_money))
return fitnesses定義した評価関数や、deapに内蔵されている選択、交叉、突然変異の関数をtoolboxに登録します。また、個体の深さ(合成関数の処理規模だと思ってください)に制限を設けます。
選択方法をtools.selAutomaticEpsilonLexicaseに設定します。
toolbox.register("evaluate", eval_individual)
toolbox.register("select", tools.selAutomaticEpsilonLexicase)
toolbox.register("mate", gp.cxOnePoint)
toolbox.register("expr_mut", gp.genFull, min_=0, max_=2)
toolbox.register('mutate', gp.mutUniform, expr=toolbox.expr_mut, pset=pset)
toolbox.decorate("mate", gp.staticLimit(key=operator.attrgetter("height"), max_value=12))
toolbox.decorate("mutate", gp.staticLimit(key=operator.attrgetter("height"), max_value=12))準備が揃ったので、学習を始めます。
個体数500、世代数200で学習します。
crossover_rate = 0.8
mutate_rate = 0.2
max_generation = 200
population_size = 500
random.seed(318)
population = toolbox.population(n=population_size)
halloffame = tools.HallOfFame(1)
stats_fit = tools.Statistics(lambda ind: np.sum(ind.fitness.values))
height_size = tools.Statistics(lambda ind: ind.height)
mstats = tools.MultiStatistics(fitness=stats_fit, height=height_size)
mstats.register("avg", np.mean)
mstats.register("std", np.std)
mstats.register("min", np.min)
mstats.register("max", np.max)
logbook = tools.Logbook()
logbook.header = ['gen', 'nevals'] + mstats.fields
# Evaluate the individuals with an invalid fitness
invalid_ind = [ind for ind in population if not ind.fitness.valid]
fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
for ind, fit in zip(invalid_ind, fitnesses):
ind.fitness.values = fit
if halloffame is not None:
halloffame.update(population)
record = mstats.compile(population)
logbook.record(gen=0, nevals=len(invalid_ind), **record)
print(logbook.stream)
elite = None
# Begin the generational process
for gen in range(1, max_generation + 1):
# Select the next generation individuals
best = toolbox.clone(max(population, key=lambda ind:np.sum(ind.fitness.values)))
if elite == None or np.sum(elite.fitness.values) < np.sum(best.fitness.values):
elite = best
offspring = toolbox.select(population, len(population)-1)
offspring.append(toolbox.clone(elite))
offspring = [toolbox.clone(ind) for ind in population]
# Apply crossover and mutation on the offspring
for i in range(1, len(offspring), 2):
if random.random() < crossover_rate:
offspring[i - 1], offspring[i] = toolbox.mate(offspring[i - 1], offspring[i])
del offspring[i - 1].fitness.values, offspring[i].fitness.values
for i in range(len(offspring)):
if random.random() < mutate_rate:
offspring[i], = toolbox.mutate(offspring[i])
del offspring[i].fitness.values
# Evaluate the individuals with an invalid fitness
invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
for ind, fit in zip(invalid_ind, fitnesses):
ind.fitness.values = fit
# Update the hall of fame with the generated individuals
if halloffame is not None:
halloffame.update(offspring)
# Replace the current population by the offspring
population[:] = offspring
# Append the current generation statistics to the logbook
record = mstats.compile(population)
logbook.record(gen=gen, nevals=len(invalid_ind), **record)
print(logbook.stream)実行するとログが流れます。終わるまで待ちます。
fitness height
------------------------------------------------------- ------------------------------------------------
gen nevals avg gen max min nevals std avg gen max min nevals std
...
198 430 -74512.8 198 0 -138788 430 29197.8 3.312 198 12 0 430 2.2268
199 416 -74042.3 199 0 -138788 416 28761.3 3.33 199 12 0 416 2.22915
200 413 -73816.3 200 0 -138788 413 29875 3.35 200 12 0 413 2.27057 一番良かった解候補の利益を算出してみます。
best_ind = toolbox.clone(max(population+[elite], key=lambda ind:np.sum(ind.fitness.values)))
def calc_win_money(X, Y):
inputs = {f'x_{col}':X.loc[:, col].values for col in X.columns}
func = toolbox.compile(expr=best_ind)
ret = func(**inputs)
bet = np.zeros_like(ret)
bet[ret > 0] = 1
win_money = bet*(Y['hit']*Y['payout'] - 1.0)
return np.sum(win_money)
print(f'train win money: {calc_win_money(train_X, train_Y)}')
print(f'train win money: {calc_win_money(valid_X, valid_Y)}')
print(f'train win money: {calc_win_money(test_X, test_Y)}')学習データは利益がでましたが、バリデーション、テストデータは損失がでました。あまりうまくいかないようです。
前回の素のGPよりも良さそうですが、結局テストデータの利益はでませんでした。
train win money: 858.4000000000007
train win money: -1394.6999999999998
train win money: -1892.3