[Ubiquant Market Prediction]シンプルなDAE(Denoising Autoencoder)

こんにちは、tmrtjと申します。

私は現在株式会社音圧爆上げくんという会社に所属しており、プロKagglerとして活動しています。


私は現在プロKagglerとして、業務の一環で、現在Kaggleで行われているコンペの一つである[Ubiquant Market Prediction](https://www.kaggle.com/c/ubiquant-market-prediction)に参加しています。



また、Kaggle上にNoteBookとして同じ内容の記事を英語で投稿しています。
リンクは[こちら](https://www.kaggle.com/tmrtj9999/simple-dae-denoising-autoencoder)です。

#なぜDAEを使うのか
一般に金融市場のデータは過学習しやすいことが知られています。そのため金融市場の予測に機械学習を適用する場合には過学習対策に力を注ぐことになります。

例えばKaggleで２０１6年に開催された[金融市場モデリングコンペ](https://www.kaggle.com/c/two-sigma-financial-modeling)では上位入賞者の解法としてExtra Treesという手法が使われています。これはランダムフォレストに似た手法で、ランダムフォレストでは複数の決定木を組み合わせて予測しますが、その代わりにランダムに特徴量を選んだ木を作って、それらのバギングによる予測をする手法です。個々の木における特徴量の選び方がランダムなので、ランダムフォレストと比べてより過学習対策を重くみた手法と言えます。

このように金融データ分析では過学習対策が重要ですが、今回はそのための一つの手法としてDAE(Denosing Autoencoder)を紹介します。これはAutoencoderにノイズを加えて特徴抽出を行う手法で、従来のAutoencoderとは違い、ノイズを乗せているため、よりロバストな特徴量が抽出できることが期待できます。DAEには様々なタイプがありますが、今回はガウスノイズと呼ばれる、取る値が正規分布に従うノイズを入力に足し合わせたDAEを紹介します。



#実装

まずはtrainデータを読み込んでDataFrameにします。
```reduce_mem_usage(df)```　はメモリを節約できるような変換をかける関数です。

```
import os
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

import lightgbm as lgbm
from lightgbm import *

df = pd.read_parquet('../input/ubiquant-parquet/train_low_mem.parquet')

def reduce_mem_usage(df):
  
    start_mem = df.memory_usage().sum() / 1024**2
    print('Memory usage of dataframe is {:.2f} MB'.format(start_mem))
    
    for col in df.columns:
        col_type = df[col].dtype
        
        if col_type != object:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == 'int':
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64)  
            else:
                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)
        else:
            df[col] = df[col].astype('category')

    end_mem = df.memory_usage().sum() / 1024**2
    
    print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))
    print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))
     
    return df

df = reduce_mem_usage(df)
```

<br>

trainデータを特徴量データのX_0, X_1, ... , X_4と、targetデータのY_0, Y_1, ... , Y_4に時系列順に分割します。


```
from sklearn.model_selection import KFold, train_test_split


features = [f'f_{i}' for i in range(300)]
target = 'target'

X_012, X_34, Y_012, Y_34 = train_test_split(df[features], df[target], train_size=0.6, shuffle=False)

X_3, X_4, Y_3, Y_4 = train_test_split(X_34, Y_34, train_size=0.5, shuffle=False)

X_0, X_12, Y_0, Y_12 = train_test_split(X_012, Y_012, train_size=0.33, shuffle=False)

X_1, X_2, Y_1, Y_2 = train_test_split(X_12, Y_12, train_size=0.5, shuffle=False)

df = [[]]

X_0 = X_0.reset_index()
X_1 = X_1.reset_index()
X_2 = X_2.reset_index()
X_3 = X_3.reset_index()
X_4 = X_4.reset_index()

Y_0 = Y_0.reset_index()
Y_1 = Y_1.reset_index()
Y_2 = Y_2.reset_index()
Y_3 = Y_3.reset_index()
Y_4 = Y_4.reset_index()

del X_0['index']
del X_1['index']
del X_2['index']
del X_3['index']
del X_4['index']
del Y_0['index']
del Y_1['index']
del Y_2['index']
del Y_3['index']
del Y_4['index']
```
<br>
<br>
まずはDAEを訓練します。X_0, Y_0 で訓練します。validation dataとしてはX_1, Y_1を使用します。

```
import keras
from keras import layers
from keras.layers import Input, Dense, BatchNormalization
from keras.layers.noise import GaussianNoise
from keras.models import Model, load_model
from keras.callbacks import EarlyStopping, TensorBoard, ModelCheckpoint

# EaelyStopping
early_stopping =  EarlyStopping(
                            monitor='val_loss',
                            min_delta=0.0,
                            patience=3,
)

input_dim = X_0.shape[1]
encoding_dim = 100

input_layer = Input(shape=(input_dim, ))
encoder = GaussianNoise(0.1)(input_layer)
encoder = Dense(encoding_dim, activation="tanh")(encoder)
encoder = BatchNormalization()(encoder)
encoder = Dense(10, activation="relu")(encoder)

decoder = Dense(10, activation='tanh')(encoder)
decoder = Dense(encoding_dim, activation='relu')(decoder)
output_layer = Dense(1, activation='tanh')(decoder)

autoencoder = Model(inputs=input_layer, outputs=output_layer)

autoencoder.compile(optimizer="adam", loss='mean_squared_error', metrics=["mse"])

modelCheckpoint = ModelCheckpoint(filepath = 'XXX.h5',
                                  monitor='val_loss',
                                  verbose=1,
                                  save_best_only=True,
                                  save_weights_only=False,
                                  mode='min',
                                  save_freq=1)


autoencoder.fit(X_0, Y_0, batch_size=1440, epochs=11890,
          validation_data=(X_1, Y_1),
          callbacks=[early_stopping] # CallBacks
            )
```
<br>

次に学習が終わったら、訓練したDAEをX_2~X_4に適用して、特徴量を作ります。

```
X_0 = [[]]
X_1 = [[]]
Y_0 = []
Y_1 = []
df = [[]]


dae= Model(input_layer, encoder)

encoding_X_2 = dae.predict(X_2)
dae_X_2 = pd.DataFrame(encoding_X_2)
encoding_X_3 = dae.predict(X_3)
dae_X_3 = pd.DataFrame(encoding_X_3)
encoding_X_4 = dae.predict(X_4)
dae_X_4 = pd.DataFrame(encoding_X_4)

X_2 = pd.concat([X_2, dae_X_2], axis=1)
X_3 = pd.concat([X_3, dae_X_3], axis=1)
X_4 = pd.concat([X_4, dae_X_4], axis=1)
```
例えば、X_2はDAE使用後はこのようになっています。

```
print(X_2.head(4))
```

**結果**

```
f_0       f_1       f_2       f_3       f_4       f_5       f_6  \
0 -0.930176 -0.017868  0.228271 -0.281250 -0.250000 -2.283203 -2.091797   
1  0.246582 -0.548340  0.625000 -0.350342 -0.395508  1.147461  0.520020   
2  0.943848 -0.902344  1.949219 -0.340332 -0.495850  0.558594 -0.353271   
3  0.645508  0.689453  1.122070 -0.165283 -0.421387  1.242188  0.900879   

        f_7       f_8       f_9  ...         0         1         2         3  \
0 -0.971680  0.202271 -1.843750  ...  0.000000  0.000000  0.000000  1.004568   
1 -0.496826  0.202271 -1.804688  ...  1.464177  0.143047  2.361402  0.000000   
2  1.024414  0.202271 -0.614746  ...  1.254612  0.807103  2.559844  0.000000   
3 -0.653320  0.202271  0.752441  ...  0.686005  0.000000  0.090134  1.171950   

          4         5         6         7         8         9  
0  0.000000  0.000000  0.000000  0.000000  4.065059  0.000000  
1  0.000000  0.112643  1.115756  0.000000  0.160429  0.000000  
2  2.319442  0.142256  2.297849  0.226917  0.000000  0.148446  
3  0.615372  0.000000  2.083543  0.000000  0.462596  0.054800  

[4 rows x 310 columns]
```

これらの'0'~'9'がDAEによって新たに増えた特徴量です。
<br>
<br>








次にLightGBMを使ってみます。
trainデータにはX_2, Y_2を使い、valデータにはX_3, Y_3を、testデータにはX_4, Y_4を使います。


```
import warnings
import numpy as np
import lightgbm as lgb
from scipy.stats import pearsonr

warnings.simplefilter('ignore')

lgb_train = lgb.Dataset(X_2, Y_2)
lgb_eval = lgb.Dataset(X_3, Y_3, reference=lgb_train)

params = {'seed': 1,
          'verbose' : -1,
           'objective': "regression",
           'learning_rate': 0.02,
           'bagging_fraction': 0.2,
           'bagging_freq': 1,
           'feature_fraction': 0.3,
           'max_depth': 5,
           'min_child_samples': 50,
           'num_leaves': 64}
        
        
gbm = lgb.train(params,
                lgb_train,
                num_boost_round=1000,
                valid_sets=lgb_eval,
                #verbose_eval=False,
                early_stopping_rounds=10,
                )


Y_pred = gbm.predict(X_4, num_iteration=gbm.best_iteration)
l_4 = Y_4['target'].values.tolist()
score_tuple = pearsonr(l_4, Y_pred)
score = score_tuple[0]
print(f"Validation Pearsonr score : {score:.4f}")
```

<br>
LightGBMの学習が終わったら特徴量重要度を表示させます

```
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report


feature = gbm.feature_importance(importance_type='gain')


f = pd.DataFrame({'number': range(0, len(feature)),
             'feature': feature[:]})
f2 = f.sort_values('feature',ascending=False)

#features' name
label = X_2.columns[0:]

#feature rank
indices = np.argsort(feature)[::-1]

for i in range(len(feature)):
    print(str(i + 1) + "   " + str(label[indices[i]]) + "   " + str(feature[indices[i]]))
```

<br>
<br>

**結果**

```
1   5   7196.480256080627
2   3   4111.950593948364
3   f_231   3353.084394454956
4   0   3047.538405418396
5   f_119   2898.0525064468384
6   9   2198.6642055511475
7   8   1893.9976024627686
8   f_74   1485.8169984817505
9   f_276   1391.6077098846436
10   f_251   1381.2403945922852

// ... 略
```
これらの'0'~'9'がDAEによって新たに増えた特徴量です。
特徴量重要度を確認すると、確かにDAEによって増やした特徴量が効いていることがわかります。
<br>

#さいごに
株式会社音圧爆上げくんではプロKagglerを募集しています。
ご興味のある方はぜひ以下のリンクをご覧ください。
[Wantedlyリンク](https://www.wantedly.com/projects/748285)