{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Introduction\n\nIn this notebook I will share main ideas that helped me to get 130th place in LANL earthquake competition. Further I will:\n\n1. Give a brief note on how I calculated features;\n2. Share my CV-strategy;\n3. Take a look at very simple models which led me to the main insight;\n4. Outline my model;\n5. Explain how I selected features for my models.\n\nTo have an idea about the competition itself, please, go and look at this [beautiful kernel](https://www.kaggle.com/allunia/shaking-earth)."},{"metadata":{"trusted":true},"cell_type":"code","source":"import warnings\nwarnings.filterwarnings(\"ignore\")\n\nimport pandas as pd\nimport numpy as np\n\nfrom os.path import join\nfrom os import listdir\nimport os\n\nimport matplotlib\nimport matplotlib.pyplot as plt\n\nfrom tqdm import tqdm_notebook as tqdm\nfrom sklearn.metrics import roc_curve, auc\n\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.svm import NuSVR\nimport lightgbm as lgb\n\nfrom sklearn.metrics import mean_absolute_error as mae\nfrom itertools import combinations\n\nBLUE = '#3CBEC5'\nRED = '#EF758A'\nGREEN = '#19B278'\n\nfrom IPython.display import Image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 1. Feature generation\n\nAll of the features I used were calculated for 150K points intervals, however, I also used overlapping of 20% (which is 30K points):"},{"metadata":{"trusted":true},"cell_type":"code","source":"Image(\"../input/feature-generation-frame/feature_generation_framework.png\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"You can see the groups of intervals, I will explain how to deal with them when building models later on in CV-strategy section. Ideas for the features to generate were taken from the following kernels:\n\n1. [Rolling statistics](https://www.kaggle.com/artgor/even-more-features)\n2. [Magnitude statistics from FFT](https://www.kaggle.com/vettejeep/masters-final-project-model-lb-1-392)\n3. [Bandpassed feautres](https://www.kaggle.com/vettejeep/masters-final-project-model-lb-1-392)\n\nTo summarise:\n\n* There were around **700 features** generated.\n* Approximately **20K train records** were obtained from overlapping approach."},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/processed/train_data.csv\")\ntest = pd.read_csv(\"../input/processed/test_data.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below you can find first rows of the train and test data sets."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Train data set:\", train.shape)\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Test data set:\")\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. CV-strategy\n\nAs you may remember, I calculated the features for intervals that were overlapping. If one of the overlapping intervals will be used in the model training and the second one in the model validation it may cause overfitting. To remove this effect we will split our intervals in groups of 5 consecutive intervals (<font color=\"blue\">\"blue\"</font> and <font color=\"red\">\"red\"</font> colors on the picture from the feature generation section is exactly about this) and take the whole group either for training or for validation. However, we will still have overlapping groups, and for that we will remove all <font color=\"red\">\"red\"</font> groups which gives us only non-overlapping <font color=\"blue\">\"blue\"</font> groups.\n\nThis strategy doubles the training data, nevertheless you can take different sizes of <font color=\"red\">\"red\"</font> and <font color=\"blue\">\"blue\"</font> groups to increase the size of the training data. Taking 20 consecutive intervals for <font color=\"blue\">\"blue\"</font> groups and 5 for <font color=\"red\">\"red\"</font> will give you three times more data than you might have without overlapping technique.\n\nNext step is to split our groups into folds. I did it randomly to obtain 10 folds, where each fold has some data from every earthquake.\n\n***Remark: taking larger blue groups will give us more data, but since we take intervals from one group together, we will have larger time intervals from the same quake. If you take extremely big blue groups you might end up with somewhat close to quake-CV, where each fold is one quake.***\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# split groups into blue and red\ntrain['fold'] = np.arange(len(train)) // 10\ntrain.loc[(np.arange(len(train)) % 10 > 4), 'fold'] = -1\n# remove red groups\ntrain = train[train.fold != -1].reset_index().drop('index', axis=1)\n\n# get 10 folds\nN = train.fold.max()\nK = 10\nk = N // K\nnp.random.seed(0)\npermutation = np.random.choice(N, N, replace=False)\nfolds = []\nfor i in range(K):\n    fold = permutation[i*k:(i+1)*k]\n    folds.append(fold)\n    tmp = train[~np.isin(train.fold, fold)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Simple models + Crucial insight\n\n1. Let's start with very simple linear models using 10-fold-CV. I will show only some of them:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_linear_model(feature):\n    train['prediction'] = 0\n    # building out-of-fold predictions\n    for i in range(K):\n        val_fold = folds[i]\n        cv_train = train[(~np.isin(train.fold, val_fold))]\n        cv_val = train[np.isin(train.fold, val_fold)]\n        feature_train = cv_train[feature].values\n        target_train = cv_train[TARGET].values\n        feature_val = cv_val[feature].values\n        target_val = cv_val[TARGET].values\n        \n        model = LinearRegression()\n        model.fit(feature_train.reshape(-1,1), target_train)\n\n        val_predictions = model.predict(feature_val.reshape(-1,1))        \n        train.loc[np.isin(train.fold, val_fold), 'prediction'] = val_predictions\n\n    # building model in full training data\n    cv_train = train.copy()\n    feature_train = cv_train[feature].values\n    target_train = cv_train[TARGET].values\n\n    model = LinearRegression()\n    model.fit(feature_train.reshape(-1,1), target_train)\n\n    # evaluating quality\n    score = mae(train[TARGET].values, train.prediction.values)\n    \n    # plotting\n    fig, ax = plt.subplots(figsize=(18,6))\n    \n    plt.subplot(121)\n    plt.plot(train[TARGET].values, c=RED, label='true')\n    plt.plot(train.prediction.values, c=BLUE, alpha=0.5, label='prediction')\n    plt.xlabel('time')\n    plt.ylabel('time_to_failure')\n    plt.legend()\n    plt.grid()\n    \n    plt.subplot(122)\n    plt.scatter(y=train[TARGET].values, x=feature_train, c=BLUE, alpha=0.1, label='true')\n    a = np.quantile(feature_train, 0.05)\n    b = np.quantile(feature_train, 0.95)\n    x = np.linspace(a, b, 10)\n    features = np.vstack([x]).T\n    preds = model.predict(features.reshape(-1,1))\n    plt.plot(features, preds, c=RED, alpha=1., label='linear model')\n    plt.xlabel(feature)\n    plt.ylabel('time_to_failure')\n    plt.legend()\n    plt.grid()\n    \n    print(\"MAE =\", score)\n    plt.show()\n    return score\n\nTARGET = 'time_to_failure'\nfor feature in train.columns[:-3]:\n    build_linear_model(feature)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some of the models are very good, for example, with linear model based on feature **roll_100_std_percentile_5** (5% quantile of rolling standard deviation with a window of size 100) we can achieve MAE=2.19. Moreover, if we take a closer look on the right plots, we can see that there are two different patterns: one for **time_to_failure < 0.32** and one for **time_to_failure > 0.32**.\n\nThis separation comes from the fact that nearly at the same time (0.29-0.32 seconds to failure) signal shows extremely high values, and probably, there is a different nature of the experiment before and after this point."},{"metadata":{"trusted":true},"cell_type":"code","source":"# 0.32 seconds separation\nfig, ax = plt.subplots(figsize=(20,10))\n\ntarget = train[TARGET].values[700:900]\nfeature = train['roll_100_std_percentile_5'].values[700:900]\nfeature = (feature - min(feature)) / (max(feature) - min(feature)) * 15\n\nplt.subplot(211)\nplt.plot(target, color=RED, label='time_to_failure')\nplt.plot(feature, color=BLUE, label='roll_100_std_percentile_5 (normalized)')\nfor split in [115]:\n    plt.axvline(x=split, color=GREEN, linestyle='--', label='0.32 seconds to failure')\n\nplt.title('0.32 seconds separation')\nplt.xlabel('record')\nplt.ylabel('value')\nplt.grid()\nplt.legend()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try to build separate models for records with **time_to_failure < 0.32** and for **time_to_failure > 0.32**."},{"metadata":{"trusted":true},"cell_type":"code","source":"feature = 'roll_100_std_percentile_5'\n\nfeature_values = train[feature].values\nthreshold = .32\ntarget = train[TARGET]\nindex1 = (target > threshold)\nindex2 = (target <= threshold)\n        \n# build simple model for points before threshold \nmodel1 = LinearRegression()\nmodel1.fit(feature_values[index1].reshape(-1,1), train[TARGET].values[index1])\n# build simple model for points after threshold \nmodel2 = LinearRegression()\nmodel2.fit(feature_values[index2].reshape(-1,1), train[TARGET].values[index2])\n\n# evaluate quality\npreds1 = model1.predict(feature_values.reshape(-1,1))\npreds1[preds1 < 0] = 0\npreds2 = model2.predict(feature_values.reshape(-1,1))\npreds2[preds2 < 0] = 0\npreds = model1.predict(feature_values.reshape(-1,1))\npreds[index2] = preds2[index2]\npreds[preds < 0] = 0\n\nscore = mae(train[TARGET].values, preds)\nprint(\"MAE =\", score)\n\n# plot\nfig, ax = plt.subplots(figsize=(18,6))\n\nplt.subplot(131)\nplt.plot(target, c=RED, label='true')\nplt.plot(preds, c=BLUE, alpha=0.5, label='prediction')\nplt.title('Two linear models')\nplt.xlabel('time')\nplt.ylabel('time_to_failure')\nplt.legend()\nplt.grid()\n\nplt.subplot(132)\nplt.scatter(y=target[index1], x=feature_values[index1], c=BLUE, alpha=0.1, label='true')\na = np.quantile(feature_values, 0.05)\nb = np.quantile(feature_values, 0.95)\nx = np.linspace(a, b, 10)\nfeatures = np.vstack([x]).T\npreds1 = model1.predict(features.reshape(-1,1))\nplt.plot(features, preds1, c=RED, alpha=1., label='linear model')\nplt.title('More than {} seconds to quake'.format(threshold))\nplt.xlabel(feature)\nplt.ylabel('time_to_failure')\nplt.legend()\nplt.grid()\n\nplt.subplot(133)\nplt.scatter(y=target[index2], x=feature_values[index2], c=BLUE, alpha=0.1, label='true')\na = np.quantile(feature_values, 0.05)\nb = np.quantile(feature_values, 0.95)\nx = np.linspace(a, b, 10)\nfeatures = np.vstack([x]).T\npreds2 = model2.predict(features.reshape(-1,1))\nplt.plot(features, preds2, c=RED, alpha=1., label='linear model')\nplt.title('Less than {} seconds to quake'.format(threshold))\nplt.xlabel(feature)\nplt.ylabel('time_to_failure')\nplt.legend()\nplt.grid()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we can also see that data after 0.32 seconds threshold behaves in non-linear way. Let's try to build a model based on 1/x transformation (in fact, transformation is 0.33-1/x) of the target."},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_complex_linear_model(feature):\n    feature_values = train[feature].values\n    threshold = 0.32\n    target = train[TARGET]\n    index1 = (target > threshold)\n    index2 = (target <= threshold)\n\n    # build simple model for points before threshold \n    model1 = LinearRegression()\n    model1.fit(feature_values[index1].reshape(-1,1), train[TARGET].values[index1])\n    # build simple model for points after threshold \n    model2 = LinearRegression()\n    model2.fit(feature_values[index2].reshape(-1,1), 1/(threshold + 0.01 - train[TARGET].values[index2]))\n    \n    # evaluate quality\n    preds1 = model1.predict(feature_values.reshape(-1,1))\n    preds1[preds1 < 0] = 0\n    preds2 = threshold + 0.01 - 1/model2.predict(feature_values.reshape(-1,1))\n    preds2[preds2 < 0] = 0\n\n    preds = model1.predict(feature_values.reshape(-1,1))\n    preds[index2] = preds2[index2]\n    preds[preds < 0] = 0\n\n    score = mae(target, preds)\n    print(\"MAE =\", score)\n\n    # plot\n    fig, ax = plt.subplots(figsize=(18,6))\n\n    plt.subplot(131)\n    plt.plot(target, c=RED, label='true')\n    plt.plot(preds, c=BLUE, alpha=0.5, label='prediction')\n    plt.title('Linear and non-linear model')\n    plt.xlabel('time')\n    plt.ylabel('time_to_failure')\n    plt.legend()\n    plt.grid()\n\n    plt.subplot(132)\n    plt.scatter(y=target[index1], x=feature_values[index1], c=BLUE, alpha=0.1, label='true')\n    a = np.quantile(feature_values, 0.01)\n    b = np.quantile(feature_values, 0.99)\n    x = np.linspace(a, b, 1000)\n    features = np.vstack([x]).T\n    preds1 = model1.predict(features.reshape(-1,1))\n    plt.plot(features, preds1, c=RED, alpha=1., label='linear model')\n    plt.title('More than {} seconds to quake'.format(threshold))\n    plt.xlabel(feature)\n    plt.ylabel('time_to_failure')\n    plt.legend()\n    plt.grid()\n\n    plt.subplot(133)\n    plt.scatter(y=target[index2], x=feature_values[index2], c=BLUE, alpha=0.1, label='true')\n    a = np.quantile(feature_values, 0.1)\n    b = np.quantile(feature_values, 0.999)\n    x = np.linspace(a, b, 1000)\n    features = np.vstack([x]).T\n    preds2 = threshold + 0.01 - 1/model2.predict(features.reshape(-1,1))\n    plt.plot(features, preds2, c=RED, alpha=1., label='linear model')\n    plt.title('Less than {} seconds to quake'.format(threshold))\n    plt.xlabel(feature)\n    plt.ylabel('time_to_failure')\n    plt.legend()\n    plt.grid()\n\n    plt.show()\n    \nplot_complex_linear_model(feature)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Quality decreased a little, but the plot for the data after 0.32 seconds threshold is much nicer.\n\nAll in all, we can see a huge improvement in the model from **MAE=2.19** to **MAE=1.96**. \n\nA careful reader will say that we don't know in advance when the signal is before or after this 0.32 seconds threshold. Correct! So we need to build an extra classification model."},{"metadata":{},"cell_type":"markdown","source":"## 4. Final model outline\n\nThe final model consists of:\n\n1. Classifier built by LightGBM to define if the signal is before or after 0.32 seconds threshold.\n2. LightGBM model built only on data with time_to_failure > 0.32\n3. NuSVR model based on roll_100_std_percentile_5 feature built only on data with time_to_failure < 0.32\n\nIn the end we select an optimal threshold for classifier to choose wich model to choose LightGBM or NuSVR."},{"metadata":{},"cell_type":"markdown","source":"#### Classifier"},{"metadata":{"trusted":true},"cell_type":"code","source":"train['classifier_32'] = 0\ntrain['target'] = train.time_to_failure < 0.32\nTARGET = 'target'\nfeatures = [\n 'mag_freq0_percentile_50',\n 'bp5_mad',\n 'bp1_interquantile_range_5_95',\n 'evol_bp2_interquantile_range_10_90',\n 'roll_4096_mean_std',\n 'roll_1000_std_percentile_5',\n 'bp0_kurtosis',\n 'bp1_percentile_10',\n 'evol_bp3_percentile_75'\n]\nparams = {\n    'num_leaves': 32,\n    'max_bin': 63,\n    'min_data_in_leaf': 50,\n    'learning_rate': 0.05,\n    'min_sum_hessian_in_leaf': 0.01,\n    'bagging_fraction': 0.8,\n    'bagging_freq': 1,\n    'feature_fraction': 1,\n    'min_gain_to_split': 0.02,\n    'max_depth': 6,\n    'save_binary': True, \n    'seed': 0,\n    'feature_fraction_seed': 0,\n    'bagging_seed': 0,\n    'drop_seed': 0,\n    'data_random_seed': 0,\n    'objective': 'binary',\n    'boosting_type': 'gbdt',\n    'verbose': 1,\n    'boost_from_average': True,\n    'metric': 'auc',\n    'is_unbalance': True\n}\n\nnum_boost_round=100\n\nfor i in range(K):\n    val_fold = folds[i]\n    cv_train = train[(~np.isin(train.fold, val_fold))]\n    cv_val = train[np.isin(train.fold, val_fold)]\n\n    xg_train = lgb.Dataset(cv_train[features].values, label=cv_train[TARGET].values)\n    xg_val = lgb.Dataset(cv_val[features].values, label=cv_val[TARGET].values)  \n\n    model = lgb.train(params, xg_train, num_boost_round=num_boost_round, verbose_eval=0)\n    val_predictions = model.predict(cv_val[features].values, num_iteration=num_boost_round)\n    train.loc[np.isin(train.fold, val_fold), 'classifier_32'] = val_predictions\n\nxg_train = lgb.Dataset(train[features].values, label=train[TARGET].values)\nmodel = lgb.train(params, xg_train, num_boost_round=num_boost_round)\ntest['classifier_32'] = model.predict(test[features].values, num_iteration=num_boost_round)\n\nfpr, tpr, _ = roc_curve(train['target'], train['classifier_32'])\nroc_auc = auc(fpr, tpr)\n\nplt.figure(figsize=(8,8))\nlw = 2\nplt.plot(fpr, tpr, color=BLUE,\n         lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)\nplt.plot([0, 1], [0, 1], color=RED, lw=lw, linestyle='--')\nplt.xlim([0.0, 1.0])\nplt.ylim([0.0, 1.05])\nplt.xlabel('False Positive Rate')\nplt.ylabel('True Positive Rate')\nplt.title('Classifier ROC curve')\nplt.legend(loc=\"lower right\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### LightGBM model"},{"metadata":{"trusted":true},"cell_type":"code","source":"features = [\n 'roll_2048_std_percentile_10',\n 'roll_4096_mean_percentile_95',\n 'bp6_percentile_90',\n 'bp2_interquantile_range_20_80',\n 'bp5_percentile_99',\n 'bp7_interquantile_range_10_90',\n 'bp4_percentile_80',\n 'roll_100_mean_mean',\n 'roll_1000_mean_percentile_99',\n 'roll_100_std_percentile_5',\n 'bp7_mad',\n 'evol_bp5_min',\n 'evol_bp7_max'\n]\n\nparams = {'num_leaves': 5,\n 'max_bin': 63,\n 'min_data_in_leaf': 20,\n 'learning_rate': 0.05,\n 'min_sum_hessian_in_leaf': 0.008,\n 'bagging_fraction': 0.9,\n 'bagging_freq': 1,\n 'feature_fraction': 1.,\n 'min_gain_to_split': 0.5,\n 'max_depth': 4,\n 'save_binary': True,\n 'seed': 0,\n 'feature_fraction_seed': 0,\n 'bagging_seed': 0,\n 'drop_seed': 0,\n 'data_random_seed': 0,\n 'objective': 'huber',\n 'boosting_type': 'gbdt',\n 'verbose': 1,\n 'metric': 'mae',\n 'boost_from_average': True,\n}\n\nnum_rounds = 2000\nearly_stopping = 50\nTARGET = 'time_to_failure'\n\nscores_train = []\nscores_val = []\ncount = 0\ntest['lgbm_prediction'] = 0\n\nfor i in range(K):\n    val_fold = folds[i]\n\n    cv_train = train[(~np.isin(train.fold, val_fold)) & (train.time_to_failure > 0.32)]\n    cv_val = train[np.isin(train.fold, val_fold) & (train.time_to_failure > 0.0)]\n\n    xg_train = lgb.Dataset(cv_train[features].values, label=cv_train[TARGET].values)\n    model = lgb.train(params, xg_train, num_rounds, verbose_eval=0)\n\n    train_preds = model.predict(cv_train[features].values, num_iteration=num_rounds)\n    val_preds = model.predict(cv_val[features].values, num_iteration=num_rounds)\n\n    train.loc[np.isin(train.fold, val_fold), 'lgbm_prediction'] = val_preds\n\n    scores_train.append(mae(cv_train[TARGET].values, train_preds))\n    scores_val.append(mae(cv_val[TARGET].values, val_preds))\n\n    count += 1\n    test['lgbm_prediction'] += model.predict(test[features].values, num_iteration=num_rounds)\n\ntest['lgbm_prediction'] = test['lgbm_prediction'] / count\n\nprint(\"MAE =\", mae(train.time_to_failure[:-5], train.lgbm_prediction[:-5]))\n\n# plotting\nfig, ax = plt.subplots(figsize=(18,6))\n\nplt.subplot(111)\nplt.plot(train[TARGET].values, c=RED, label='true')\nplt.plot(train.lgbm_prediction.values, c=BLUE, alpha=0.5, label='prediction')\nplt.title('LightGBM predictions')\nplt.xlabel('time')\nplt.ylabel('time_to_failure')\nplt.legend()\nplt.grid()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### NuCVR model"},{"metadata":{"trusted":true},"cell_type":"code","source":"feature = 'roll_100_std_percentile_5'\nscores_train = []\nscores_val = []\ncount = 0\ntest['nusvr_prediction'] = 0\nthreshold = 0.32\nfor i in range(K):\n    val_fold = folds[i]\n\n    cv_train = train[(~np.isin(train.fold, val_fold)) & (train.time_to_failure <= threshold)]\n    cv_val = train[np.isin(train.fold, val_fold) & (train.time_to_failure > 0.0)]\n    \n    feature_train = cv_train[feature].values\n    target_train = cv_train[TARGET].values\n    feature_val = cv_val[feature].values\n    target_val = cv_val[TARGET].values\n    feature_test = test[feature].values\n    \n    params = {\n     'C': 30,\n     'gamma': 0.05,\n     'kernel': 'rbf',\n     'nu': 1,\n    }\n\n    model = NuSVR(**params)\n    model.fit(feature_train.reshape(-1,1), 1/(threshold + 0.01 - target_train))\n\n    # evaluate quality\n    train_preds = threshold + 0.01 - 1/model.predict(feature_train.reshape(-1,1))\n    train_preds[train_preds < 0] = 0\n    val_preds = threshold + 0.01 - 1/model.predict(feature_val.reshape(-1,1))\n    val_preds[val_preds < 0] = 0\n    test_preds = threshold + 0.01 - 1/model.predict(feature_test.reshape(-1,1))\n    test_preds[test_preds < 0] = 0\n    \n    train.loc[np.isin(train.fold, val_fold), 'nusvr_prediction'] = val_preds\n\n    scores_train.append(mae(cv_train[TARGET].values, train_preds))\n    scores_val.append(mae(cv_val[TARGET].values, val_preds))\n    \n    count += 1\n    test['nusvr_prediction'] += test_preds\n\ntest['nusvr_prediction'] = test['nusvr_prediction'] / count\n\n# plotting\nfig, ax = plt.subplots(figsize=(18,6))\n\nplt.subplot(111)\nplt.plot(train[TARGET].values, c=RED, label='true')\nplt.plot(train.nusvr_prediction.values, c=BLUE, alpha=0.5, label='prediction')\nplt.title('NuSVR predictions')\nplt.xlabel('time')\nplt.ylabel('time_to_failure')\nplt.legend()\nplt.grid()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Mixing models\n\nLet's define the best threshold for classifier."},{"metadata":{"trusted":true},"cell_type":"code","source":"data = []\nfor x in np.linspace(0.,1.,21):\n    def mix(row):\n        t = x\n        if row.classifier_32 < t:\n            return row.lgbm_prediction\n        else:\n            return row.nusvr_prediction\n    train['final_prediction'] = train.apply(mix, axis=1)\n    data.append([np.round(x,2), mae(train.time_to_failure[:-5], train.final_prediction[:-5])])\npd.DataFrame(data=data, columns=['threshold', 'MAE'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def mix(row):\n    t = 0.85\n    if row.classifier_32 < t:\n        return row.lgbm_prediction\n    else:\n        return row.nusvr_prediction\n    \ntrain['final_prediction'] = train.apply(mix, axis=1)\ntest['time_to_failure'] = test.apply(mix, axis=1)\nprint(\"MAE =\", mae(train.time_to_failure[:-5], train.final_prediction[:-5]))\n\n# plotting\nfig, ax = plt.subplots(figsize=(18,6))\n\nplt.subplot(111)\nplt.plot(train[TARGET].values, c=RED, label='true')\nplt.plot(train.final_prediction.values, c=BLUE, alpha=0.5, label='prediction')\nplt.title('Final predictions')\nplt.xlabel('time')\nplt.ylabel('time_to_failure')\nplt.legend()\nplt.grid()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, we achieved a model with MAE=1.915."},{"metadata":{},"cell_type":"markdown","source":"#### Preparing submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv')\nsubmission = submission.drop('time_to_failure', axis=1)\npreds = test[['time_to_failure', 'seg_id']]\nsubmission = submission.merge(preds, how='left')\nsubmission.to_csv('submission_nusvr_lgbm_classifier.csv', index=False)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5. Small note on feature selection\n\nTo select features for models I used a very greedy forward feature selection method. To select features for LightGBM I used bigger learning rate to evaluate models faster.\n\n**Thanks for reading!!!**"}],"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}