kaggleの雑記

kaggleのカーネルからアンサンブル結果の評価・可視化方法を学ぶ

コンペについて

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

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

今回読み込むカーネル

marcelo morenosさんのFE and ensemble MLP and LGBMを読み込もうと思います。

今回は特にアンサンブル学習を評価するコードが気になりました。ここの読み込みに力を入れて、その他は割愛気味とします。

データの読み込み

# read data
data = pd.read_csv('../input/data-without-drift/train_clean.csv')
data.head()

ここでdriftを取り除いたデータを参照しているようです。

どこから取得したデータか書いてありませんでしたが、ion switching コンペでは

https://www.kaggle.com/friedchips/clean-removal-of-data-drift/output

がトレンド情報を取り除いたカーネルとして有名なので、これだということにしておきます。

特徴量生成前準備

勾配の取得

def calc_gradients(s, n_grads=4):
    '''
    Calculate gradients for a pandas series. Returns the same number of samples
    '''
    grads = pd.DataFrame()
    
    g = s.values
    for i in range(n_grads):
        g = np.gradient(g)
        grads['grad_' + str(i+1)] = g
        
    return grads

lowバターワースフィルタを追加

def calc_low_pass(s, n_filts=10):
    '''
    Applies low pass filters to the signal. Left delayed and no delayed
    '''
    wns = np.logspace(-2, -0.3, n_filts)
    
    low_pass = pd.DataFrame()
    x = s.values
    for wn in wns:
        b, a = signal.butter(1, Wn=wn, btype='low')
        zi = signal.lfilter_zi(b, a)
        low_pass['lowpass_lf_' + str('%.4f' %wn)] = signal.lfilter(b, a, x, zi=zi*x[0])[0]
        low_pass['lowpass_ff_' + str('%.4f' %wn)] = signal.filtfilt(b, a, x)
        
    return low_pass

バターワースフィルタには詳しくはないが、周波数の整形処理の一種類の様子

Wikipediaを参照

Highバターワースフィルタを追加

def calc_high_pass(s, n_filts=10):
    '''
    Applies high pass filters to the signal. Left delayed and no delayed
    '''
    wns = np.logspace(-2, -0.1, n_filts)
    
    high_pass = pd.DataFrame()
    x = s.values
    for wn in wns:
        b, a = signal.butter(1, Wn=wn, btype='high')
        zi = signal.lfilter_zi(b, a)
        high_pass['highpass_lf_' + str('%.4f' %wn)] = signal.lfilter(b, a, x, zi=zi*x[0])[0]
        high_pass['highpass_ff_' + str('%.4f' %wn)] = signal.filtfilt(b, a, x)
        
    return high_pass

移動窓ごとに各数値処理

def calc_roll_stats(s, windows=[10, 50, 100, 500, 1000]):
    '''
    Calculates rolling stats like mean, std, min, max...
    '''
    roll_stats = pd.DataFrame()
    for w in windows:
        roll_stats['roll_mean_' + str(w)] = s.rolling(window=w, min_periods=1).mean()
        roll_stats['roll_std_' + str(w)] = s.rolling(window=w, min_periods=1).std()
        roll_stats['roll_min_' + str(w)] = s.rolling(window=w, min_periods=1).min()
        roll_stats['roll_max_' + str(w)] = s.rolling(window=w, min_periods=1).max()
        roll_stats['roll_range_' + str(w)] = roll_stats['roll_max_' + str(w)] - roll_stats['roll_min_' + str(w)]
        roll_stats['roll_q10_' + str(w)] = s.rolling(window=w, min_periods=1).quantile(0.10)
        roll_stats['roll_q25_' + str(w)] = s.rolling(window=w, min_periods=1).quantile(0.25)
        roll_stats['roll_q50_' + str(w)] = s.rolling(window=w, min_periods=1).quantile(0.50)
        roll_stats['roll_q75_' + str(w)] = s.rolling(window=w, min_periods=1).quantile(0.75)
        roll_stats['roll_q90_' + str(w)] = s.rolling(window=w, min_periods=1).quantile(0.90)
    
    # add zeros when na values (std)
    roll_stats = roll_stats.fillna(value=0)
             
    return roll_stats

今回のコンペで頻出しているよくあるラグ特徴量

kaggleのお手本のような時系列特徴量生成【Ion Switching】コンペについて https://www.kaggle.com/c/liverpool-ion-switching リバプールの大...

指数平滑移動平均

def calc_ewm(s, windows=[10, 50, 100, 500, 1000]):
    '''
    Calculates exponential weighted functions
    '''
    ewm = pd.DataFrame()
    for w in windows:
        ewm['ewm_mean_' + str(w)] = s.ewm(span=w, min_periods=1).mean()
        ewm['ewm_std_' + str(w)] = s.ewm(span=w, min_periods=1).std()
        
    # add zeros when na values (std)
    ewm = ewm.fillna(value=0)
        
    return ewm

指数平滑移動平均は、過去の価格よりも直近の価格になるほど比重を置いて計算された平均値です。

(https://www.smbcnikko.co.jp/terms/japan/si/J0628.html)

特徴量生成関数の実行(および分割)

def add_features(s):
    '''
    All calculations together
    '''
    
    gradients = calc_gradients(s)
    low_pass = calc_low_pass(s)
    high_pass = calc_high_pass(s)
    roll_stats = calc_roll_stats(s)
    ewm = calc_ewm(s)
    
    return pd.concat([s, gradients, low_pass, high_pass, roll_stats, ewm], axis=1)


def divide_and_add_features(s, signal_size=500000):
    '''
    Divide the signal in bags of "signal_size".
    Normalize the data dividing it by 15.0
    '''
    # normalize
    s = s/15.0
    
    ls = []
    for i in tqdm(range(int(s.shape[0]/signal_size))):
        sig = s[i*signal_size:(i+1)*signal_size].copy().reset_index(drop=True)
        sig_featured = add_features(sig)
        ls.append(sig_featured)
    
    return pd.concat(ls, axis=0)

add_featuresは今までの特徴量生成関数を実行してきただけのようです。

ただしdivide_and_add_featuresはsignal_sizeでtrainデータを幾つの単位で分割するかを決めることができようです。

訓練データと検証データに分割

# Get train and test data
x_train, x_test, y_train, y_test = train_test_split(df.values, data['open_channels'].values, test_size=0.2)

del data, df
print('x_train.shape=', x_train.shape)
print('x_test.shape=', x_test.shape)
print('y_train.shape=', y_train.shape)
print('y_test.shape=', y_test.shape)

testという名称が付けられていますが、kaggleで与えられた後のアンサンブル評価で利用する様子なので、検証と言えば検証データだと思っております。。ただ、こういうケースではテストデータと呼んだ方がいいのかもしれません。。

open_channelsの頻度をプロット

def get_class_weight(classes, exp=1):
    '''
    Weight of the class is inversely proportional to the population of the class.
    There is an exponent for adding more weight.
    '''
    hist, _ = np.histogram(classes, bins=np.arange(12)-0.5)
    class_weight = hist.sum()/np.power(hist, exp)
    
    return class_weight

class_weight = get_class_weight(y_train)
print('class_weight=', class_weight)
plt.figure()
plt.title('classes')
plt.hist(y_train, bins=np.arange(12)-0.5)
plt.figure()
plt.title('class_weight')
plt.bar(np.arange(11), class_weight)
plt.title('class_weight')

y_trainはopen_channelsに相当し、チャンネル数がいくつ空いているか、0~10まであるその頻度をプロットしています。

後でLightGBM.Datasetのweightの引数に利用しています。

https://lightgbm.readthedocs.io/en/latest/pythonapi/lightgbm.Dataset.html

https://lightgbm.readthedocs.io/en/latest/pythonapi/lightgbm.Dataset.html

https://lightgbm.readthedocs.io/en/latest/pythonapi/lightgbm.Dataset.html

KerasによるMLPの実装

def create_mpl(shape):
    '''
    Returns a keras model
    '''
    
    X_input = L.Input(shape)
    
    X = L.Dense(150, activation='relu')(X_input)
    X = L.Dense(150, activation='relu')(X)
    X = L.Dense(125, activation='relu')(X)
    X = L.Dense(75, activation='relu')(X)
    X = L.Dense(50, activation='relu')(X)
    X = L.Dense(25, activation='relu')(X)
    X = L.Dense(11, activation='softmax')(X)
    
    model = Model(inputs=X_input, outputs=X)
    
    return model


mlp = create_mpl(x_train[0].shape)
mlp.compile(optimizer='adam', loss='sparse_categorical_crossentropy', metrics=['sparse_categorical_accuracy'])
print(mlp.summary())

Deep Learningによる実装ですRelu関数とSoftmax関数、オプティマイザはadamです。

sparse_categorical_crossentropyは後で調べます。

MLPの学習

# fit the model
mlp.fit(x=x_train, y=y_train, epochs=30, batch_size=1024, class_weight=class_weight)

学習の進行結果をplotする

# plot history
plt.figure(1)
plt.plot(mlp.history.history['loss'], 'b', label='loss')
plt.xlabel('epochs')
plt.legend()
plt.figure(2)
plt.plot(mlp.history.history['sparse_categorical_accuracy'], 'g', label='sparse_categorical_accuracy')
plt.xlabel('epochs')
plt.legend()

学習が正しく進んでいることがわかります。

検証データのMLP予測結果を変数に格納

# predict on test
mlp_pred = mlp.predict(x_test)
print('mlp_pred.shape=', mlp_pred.shape)

LightGBMの学習用のデータとパラメータを確定

# build model
# lgb.Datasetを利用してLightGBM用の学習データオブジェクトに変更
dataset = lgb.Dataset(x_train, label=y_train, weight=class_weight[y_train])
params = {'learning_rate': 0.1,
          'objective': 'multiclass',
          'num_class': 11,
          'metric': 'multi_logloss',
          'verbose': 1}

LightGBMの学習

# fit the model
print('Training LGBM...')
gbc = lgb.train(params, dataset, 200, verbose_eval=10)
print('LGBM trained!')

検証データのLightGBM予測を変数に格納

# predict on test
gbc_pred = gbc.predict(x_test)
print('gbc_pred.shape=', gbc_pred.shape)

今回の本題:アンサンブルデータを評価する

アンサンブルをMLPとLightGBMで1%毎変えていき、最適な割合を探索する

# lists for keep results
f1s = []
alphas = []

# loop for every alpha
# tdqm関数をfor分のイテレータを包むと、プログレスバーを表示する
for alpha in tqdm(np.linspace(0,1,101)): # linspace: 0~1の間を101等分する
    # mlpとgbc(lgbm)の配分を変えている
    y_pred = alpha*mlp_pred + (1 - alpha)*gbc_pred
    # それぞれの配分毎の結果をy_testでスコアリングして、f1として残していく
    f1 = f1_score(y_test, np.argmax(y_pred, axis=-1), average='macro')
    f1s.append(f1)
    alphas.append(alpha)

# convert to numpy arrays
f1s = np.array(f1s)
alphas = np.array(alphas)

# get best_alpha f1sが最大だった時のindexをalphasの中から取得している
best_alpha = alphas[np.argmax(f1s)]

print('best_f1=', f1s.max())
print('best_alpha=', best_alpha)

tdqmは進捗状況を確認できるプログレスバーを表示させる

https://qiita.com/pontyo4/items/76145cb10e030ad8186a

割合ごとにどれくらいスコアに影響しているかplotする

plt.plot(alphas, f1s)
plt.title('f1_score for ensemble')
plt.xlabel('alpha')
plt.ylabel('f1_score')

alphaが0の時にLightGBMが10割になっていることを示している。

0.5で突発的に変わっているけどこんなもの??

正答率をプロットする

# Thanks to https://www.kaggle.com/marcovasquez/basic-nlp-with-tensorflow-and-wordcloud
def plot_cm(y_true, y_pred, title):
    figsize=(16,16)
    y_pred = y_pred.astype(int)
    # confusion_matrixを利用する、ラベルはy_trueから取得
    cm = confusion_matrix(y_true, y_pred, labels=np.unique(y_true))
    # confusion_matrixの合計値を取得し、全体を割り100を掛けている。つまり%表記にしている
    cm_sum = np.sum(cm, axis=1, keepdims=True)
    cm_perc = cm / cm_sum.astype(float) * 100
    annot = np.empty_like(cm).astype(str)
    nrows, ncols = cm.shape
    # それぞれにスコアリング結果を入れている
    for i in range(nrows):
        for j in range(ncols):
            c = cm[i, j]
            p = cm_perc[i, j]
            if i == j:
                s = cm_sum[i]
                annot[i, j] = '%.1f%%\n%d/%d' % (p, c, s)
            elif c == 0:
                annot[i, j] = ''
            else:
                annot[i, j] = '%.1f%%\n%d' % (p, c)
    # plotのためにDataFrameに格納している、indexとカラム名を決められる
    cm = pd.DataFrame(cm, index=np.unique(y_true), columns=np.unique(y_true))
    cm.index.name = 'Actual'
    cm.columns.name = 'Predicted'
    fig, ax = plt.subplots(figsize=figsize)
    plt.title(title)
    sns.heatmap(cm, cmap='viridis', annot=annot, fmt='', ax=ax)
f1_mlp = f1_score(y_test, np.argmax(mlp_pred, axis=-1), average='macro')
plot_cm(y_test, np.argmax(mlp_pred, axis=-1), 'Only MLP \n f1=' + str('%.4f' %f1_mlp))

これはMLPだけを利用した時の画像である。

横軸が予測、縦軸が実際の値でありその正答率が%表記で記されている。

open_channelsが0の時の正答率、1の時の正答率が最も高く、数が大きくなっていくほど正答率は低くなっている。

LightGBMだけの時のマトリクスをプロット

f1_gbc = f1_score(y_test, np.argmax(gbc_pred, axis=-1), average='macro')
plot_cm(y_test, np.argmax(gbc_pred, axis=-1), 'Only GBC \n f1=' + str('%.4f' %f1_gbc))

MLPとLightGBMを50%でアンサンブルさせた時のプロット

y_pred = best_alpha*mlp_pred + (1 - best_alpha)*gbc_pred
f1_ens = f1_score(y_test, np.argmax(y_pred, axis=-1), average='macro')
plot_cm(y_test, np.argmax(y_pred, axis=-1), 'Ensemble \n f1=' + str('%.4f' %f1_ens))

アンサンブルをさせるとopne_channelsが10の時などに大きな改善が見られます。

open_channels=10はweightの計算の時からもわかるように、そもそも評価する母数が他のopen_channels < 10に対して少ないです。

そのため一つ当たりの誤答が大きく影響してしまうのでしょう。アンサンブルが50%の時に跳ね上がったことも、それが原因かと思います。

提出

def submit_result(mlp, gbc, alpha):
    
    print('Reading data...')
    data = pd.read_csv('../input/data-without-drift/test_clean.csv')
    
    print('Feature engineering...')
    df = divide_and_add_features(data['signal'])
    
    print('Predicting MLP...')
    mlp_pred = mlp.predict(df.values)
    
    print('Predicting GBC...')
    gbc_pred = gbc.predict(df.values)
    
    print('Ensembling...')
    y_pred = alpha*mlp_pred + (1 - alpha)*gbc_pred
    y_pred = np.argmax(y_pred, axis=-1)
    
    print('Writing submission...')
    submission = pd.DataFrame()
    submission['time'] = data['time']
    submission['open_channels'] = y_pred
    submission.to_csv('submission.csv', index=False, float_format='%.4f')
    
    print('Submission finished!')
submit_result(mlp, gbc, best_alpha)
ABOUT ME
hirayuki
今年で社会人3年目になります。 日々体当たりで仕事を覚えています。 テーマはIT・教育です。 少しでも技術に親しんでもらえるよう、noteで4コマ漫画も書いています。 https://note.mu/hirayuki