{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Summary\n* LightGBM regressor.\n* A customized loss function (a weighted combo of fair+huber+RMSE+MAE).\n* Testing a bootstrapping approach by adding an estimate of `time_after_failure` to the features, to which the true `time_to_failure` is highly correlated to. \n* Conclusion: bootstrapping will lead to leakage and deteriorate our model. \n\nReference:\n* [My tries to find magic features](https://www.kaggle.com/scaomath/lanl-earthquakes-try-to-find-magic-features)\n* [BigIronSphere's data augmentation](https://www.kaggle.com/bigironsphere/basic-data-augmentation-feature-reduction)\n* [Andrew's \"Even more features\" notebook](https://www.kaggle.com/artgor/even-more-features)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\npaths = os.listdir(\"../input\")\nprint(paths)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set()\n\nfrom tqdm import tqdm_notebook\n\nfrom lightgbm import LGBMRegressor\nimport lightgbm as lgb\nfrom sklearn.model_selection import KFold\nfrom scipy import stats\n\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_features_ = pd.read_csv('../input/lanl-earthquake-nonmagic-features/train_X.csv')\ntrain_y_ = pd.read_csv('../input/lanl-earthquake-nonmagic-features/train_y.csv')\ntrain_y_extra_ = pd.read_csv('../input/lanl-earthquake-nonmagic-features/train_y_extra.csv')\ntest_features_ = pd.read_csv('../input/lanl-earthquake-nonmagic-features/test_X.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = train_features_.copy()\nX_test = test_features_.copy()\ny_train = train_y_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(X_train['mean'],label='Training samples')\nsns.distplot(X_test['mean'],label='Testing samples')\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(X_train['fft_skew'],label='Training samples')\nsns.distplot(X_test['fft_skew'],label='Testing samples')\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# time_after_failure"},{"metadata":{"trusted":true},"cell_type":"code","source":"rows = 150_000\ny_train_new = pd.DataFrame(columns=['taf'], dtype=np.float64, index = y_train.index)\ny_train_new['taf'] = train_y_extra_['time_after_failure']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# most correlated features with time_after_failure\ncols = list(np.abs(X_train.corrwith(y_train_new['taf'])).sort_values(ascending=False).head(50).index)\ncols[::5]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_plot = ['num_peaks_5',\n             'std_0_to_10',\n             'fft_100_roll_std_70',\n             'std_0_to_10',\n             'iqr',\n             'energy_spectra_9306hz']\n_, ax1 = plt.subplots(3, 2, figsize=(20, 12))\nax1 = ax1.reshape(-1)\n\nfor i, col in enumerate(cols_plot):\n    ax1[i].plot(X_train[col], color='blue')\n    ax1[i].set_title(col)\n    ax1[i].set_ylabel(col, color='b')\n\n    ax2 = ax1[i].twinx()\n    ax2.plot(y_train_new['taf'], color='g', linewidth=2)\n    ax2.set_ylabel('time_after_failure', color='g')\n    ax2.legend([col, 'time_after_failure'], loc= 'upper right')\n    ax2.grid(False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# LightGBM model"},{"metadata":{"trusted":true},"cell_type":"code","source":"coef = [0.35, 0.5, 0.05, 0.1]\n\ndef custom_objective(y_true, y_pred):\n    \n    # fair\n    c = 0.5\n    residual = y_pred - y_true\n    grad = c * residual /(np.abs(residual) + c)\n    hess = c ** 2 / (np.abs(residual) + c) ** 2\n    \n    # huber\n    h = 1.2  #h is delta in the formula\n    scale = 1 + (residual / h) ** 2\n    scale_sqrt = np.sqrt(scale)\n    grad_huber = residual / scale_sqrt\n    hess_huber = 1 / scale / scale_sqrt\n\n    # rmse grad and hess\n    grad_rmse = residual\n    hess_rmse = 1.0\n\n    # mae grad and hess\n    grad_mae = np.array(residual)\n    grad_mae[grad_mae > 0] = 1.\n    grad_mae[grad_mae <= 0] = -1.\n    hess_mae = 1.0\n\n    return coef[0] * grad + coef[1] * grad_huber + coef[2] * grad_rmse + coef[3] * grad_mae, \\\n           coef[0] * hess + coef[1] * hess_huber + coef[2] * hess_rmse + coef[3] * hess_mae\n\n\ndef logcosh_objective(y_true, y_pred):\n    d = y_pred - y_true \n    grad = np.tanh(d)/y_true\n    hess = (1.0 - grad*grad)/y_true\n    return grad, hess\n\n\ndef huber_objective(y_true, y_pred):\n    d = y_pred - y_true\n    h = 1.2  #h is delta in the formula\n    scale = 1 + (d / h) ** 2\n    scale_sqrt = np.sqrt(scale)\n    grad = d / scale_sqrt\n    hess = 1 / scale / scale_sqrt\n    return grad, hess","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"params = {'bagging_fraction': 0.71,\n         'boosting': 'gbdt',\n         'feature_fraction': 0.94,\n         'lambda_l1': 3.131216244016188,\n         'lambda_l2': 2.4124061313905836,\n         'learning_rate': 0.049886848207269734,\n         'max_bin': 193,\n         'max_depth': 15,\n         'metric': 'MAE',\n         'min_data_in_bin': 167,\n         'min_data_in_leaf': 62,\n         'min_gain_to_split': 2.07,\n         'num_leaves': 38,\n         'objective': custom_objective,\n         'subsample': 0.9133120405819966}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_prediction(X, y, X_test=None, \n                   n_fold=5, random_state=1127, eval_metric='mae',\n                   verbose=0, early_stopping_rounds=1000):\n    '''\n    X: dataframe\n    y: series or dataframe\n    params: global variable\n    '''\n    folds = KFold(n_splits=n_fold, shuffle=True, random_state=random_state)\n    n_fold_disp = n_fold // 2\n    \n    if X_test is None:\n        X_test = X\n    \n    pred_train = np.zeros(len(X))\n    pred_oof = np.zeros(len(X))\n    pred_test = np.zeros(len(X_test))\n\n    for fold_, (trn_idx, val_idx) in enumerate(folds.split(X,y)):\n        \n        if (fold_% n_fold_disp == 0 or fold_ == n_fold-1) and verbose > 0:\n            print(\"Fold {}\".format(fold_))\n\n        X_tr, X_val = X.iloc[trn_idx], X.iloc[val_idx]\n        y_tr, y_val = y.values[trn_idx].reshape(-1), y.values[val_idx].reshape(-1)\n\n        model = lgb.LGBMRegressor(**params, n_estimators = 10000, n_jobs = -1)\n        model.fit(X_tr, y_tr, \n                  eval_set=[(X_tr, y_tr), (X_val, y_val)], \n                  eval_metric=eval_metric,\n                  verbose=verbose, \n                  early_stopping_rounds=early_stopping_rounds)\n\n        #predictions\n        pred_train += model.predict(X, num_iteration=model.best_iteration_) / folds.n_splits\n        pred_test += model.predict(X_test, num_iteration=model.best_iteration_) / folds.n_splits\n        pred_oof[val_idx] = model.predict(X_val, num_iteration=model.best_iteration_)\n        \n    return pred_train, pred_test, pred_oof, model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Trial run on training data using selected features\n\n* Compute `time_to_failure`'s estimate based on `X_train`.\n* Estimate `time_after_failure` for training set.\n* Adding `time_after_failure`'s estimate into `X_train`.\n* Compute `time_to_failure` for training set again."},{"metadata":{"trusted":true},"cell_type":"code","source":"features_select = ['abs_min',\n             'abs_q01',\n             'abs_q05',\n             'abs_trend',\n             'autocorrelation_10',\n             'autocorrelation_1000',\n             'autocorrelation_500',\n             'av_change_abs_roll_mean_10',\n             'av_change_abs_roll_std_10',\n             'av_change_abs_roll_std_100',\n             'av_change_abs_roll_std_1000',\n             'avg_first_10000',\n             'avg_last_10000',\n             'avg_last_5000',\n             'c3_5',\n             'c3_500',\n             'classic_sta_lta1_mean',\n             'classic_sta_lta4_mean',\n             'classic_sta_lta5_mean',\n             'count_big_50000_less_threshold_5',\n             'energy_spectra_2640hz',\n             'energy_spectra_9306hz',\n             'energy_spectra_norm_109306hz_denoised',\n             'energy_spectra_norm_129306hz',\n             'energy_spectra_norm_135973hz_denoised',\n             'energy_spectra_norm_142640hz_denoised',\n             'energy_spectra_norm_149306hz_denoised',\n             'energy_spectra_norm_15973hz_denoised',\n             'energy_spectra_norm_169306hz_denoised',\n             'energy_spectra_norm_22640hz_denoised',\n             'energy_spectra_norm_29306hz_denoised',\n             'energy_spectra_norm_35973hz_denoised',\n             'energy_spectra_norm_42640hz_denoised',\n             'energy_spectra_norm_49306hz_denoised',\n             'energy_spectra_norm_55973hz_denoised',\n             'energy_spectra_norm_62640hz_denoised',\n             'energy_spectra_norm_89306hz_denoised',\n             'energy_spectra_norm_9306hz_denoised',\n             'energy_spectra_norm_95973hz_denoised',\n             'fft_1000_roll_std_20',\n             'fft_1000_roll_std_25',\n             'fft_1000_roll_std_70',\n             'fft_100_roll_std_1',\n             'fft_100_roll_std_20',\n             'fft_100_roll_std_70',\n             'fft_100_roll_std_75',\n             'fft_10_roll_std_75',\n             'fft_mean_change_rate',\n             'fft_min',\n             'fft_min_roll_mean_100',\n             'fft_min_roll_std_100',\n             'fft_skew',\n             'fft_skew_first_50000',\n             'fft_spkt_welch_density_100',\n             'fft_spkt_welch_density_5',\n             'fft_spkt_welch_density_50',\n             'fft_time_rev_asym_stat_10',\n             'fft_time_rev_asym_stat_100',\n             'iqr',\n             'kstat_3',\n             'mad',\n             'max_first_5000',\n             'max_last_10000',\n             'max_roll_mean_1000',\n             'max_to_min',\n             'mean_change_abs',\n             'med',\n             'min_last_10000',\n             'min_roll_std_100',\n             'num_crossings',\n             'num_peaks_10',\n             'num_peaks_5',\n             'q01_roll_mean_100',\n             'q01_roll_std_10',\n             'q01_roll_std_1000',\n             'q05_roll_mean_100',\n             'q05_roll_std_1000',\n             'q95',\n             'q95_roll_mean_10',\n             'q95_roll_mean_100',\n             'q95_roll_std_1000',\n             'q99_roll_mean_1000',\n             'skew',\n             'std_0_to_10',\n             'std_first_5000',\n             'std_neg_10_to_0',\n             'std_neg_2_to_2',\n             'std_neg_5_to_5',\n             'std_roll_mean_1000']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_select = X_train[features_select].copy()\nX_test_select = X_test[features_select].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ttf_pred_orig, ttf_test_orig, ttf_pred_oof, model1= get_prediction(X_train_select, y_train,\n                                        X_test = X_test_select)\n\nprint(\"Training MAE for time_to_failure is {:.7f}\"\\\n          .format(np.abs(y_train['time_to_failure'] - ttf_pred_orig).mean())) \nprint(\"OOF MAE for time_to_failure is {:.7f}\"\\\n          .format(np.abs(y_train['time_to_failure'] - ttf_pred_oof).mean())) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,4))\nsns.distplot(y_train['time_to_failure'] , color=\"red\", label=\"y_train\")\nsns.distplot(pd.DataFrame(ttf_test_orig)[0], color=\"skyblue\", label=\"lgb OOF pred for y_train\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"taf_pred, _, taf_pred_oof, model2 = get_prediction(X_train_select, y_train_new)\n    \nprint(\"Training MAE for time_after_failure is {:.7f}\"\\\n      .format((y_train_new['taf'] - taf_pred).abs().mean()))\nprint(\"OOF MAE for time_after_failure is {:.7f}\"\\\n          .format(np.abs(y_train_new['taf'] - taf_pred_oof).mean())) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Checking the correlation of the `time_after_failure` (including the estimate by LGB regressor) and the `time_to_failure`."},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train['time_to_failure'].corr(y_train_new['taf'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train['time_to_failure'].corr(pd.Series(taf_pred))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"_, ax = plt.subplots(2,1, figsize=(12,8))\n\n\nax[0].plot(taf_pred_oof, color='r', label='taf prediction')\nax[0].plot(y_train_new.values, color='b', label='time_after_failure', linewidth = 2)\nax[0].set_ylabel('taf', color='orange')\nax[0].autoscale(axis='x',tight=True)\nax[0].set_title(\"OOF LightGBM prediction vs TAF\")\nax[0].legend(loc='best');\n\nax[1].plot(ttf_pred_oof, color='orange', label='ttf prediction')\nax[1].plot(y_train['time_to_failure'], color='b', label='time_to_failure', linewidth = 2)\nax[1].set_ylabel('ttf', color='r')\nax[1].autoscale(axis='x',tight=True)\nax[1].set_title(\"OOF LightGBM prediction vs TTF\")\nax[1].legend(loc='best');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualization of important features"},{"metadata":{},"cell_type":"markdown","source":"Not surprising that the most important ones are all FFT-based features."},{"metadata":{"trusted":true},"cell_type":"code","source":"important_feature_index = np.argsort(model1.feature_importances_)[::-1]\ncols = X_train.columns[important_feature_index[:100]]\ncols[:20]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_top8 = cols[:8]\n_, ax1 = plt.subplots(4, 2, figsize=(20, 20))\nax1 = ax1.reshape(-1)\n\nfor i, col in enumerate(cols_top8):\n    ax1[i].plot(X_train[col], color='blue')\n    ax1[i].set_title(col)\n    ax1[i].set_ylabel(col, color='b')\n\n    ax2 = ax1[i].twinx()\n    ax2.plot(y_train['time_to_failure'], color='g', linewidth=2)\n    ax2.set_ylabel('time_to_failure', color='g')\n    ax2.legend([col, 'time_to_failure'], loc= 'upper right')\n    ax2.grid(False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Adding `time_after_failure` estimate to training"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train['taf_estimate'] = taf_pred_oof\nttf_pred,_,ttf_pred_oof, model3 = get_prediction(X_train, y_train)\n\nprint(\"\\nTraining MAE for time_to_failure is {:.7f}\"\\\n      .format(np.abs(y_train['time_to_failure'] - ttf_pred).mean()))\nprint(\"OOF MAE for time_to_failure is {:.7f}\"\\\n          .format(np.abs(y_train['time_to_failure'] - ttf_pred_oof).mean())) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"important_feature_index = np.argsort(model3.feature_importances_)[::-1]\nX_train.columns[important_feature_index[:10]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is not surprising after adding `time_after_failure` estimate to the features, it becomes the most important one..."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,4))\nsns.distplot(y_train['time_to_failure'] , color=\"red\", label=\"y_train\")\nsns.distplot(pd.DataFrame(ttf_pred_oof)[0], color=\"skyblue\", label=\"lgb pred for y_train\")\nplt.title(\"TTF vs LGB OOF prediction after 1 iteration of bootstrap\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Bootstrapping effectivity test\n\nReload `X_train` and `X_test`, and do a manual split `X_cv` to test whether the bootstrapping works or not."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = X_train_select.copy()\nX_test = X_test_select.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_bsp_train = X_train[1200:]\nX_bsp_cv = X_train[:1200]\ny_bsp_train = y_train[1200:]\ny_bsp_cv = y_train[:1200]\ny_bsp_train_new = y_train_new[1200:]\ny_bsp_cv_new = y_train_new[:1200]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 6\nrandom_state = 802","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Step 0: get reference time_to_failure testing error"},{"metadata":{"trusted":true},"cell_type":"code","source":"ttf_est_train, ttf_est_cv, ttf_est_oof, _ = get_prediction(X_bsp_train, \n                                            y_bsp_train, \n                                            X_test=X_bsp_cv, \n                                            n_fold=n_fold,\n                                            random_state=random_state)\n\nprint(\"TTF CV MAE: {:.7f}\"\\\n      .format(np.abs(y_bsp_cv['time_to_failure'] - ttf_est_cv).mean())) \nprint(\"TTF OOF MAE: {:.7f}\"\\\n      .format(np.abs(y_bsp_train['time_to_failure'] - ttf_est_oof).mean()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Step 1: estimate time_after_failure for CV and train"},{"metadata":{"trusted":true},"cell_type":"code","source":"taf_est_train, taf_est_cv, taf_est_oof, _ = get_prediction(X_bsp_train, \n                                             y_bsp_train_new, \n                                             X_test=X_bsp_cv, \n                                             n_fold=n_fold,\n                                             random_state=random_state)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Step 2: estimate time_to_failure for test and train"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_bsp_train = X_bsp_train.assign(taf_estimate = taf_est_oof);\nX_bsp_cv = X_bsp_cv.assign(taf_estimate = taf_est_cv);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ttf_est_train, ttf_est_cv, ttf_est_oof, model_trial  = get_prediction(X_bsp_train, y_bsp_train,\n                                                X_test=X_bsp_cv, \n                                                n_fold=n_fold,random_state=random_state)\n\nprint(\"TTF CV MAE after adding TAF to training features: {:.7f}\"\\\n      .format(np.abs(y_bsp_cv['time_to_failure'] - ttf_est_cv).mean())) \nprint(\"TTF OOF MAE after adding TAF to training features: {:.7f}\"\\\n      .format(np.abs(y_bsp_train['time_to_failure'] - ttf_est_oof).mean())) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature importance for CV set"},{"metadata":{"trusted":true},"cell_type":"code","source":"from eli5 import show_weights\nfrom eli5.sklearn import PermutationImportance\nfrom IPython.display import display","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"perm = PermutationImportance(model_trial, random_state=1).fit(X_bsp_cv,y_bsp_cv)\ndisplay(show_weights(perm, feature_names = X_bsp_cv.columns.tolist()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Step 3: iteratively updating ttf and taf like a leapfrog scheme\n(This is problematic) Now we can add the estimate of `time_to_failure` based on the `time_after_failure` estimate to the training features. When estimating `time_to_failure`, we drop `time_to_failure`'s estimate feature, and vice versa for `time_after_failure`. In this way, we always estimate one another in a bootstrapping fashion."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_bsp_train['ttf_estimate'] = ttf_est_oof\nX_bsp_cv['ttf_estimate'] = ttf_est_cv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_bootstrap = 2\ncoeff_bsp = [0.4, 0.6] # coeff for [old, new]\n\nfor b in tqdm_notebook(range(num_bootstrap)):\n    \n    print(\"\\nBootstrapping iteration {0}.\".format(b+1))\n    \n    # reset time_after_failure values\n    random_state = b\n    taf_est_train, taf_est_cv, taf_est_oof, model_taf = \\\n    get_prediction(X_bsp_train.drop(columns=['taf_estimate']), \n                   y_bsp_train_new, \n                   X_test=X_bsp_cv.drop(columns=['taf_estimate']), \n                   n_fold=n_fold,random_state=random_state)\n    \n    print(\"TAF training OOF MAE: {:.7f}\"\\\n      .format(np.abs(y_bsp_train_new['taf'] - taf_est_oof).mean())) \n    print(\"TAF CV MAE: {:.7f}\"\\\n      .format(np.abs(y_bsp_cv_new['taf'] - taf_est_cv).mean())) \n    print(\"Mean abs change in TAF feature for CV: {:.7f}\"\\\n      .format(np.abs(X_bsp_cv['taf_estimate'] - taf_est_cv).mean())) \n    \n    # update the time_after_failure feature\n    X_bsp_train['taf_estimate'] = coeff_bsp[0]*X_bsp_train['taf_estimate'] + coeff_bsp[1]*taf_est_oof\n    X_bsp_cv['taf_estimate'] = coeff_bsp[0]*X_bsp_cv['taf_estimate'] + coeff_bsp[1]*taf_est_cv\n    \n    # reset time_to_failure value\n    ttf_est_train, ttf_est_cv, ttf_est_oof, model = \\\n    get_prediction(X_bsp_train.drop(columns=['ttf_estimate']), \n                   y_bsp_train, \n                   X_test=X_bsp_cv.drop(columns=['ttf_estimate']), \n                   n_fold=n_fold,random_state=random_state)\n    \n    print(\"TTF training OOF MAE: {:.7f}\"\\\n      .format(np.abs(y_bsp_train['time_to_failure'] - ttf_est_oof).mean())) \n    print(\"TTF CV MAE: {:.7f}\"\\\n      .format(np.abs(y_bsp_cv['time_to_failure'] - ttf_est_cv).mean())) \n    print(\"Mean abs change in TTF feature for CV: {:.7f}\"\\\n      .format(np.abs(X_bsp_cv['ttf_estimate'] - ttf_est_cv).mean())) \n    \n    # updating the time_to_failure feature in order to bootstrap time_after_failure feature\n    X_bsp_train['ttf_estimate'] = coeff_bsp[0]*X_bsp_train['ttf_estimate'] + coeff_bsp[1]*ttf_est_oof\n    X_bsp_cv['ttf_estimate'] = coeff_bsp[0]*X_bsp_cv['ttf_estimate'] + coeff_bsp[1]*ttf_est_cv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"important_feature_index = np.argsort(model.feature_importances_)[::-1]\nX_bsp_train.columns[important_feature_index[:20]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Observation\n\nEven as the second target for bootstrapping becomes the most important features, the MAE is getting worse and worse."},{"metadata":{"trusted":true},"cell_type":"code","source":"_, ax = plt.subplots(2,1, figsize=(12,8))\nax[0].plot(taf_est_cv, color='r', label='taf prediction')\nax[0].plot(y_bsp_cv_new['taf'], color='b', label='time_after_failure', linewidth = 2)\nax[0].set_ylabel('taf', color='orange')\nax[0].set_title(\"LightGBM CV prediction vs TAF\")\nax[0].legend(loc='best');\n\nax[1].plot(ttf_est_cv, color='orange', label='ttf prediction')\nax[1].plot(y_bsp_cv['time_to_failure'], color='b', label='time_to_failure', linewidth = 2)\nax[1].set_ylabel('ttf', color='r')\nax[1].set_title(\"LightGBM CV prediction vs TTF\")\nax[1].legend(loc='best');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Bootstrapping for actual testing data"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = train_features_[features_select]\nX_test = test_features_[features_select]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Final bootstrapping is performed based on {} features.\".format(X_train.shape[-1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 6\nrandom_state = 1127\n\nparams = {'bagging_fraction': 0.51,\n         'boosting': 'gbdt',\n         'feature_fraction': 0.9,\n         'lambda_l1': 2,\n         'lambda_l2': 0.03,\n         'learning_rate': 0.1,\n         'max_bin': 48,\n         'max_depth': 8,\n         'metric': 'MAE',\n         'min_data_in_bin': 57,\n         'min_data_in_leaf': 11,\n         'min_gain_to_split': 0.53,\n         'num_leaves': 83,\n         'objective': custom_objective,\n         'subsample': 0.55}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Step 1: get time_after_failure estimate\ntaf_est_train, taf_est_test, taf_est_oof, _ = get_prediction(X_train, y_train_new, \n                                               X_test=X_test, \n                                               n_fold=n_fold,random_state=random_state)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Step 2: \nX_train['taf_estimate'] = taf_est_oof\nX_test['taf_estimate'] = taf_est_test\n\nttf_est_train, ttf_est_test, ttf_est_oof, _   = get_prediction(X_train, y_train, \n                                                X_test=X_test, \n                                                n_fold=n_fold,random_state=random_state)\nttf_est_test_bsp1 = ttf_est_test\nprint(\"TTF OOF MAE after adding TAF to training features: {:.7f}\"\\\n      .format((np.abs(y_train['time_to_failure'] - ttf_est_oof)).mean())) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Step 3\nX_train['ttf_estimate'] = ttf_est_oof\nX_test['ttf_estimate'] = ttf_est_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Step 3:\nnum_bootstrap = 20\ncoeff_bsp = [0.4, 0.6] # coeff for [old, new]\n\nfor b in tqdm_notebook(range(num_bootstrap)):\n    \n    print(\"\\nBootstrapping iteration {0}.\".format(b+1))\n    random_state = b\n    # reset time_after_failure values\n    taf_est_train, taf_est_test, taf_est_oof, _ = \\\n    get_prediction(X_train.drop(columns=['taf_estimate']), \n                   y_train_new, \n                   X_test=X_test.drop(columns=['taf_estimate']), \n                   n_fold=n_fold,random_state=random_state)\n    \n    print(\"TAF training OOF MAE: {:.7f}\"\\\n      .format(np.abs(y_train_new['taf'] - taf_est_oof).mean())) \n#     print(\"Mean abs change in TAF feature for training: {:.7f}\"\\\n#       .format(np.abs(X_train['taf_estimate'] - taf_est_train).mean())) \n#     print(\"Correlation between TAF estimate and time_to_failure: {:.5f}\"\\\n#           .format(y_train['time_to_failure'].corr(pd.Series(taf_est_train))))\n    \n    # update the time_after_failure feature\n    X_train['taf_estimate'] = coeff_bsp[0]*X_train['taf_estimate'] + coeff_bsp[1]*taf_est_oof\n    X_test['taf_estimate'] = coeff_bsp[0]*X_test['taf_estimate'] + coeff_bsp[1]*taf_est_test\n    \n    # reset time_to_failure value\n    ttf_est_train, ttf_est_test, ttf_est_oof, model = \\\n    get_prediction(X_train.drop(columns=['ttf_estimate']), \n                   y_train, \n                   X_test=X_test.drop(columns=['ttf_estimate']), \n                   n_fold=n_fold,random_state=random_state)\n    \n    print(\"TTF training OOF MAE: {:.7f}\"\\\n      .format((np.abs(y_train['time_to_failure'] - ttf_est_oof)).mean())) \n#     print(\"Mean abs change in TTF feature for training: {:.7f}\"\\\n#       .format((np.abs(X_train['ttf_estimate'] - ttf_est_train)).mean())) \n#     print(\"Correlation between TTF estimate and time_after_failure: {:.5f}\"\\\n#           .format(y_train_new['taf'].corr(pd.Series(ttf_est_train))))\n    \n    # updating the time_to_failure feature in order to bootstrap time_after_failure feature\n    X_train['ttf_estimate'] = coeff_bsp[0]*X_train['ttf_estimate'] + coeff_bsp[1]*ttf_est_oof\n    X_test['ttf_estimate'] = coeff_bsp[0]*X_test['ttf_estimate'] + coeff_bsp[1]*ttf_est_test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusion:\nAfter several iteration of bootstrapping, the leakage is severe."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,4))\nsns.distplot(y_train['time_to_failure'] , color=\"red\", label=\"time_to_failure\")\nsns.distplot(pd.DataFrame(ttf_est_oof)[0], color=\"skyblue\", label=\"LGB bootstrap OOF training pred\")\n# sns.distplot(pd.DataFrame(ttf_est_train)[0], label=\"LGB bootstrap training pred\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,4))\nsns.distplot(pd.DataFrame(ttf_est_test)[0] , color=\"skyblue\", \n             label=\"Bootstrapped prediction for X_test \")\nsns.distplot(pd.DataFrame(ttf_est_test_bsp1)[0] , color=\"orange\", \n             label=\"After adding taf estimate to feature once\")\nsns.distplot(pd.DataFrame(ttf_test_orig)[0] , color=\"green\", \n             label=\"Original prediction for X_test\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 6))\nplt.plot(y_train.values, color='b', label='time_to_failure', linewidth = 2)\nplt.plot(ttf_est_oof, color='r', label='LGB estimate')\nplt.legend(loc='best')\nplt.autoscale(axis='x',tight=True)\nplt.title('TTF vs LGB OOF prediction after {0} iteration of bootstrap'.format(num_bootstrap));","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', \n                         index_col='seg_id')\nsubmission_nonbsp = submission.copy()\nsubmission_bsp1 = submission.copy() # only bootstrapping once\nsubmission['time_to_failure'] = ttf_est_test\nsubmission_bsp1['time_to_failure'] = ttf_est_test_bsp1\nsubmission_nonbsp['time_to_failure'] = ttf_test_orig\nsubmission.to_csv('submission.csv')\nsubmission_nonbsp.to_csv('submission_nonbsp.csv')\nsubmission_bsp1.to_csv('submission_bsp1.csv')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}