kaggleの雑記

kaggleのお手本のような時系列特徴量生成【Ion Switching】

コンペについて

https://www.kaggle.com/c/liverpool-ion-switching

リバプールの大学が行っている細胞研究です。観測器から取得した電気信号のデータを元にして細胞のチャンネルが開放されている瞬間を検知することが目標です。

今回読み込むカーネル

Rob MullaさんのIon Switching – 5kfold LGBM & Trackingを読み込もうと思います。

データの読み込み

train = pd.read_csv('../input/liverpool-ion-switching/train.csv')
test = pd.read_csv('../input/liverpool-ion-switching/test.csv')
ss = pd.read_csv('../input/liverpool-ion-switching/sample_submission.csv')
train['train'] = True
test['train'] = False
tt = pd.concat([train, test], sort=False).reset_index(drop=True)
tt['train'] = tt['train'].astype('bool')

このカーネルではこのタイミングでttと言う一つのDataFrameにまとめています。trainデータにはTrue、テストデータにはFalseをつけて最後の予測のタイミングで分割する仕組みです。

ハイパーパラメータの定義

# ここはちょっとわかりませんでした。
MODEL = os.path.basename(__file__).split('.')[0]

TARGET = 'open_channels'

# LightGBMのパラメータです
TOTAL_FOLDS = 5
RANDOM_SEED = 529
MODEL_TYPE = 'LGBM'
LEARNING_RATE = 0.009
SHUFFLE = True
NUM_BOOST_ROUND = 500_000
EARLY_STOPPING_ROUNDS = 50
N_THREADS = -1
OBJECTIVE = 'regression'
METRIC = 'rmse'
NUM_LEAVES = 2**8+1
MAX_DEPTH = -1
FEATURE_FRACTION = 1
BAGGING_FRACTION = 1
BAGGING_FREQ = 0

トラッキングとして記録

###########
# TRACKING
###########

run_id = "{:%m%d_%H%M}".format(datetime.now())

def update_tracking(
        run_id, field, value, csv_file="./tracking.csv",
        integer=False, digits=None, nround=6,
        drop_broken_runs=False):
    """
    Tracking function for keep track of model parameters and
    CV scores. `integer` forces the value to be an int.
    """
    try:
        df = pd.read_csv(csv_file, index_col=[0])
    except:
        df = pd.DataFrame()
    if drop_broken_runs:
        df = df.dropna(subset=['1_f1'])
    if integer:
        value = round(value)
    elif digits is not None:
        value = round(value, digits)
    df.loc[run_id, field] = value  # Model number is index
    df = df.round(nround)
    df.to_csv(csv_file)

# Update Tracking
update_tracking(run_id, 'model_number', MODEL)
update_tracking(run_id, 'model_type', MODEL_TYPE)
update_tracking(run_id, 'seed', RANDOM_SEED, integer=True)
update_tracking(run_id, 'nfolds', TOTAL_FOLDS, integer=True)
update_tracking(run_id, 'lr', LEARNING_RATE)
update_tracking(run_id, 'shuffle', SHUFFLE)
update_tracking(run_id, 'boost_rounds', NUM_BOOST_ROUND)
update_tracking(run_id, 'es_rounds', EARLY_STOPPING_ROUNDS)
update_tracking(run_id, 'threads', N_THREADS)
update_tracking(run_id, 'objective', OBJECTIVE)
update_tracking(run_id, 'metric', METRIC)
update_tracking(run_id, 'num_leaves', NUM_LEAVES)
update_tracking(run_id, 'max_depth', MAX_DEPTH)
update_tracking(run_id, 'feature_fraction', FEATURE_FRACTION)
update_tracking(run_id, 'bagging_fraction', BAGGING_FRACTION)
update_tracking(run_id, 'bagging_freq', BAGGING_FREQ)

ここは自分のトライアル結果を記録するため、それをtracking.csvというファイルに残すためのコードです。

画像の様な1行のcsvが作成され、記録が残ります。あとからどんなパラメータでスコアを出せたのか、確認しやすい状態になっています。

バッチサイズに分割する

# # Include batch

# 時間でsortする。reset_indexは上でやっているので、いらないかも
tt = tt.sort_values(by=['time']).reset_index(drop=True)
# indexを振る
tt.index = ((tt.time * 10_000) - 1).values
# バッチサイズ50000で区切る
tt['batch'] = tt.index // 50_000
# バッチないのindexを作成する
tt['batch_index'] = tt.index - (tt.batch * 50_000)
# バッチスライスという形で、バッチをさらに5000で区切る
tt['batch_slices'] = tt['batch_index'] // 5_000
# zfillで右詰で左側に0がはいる。バッチとバッチスライスを文字列連結している
tt['batch_slices2'] = tt.apply(lambda r: '_'.join(
    [str(r['batch']).zfill(3), str(r['batch_slices']).zfill(3)]), axis=1)

50000のバッチサイズ毎に特徴量を生成する

# 50_000 Batch Features
# 最小値
tt['signal_batch_min'] = tt.groupby('batch')['signal'].transform('min')
# 最大値
tt['signal_batch_max'] = tt.groupby('batch')['signal'].transform('max')
# 標準偏差
tt['signal_batch_std'] = tt.groupby('batch')['signal'].transform('std')
# 平均
tt['signal_batch_mean'] = tt.groupby('batch')['signal'].transform('mean')
# 前回との差分の平均
tt['mean_abs_chg_batch'] = tt.groupby(['batch'])['signal'].transform(
    lambda x: np.mean(np.abs(np.diff(x))))
# 絶対値の最大値
tt['abs_max_batch'] = tt.groupby(
    ['batch'])['signal'].transform(lambda x: np.max(np.abs(x)))
# 絶対値の最初うち
tt['abs_min_batch'] = tt.groupby(
    ['batch'])['signal'].transform(lambda x: np.min(np.abs(x)))

# 最大値と最小値のギャップ
tt['range_batch'] = tt['signal_batch_max'] - tt['signal_batch_min']
# 最大値÷最小値
tt['maxtomin_batch'] = tt['signal_batch_max'] / tt['signal_batch_min']
# 最大値(絶対値)と最小値(絶対値)の平均
tt['abs_avg_batch'] = (tt['abs_min_batch'] + tt['abs_max_batch']) / 2

5000のバッチサイズ毎に特徴量を生成する

# 5_000 Batch Features
# 5000行から5000行に変わっただけで、生成している特徴量は同じ
tt['signal_batch_5k_min'] = tt.groupby(
    'batch_slices2')['signal'].transform('min')
tt['signal_batch_5k_max'] = tt.groupby(
    'batch_slices2')['signal'].transform('max')
tt['signal_batch_5k_std'] = tt.groupby(
    'batch_slices2')['signal'].transform('std')
tt['signal_batch_5k_mean'] = tt.groupby(
    'batch_slices2')['signal'].transform('mean')
tt['mean_abs_chg_batch_5k'] = tt.groupby(['batch_slices2'])[
    'signal'].transform(lambda x: np.mean(np.abs(np.diff(x))))
tt['abs_max_batch_5k'] = tt.groupby(['batch_slices2'])[
    'signal'].transform(lambda x: np.max(np.abs(x)))
tt['abs_min_batch_5k'] = tt.groupby(['batch_slices2'])[
    'signal'].transform(lambda x: np.min(np.abs(x)))

tt['range_batch_5k'] = tt['signal_batch_5k_max'] - tt['signal_batch_5k_min']
tt['maxtomin_batch_5k'] = tt['signal_batch_5k_max'] / tt['signal_batch_5k_min']
tt['abs_avg_batch_5k'] = (tt['abs_min_batch_5k'] + tt['abs_max_batch_5k']) / 2

シフトによるラグ特徴量

# add shifts
tt['signal_shift+1'] = tt.groupby(['batch']).shift(1)['signal']
tt['signal_shift-1'] = tt.groupby(['batch']).shift(-1)['signal']
tt['signal_shift+2'] = tt.groupby(['batch']).shift(2)['signal']
tt['signal_shift-2'] = tt.groupby(['batch']).shift(-2)['signal']

前後2つ先までのsignalを特徴量として取り入れています。

testデータとくっつけて処理しているため、trainに2行だけデータが入り込んでしまうけれども、誤差ということだろうか?そこに関しては不明です。

生成した特徴とその時の目的関数との差分?

for c in ['signal_batch_min', 'signal_batch_max',
          'signal_batch_std', 'signal_batch_mean',
          'mean_abs_chg_batch', 'abs_max_batch',
          'abs_min_batch',
          'range_batch', 'maxtomin_batch', 'abs_avg_batch',
          'signal_shift+1', 'signal_shift-1',
          'signal_batch_5k_min', 'signal_batch_5k_max',
          'signal_batch_5k_std',
          'signal_batch_5k_mean', 'mean_abs_chg_batch_5k',
          'abs_max_batch_5k', 'abs_min_batch_5k',
          'range_batch_5k', 'maxtomin_batch_5k',
          'abs_avg_batch_5k','signal_shift+2','signal_shift-2']:
    tt[f'{c}_msignal'] = tt[c] - tt['signal']

新しく生成した特徴量とsignalの差分を取得しています。

実際に学習に用いる特徴量を選別

FEATURES = ['signal',
            'signal_batch_min',
            'signal_batch_max',
            'signal_batch_std',
            'signal_batch_mean',
            'mean_abs_chg_batch',
            #'abs_max_batch',
            #'abs_min_batch',
            #'abs_avg_batch',
            'range_batch',
            'maxtomin_batch',
            'signal_batch_5k_min',
            'signal_batch_5k_max',
            'signal_batch_5k_std',
            'signal_batch_5k_mean',
            'mean_abs_chg_batch_5k',
            'abs_max_batch_5k',
            'abs_min_batch_5k',
            'range_batch_5k',
            'maxtomin_batch_5k',
            'abs_avg_batch_5k',
            'signal_shift+1',
            'signal_shift-1',
            # 'signal_batch_min_msignal',
            'signal_batch_max_msignal',
            'signal_batch_std_msignal',
            # 'signal_batch_mean_msignal',
            'mean_abs_chg_batch_msignal',
            'abs_max_batch_msignal',
            'abs_min_batch_msignal',
            'range_batch_msignal',
            'maxtomin_batch_msignal',
            'abs_avg_batch_msignal',
            'signal_shift+1_msignal',
            'signal_shift-1_msignal',
            'signal_batch_5k_min_msignal',
            'signal_batch_5k_max_msignal',
            'signal_batch_5k_std_msignal',
            'signal_batch_5k_mean_msignal',
            'mean_abs_chg_batch_5k_msignal',
            'abs_max_batch_5k_msignal',
            'abs_min_batch_5k_msignal',
            #'range_batch_5k_msignal',
            'maxtomin_batch_5k_msignal',
            'abs_avg_batch_5k_msignal',
            'signal_shift+2',
            'signal_shift-2']

print('....: FEATURE LIST :....')
print([f for f in FEATURES])

update_tracking(run_id, 'n_features', len(FEATURES), integer=True)
update_tracking(run_id, 'target', TARGET)

実際に学習に使う特徴量をここで明示しています。コメントアウトされているのはカーネルの作成者がチューニングした後だと思われます。

この様な形で学習に利用する特徴量を変更しています。

オリジナルメトリックの関数

def lgb_Metric(preds, dtrain):
    labels = dtrain.get_label()
    print(preds.shape)
    print(preds)
    preds = np.argmax(preds, axis=0)
    # score = metrics.cohen_kappa_score(labels, preds, weights = 'quadratic')
    score = f1_score(labels, preds, average='macro')
    return ('KaggleMetric', score, True)

sklearnのf1_scoreを利用して評価を行っています。F1スコアの定義は以下の様になっています。

https://en.wikipedia.org/wiki/F1_score

学習データを整頓する

# あらためて、trainデータか否かをbooleanに直す
tt['train'] = tt['train'].astype('bool')
# train = Trueのデータを変数trainに
train = tt.query('train').copy()
# train = Falseのデータを変数testに
test = tt.query('not train').copy()
# open_channelsをint型に指定
train['open_channels'] = train['open_channels'].astype(int)

#学習に使う特徴量だけを抜粋
X = train[FEATURES]
X_test = test[FEATURES]
y = train[TARGET].values

submission用にtimeだけのDataFrameを作成
sub = test[['time']].copy()
groups = train['batch']

マルチクラスと回帰を選択する

if OBJECTIVE == 'multiclass':
    NUM_CLASS = 11
else:
    NUM_CLASS = 1

多分やっていることはそのはずです。

回帰(regression)の時はNUM_CLASSを指定しては行けないと思い込んでいましたが、1を指定しておけばちゃんと動作する様ですね。

回帰の時にはparamsのNUM_CLASSは1にしておけば十分動作する

params

# define hyperparammeter (some random hyperparammeters)
params = {'learning_rate': LEARNING_RATE,
          'max_depth': MAX_DEPTH,
          'num_leaves': NUM_LEAVES,
          'feature_fraction': FEATURE_FRACTION,
          'bagging_fraction': BAGGING_FRACTION,
          'bagging_freq': BAGGING_FREQ,
          'n_jobs': N_THREADS,
          'seed': RANDOM_SEED,
          'metric': METRIC,
          'objective': OBJECTIVE,
          'num_class': NUM_CLASS
          }

Kfoldの指定・trainのsignalとopen_channelsだけのDataFrameの作成・特徴量を列にしたDataFrameの作成

kfold = KFold(n_splits=TOTAL_FOLDS, shuffle=SHUFFLE, random_state=RANDOM_SEED)

oof_df = train[['signal', 'open_channels']].copy()
fi_df = pd.DataFrame(index=FEATURES)

学習と記録

fold = 1
# groupsでバッチ構成を指定して、バッチ毎に学習させていく
for tr_idx, val_idx in kfold.split(X, y, groups=groups):
    # foldは別に変数として持っておいて、標準出力していく
    print(f'====== Fold {fold:0.0f} of {TOTAL_FOLDS} ======')
    
    抜粋したバッチを検証データとして用いて分割する
    X_tr, X_val = X.iloc[tr_idx], X.iloc[val_idx]
    y_tr, y_val = y[tr_idx], y[val_idx]
    # LightGBMに用にlgb.Datasetオブジェクトに変換
    train_set = lgb.Dataset(X_tr, y_tr)
    val_set = lgb.Dataset(X_val, y_val)

    # 最初に決めたハイパーパラメータの元学習
    model = lgb.train(params,
                      train_set,
                      num_boost_round=NUM_BOOST_ROUND,
                      early_stopping_rounds=EARLY_STOPPING_ROUNDS,
                      valid_sets=[train_set, val_set],
                      verbose_eval=50)
    # メトリックはいつでも使える様だがコメントアウトされていた
    # feval=lgb_Metric)

    # 多分類(マルチクラス)なら、そのままpredictionを使用できる
    if OBJECTIVE == 'multi_class':
        preds = model.predict(X_val, num_iteration=model.best_iteration)
        preds = np.argmax(preds, axis=1)
        test_preds = model.predict(X_test, num_iteration=model.best_iteration)
        test_preds = np.argmax(test_preds, axis=1)
    # ウォーターサーバー として都いたのであれば、np.clipで、0〜10の間の整数値に整形が必要
    elif OBJECTIVE == 'regression':
        preds = model.predict(X_val, num_iteration=model.best_iteration)
        preds = np.round(np.clip(preds, 0, 10)).astype(int)
        test_preds = model.predict(X_test, num_iteration=model.best_iteration)
        test_preds = np.round(np.clip(test_preds, 0, 10)).astype(int)
    
    # trainデータの一部を検証データとして利用して、どういう予測がされたのかを格納
    oof_df.loc[oof_df.iloc[val_idx].index, 'oof'] = preds
    # 予測結果を格納する、その時fold番号をつけることで、何fold目の予測なのか明示する
    sub[f'open_channels_fold{fold}'] = test_preds

    # oof_dfを利用して、実績と予測に対してスコアリングする
    # F1 score
    f1 = f1_score(oof_df.loc[oof_df.iloc[val_idx].index]['open_channels'],
                  oof_df.loc[oof_df.iloc[val_idx].index]['oof'],
                  average='macro')
    # rmseスコア
    rmse = np.sqrt(mean_squared_error(oof_df.loc[oof_df.index.isin(val_idx)]['open_channels'],
                                      oof_df.loc[oof_df.index.isin(val_idx)]['oof']))

    # tracking.csvに結果を記録する
    update_tracking(run_id, f'{fold}_best_iter', model.best_iteration, integer=True)
    update_tracking(run_id, f'{fold}_rmse', rmse)
    update_tracking(run_id, f'{fold}_f1', f1)
    # LightGBMの出力するfeature importance(特徴量の重要度)を記録する
    fi_df[f'importance_{fold}'] = model.feature_importance()
    print(f'Fold {fold} - validation f1: {f1:0.5f}')
    print(f'Fold {fold} - validation rmse: {rmse:0.5f}')

    fold += 1

kFoldでバッチサイズ毎に交差検証を実施する。

その結果を各検証毎にF1 ScoreとRMSEの2種類の評価基準で評価し、結果を出力及びcsvに格納する。

最後に検証仕切った後のtrainデータの検証スコアを記録する

of_f1 = f1_score(oof_df['open_channels'],
                  oof_df['oof'],
                  average='macro')
oof_rmse = np.sqrt(mean_squared_error(oof_df['open_channels'],
                                      oof_df['oof']))

oofには前バッチ区間で検証を終えた検証データに対する予測データが格納されているため、そのスコアリングおよび、tracking.csvへの格納を行う。

提出データの作成・及び記録、特徴量の影響度の記録

# DataFrame: sub にopen_channelsと言う列名が含まれる全てのカラム
s_cols = [s for s in sub.columns if 'open_channels' in s]

# open_channelsに名前が含まれる列の中央値を取得
sub['open_channels'] = sub[s_cols].median(axis=1).astype(int)
# 各検証毎のデータを残した状態で保存
sub.to_csv(f'./pred_{MODEL}_{oof_f1:0.6}.csv', index=False)
# 提出用に、'time', 'open_channels'だけの形にして、保存
sub[['time', 'open_channels']].to_csv(f'./sub_{MODEL}_{oof_f1:0.10f}.csv',
                                      index=False,
                                      float_format='%0.4f')

# oof trainの目的変数と交差検証時に予測した値が含まれているdfの保存
oof_df.to_csv(f'./oof_{MODEL}_{oof_f1:0.6}.csv', index=False)

# feature importanceを各検証毎に記録していたので、それの合計値をimportance列にまとめる
fi_df['importance'] = fi_df.sum(axis=1)
fi_df.to_csv(f'./fi_{MODEL}_{oof_f1:0.6}.csv', index=True)

# plotする
fig, ax = plt.subplots(figsize=(15, 30))
fi_df.sort_values('importance')['importance'] \
    .plot(kind='barh',
          figsize=(15, 30),
          title=f'{MODEL} - Feature Importance',
          ax=ax)
plt.savefig(f'./{MODEL}__{oof_f1:0.6}.png')

最後に

ラグ特徴量の取り扱いに、バッチサイズ毎の学習に関してお手本の様なコードを公開してくださったRob Mullaさんに改めて感謝です。

ABOUT ME
hirayuki
今年で社会人3年目になります。 日々体当たりで仕事を覚えています。 テーマはIT・教育です。 少しでも技術に親しんでもらえるよう、noteで4コマ漫画も書いています。 https://note.mu/hirayuki