{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## 1. Introduction\n<font color=#6698FF> Member 1: Carmen Chen (carmenxiuwenchen)\n    \n<font color=#6698FF>Member 2: Joeyin Lin (joeyinlin), \n\n<font color=#6698FF> Linear model:\n    \n<font color=#6698FF> Best Score: 2.50103\n    \n<font color=#6698FF>Best Public Score: 1.74125\n    \n<font color=#6698FF>Non Linear Model:\n\n<font color=#6698FF>Best Score: 2.45312\n    \n<font color=#6698FF>Best Public Score: 1.73448\n    \n<font color=#6698FF>Leaderboard rank: None since the competition closed 4 years ago","metadata":{"id":"3NnAgB-y0knj"}},{"cell_type":"markdown","source":"## 2. Data\n","metadata":{"id":"PB7FLdQn-dmK"}},{"cell_type":"code","source":"import gc\nimport os\nimport time\nimport logging\nimport datetime\nimport warnings\nimport numpy as np\nimport pandas as pd\npd.options.display.precision = 15\nfrom tqdm import tqdm, tqdm_notebook\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport xgboost as xgb\nimport lightgbm as lgb\n\nfrom scipy import stats\nfrom scipy.signal import hann\nfrom scipy.signal import hilbert\nfrom scipy.signal import convolve\nfrom sklearn.svm import NuSVR, SVR\nfrom catboost import CatBoostRegressor\nfrom sklearn.kernel_ridge import KernelRidge\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder\nfrom sklearn.metrics import confusion_matrix\nfrom sklearn.model_selection import KFold, StratifiedKFold, RepeatedKFold, train_test_split, cross_val_score\nfrom sklearn.linear_model import LogisticRegression, LinearRegression\nfrom sklearn import datasets, svm\nfrom sklearn.pipeline import Pipeline\nimport lightgbm as lgb\nwarnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:27:12.804271Z","iopub.execute_input":"2023-04-27T10:27:12.804710Z","iopub.status.idle":"2023-04-27T10:27:12.815210Z","shell.execute_reply.started":"2023-04-27T10:27:12.804660Z","shell.execute_reply":"2023-04-27T10:27:12.813835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Loading the data","metadata":{"id":"vfGikOxUAxJB"}},{"cell_type":"code","source":"train = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure':np.float64})","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:27:16.923301Z","iopub.execute_input":"2023-04-27T10:27:16.924615Z","iopub.status.idle":"2023-04-27T10:31:12.379797Z","shell.execute_reply.started":"2023-04-27T10:27:16.924550Z","shell.execute_reply":"2023-04-27T10:31:12.377572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:31:26.276027Z","iopub.execute_input":"2023-04-27T10:31:26.277308Z","iopub.status.idle":"2023-04-27T10:31:26.315470Z","shell.execute_reply.started":"2023-04-27T10:31:26.277239Z","shell.execute_reply":"2023-04-27T10:31:26.314229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:31:29.484338Z","iopub.execute_input":"2023-04-27T10:31:29.484773Z","iopub.status.idle":"2023-04-27T10:31:29.491928Z","shell.execute_reply.started":"2023-04-27T10:31:29.484737Z","shell.execute_reply":"2023-04-27T10:31:29.490937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Data Exploration","metadata":{}},{"cell_type":"code","source":"acoustic_data = train['acoustic_data'].values[::100] #samples for every 100 points of data\ntimefailure_data = train['time_to_failure'].values[::100] #samples for every 100 points of data\n\nfig, ax1 = plt.subplots(figsize=(12, 8))\nplt.title(\"Acoustic data and time to failure: 1% sampled data\")\nplt.plot(acoustic_data, color='r')\nax1.set_ylabel('acoustic data', color='r')\nplt.legend(['acoustic data'], loc=(0.01, 0.95))\nax2 = ax1.twinx()\nplt.plot(timefailure_data, color='b')\nax2.set_ylabel('time to failure', color='b')\nplt.legend(['time to failure'], loc=(0.01, 0.9))\nplt.grid(True)\n\ndel acoustic_data\ndel timefailure_data","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:31:32.370241Z","iopub.execute_input":"2023-04-27T10:31:32.371693Z","iopub.status.idle":"2023-04-27T10:31:35.743178Z","shell.execute_reply.started":"2023-04-27T10:31:32.371638Z","shell.execute_reply":"2023-04-27T10:31:35.742215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#here we plot the first 1% of the data to see what it looks like\n# 1 percent of 629145480 is approximately 6291455\nacoustic_data = train['acoustic_data'].values[:6291455] \ntimefailure_data = train['time_to_failure'].values[:6291455] \n\nfig, ax1 = plt.subplots(figsize=(12, 8))\nplt.title(\"Acoustic data and time to failure: first 1% of the data\")\nplt.plot(acoustic_data, color='r')\nax1.set_ylabel('acoustic data', color='r')\nplt.legend(['acoustic data'], loc=(0.01, 0.95))\nax2 = ax1.twinx()\nplt.plot(timefailure_data, color='b')\nax2.set_ylabel('time to failure', color='b')\nplt.legend(['time to failure'], loc=(0.01, 0.9))\nplt.grid(True)\n\ndel acoustic_data\ndel timefailure_data","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:31:58.732042Z","iopub.execute_input":"2023-04-27T10:31:58.732499Z","iopub.status.idle":"2023-04-27T10:32:01.341248Z","shell.execute_reply.started":"2023-04-27T10:31:58.732462Z","shell.execute_reply":"2023-04-27T10:32:01.340302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#6698FF> In the figure above we can see that the time to failure is slightly later than when the spike in the acoustic data occurs.","metadata":{}},{"cell_type":"markdown","source":"### Feature Extraction","metadata":{}},{"cell_type":"code","source":"#The test segments are 150000 each so we split the training data using the same dimension\nrows = 150000\nsegments = int(np.floor(train.shape[0] / rows))\nprint(segments)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:32:12.443938Z","iopub.execute_input":"2023-04-27T10:32:12.444460Z","iopub.status.idle":"2023-04-27T10:32:12.452674Z","shell.execute_reply.started":"2023-04-27T10:32:12.444414Z","shell.execute_reply":"2023-04-27T10:32:12.451130Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# From https://www.kaggle.com/code/artgor/even-more-features\ndef add_trend_feature(arr, abs_values=False):\n    idx = np.array(range(len(arr)))\n    if abs_values:\n        arr = np.abs(arr)\n    lr = LinearRegression()\n    lr.fit(idx.reshape(-1, 1), arr)\n    return lr.coef_[0]\n\ndef classic_sta_lta(x, length_sta, length_lta):\n    sta = np.cumsum(x ** 2)\n    # Convert to float\n    sta = np.require(sta, dtype=np.float)\n    # Copy for LTA\n    lta = sta.copy()\n    # Compute the STA and the LTA\n    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n    sta /= length_sta\n    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n    lta /= length_lta\n    # Pad zeros\n    sta[:length_lta - 1] = 0\n    # Avoid division by zero by setting zero values to tiny float\n    dtiny = np.finfo(0.0).tiny\n    idx = lta < dtiny\n    lta[idx] = dtiny\n    return sta / lta","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:32:15.674810Z","iopub.execute_input":"2023-04-27T10:32:15.675288Z","iopub.status.idle":"2023-04-27T10:32:15.686093Z","shell.execute_reply.started":"2023-04-27T10:32:15.675245Z","shell.execute_reply":"2023-04-27T10:32:15.684560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = pd.DataFrame(index=range(segments), dtype=np.float64)\ny_train = pd.DataFrame(index=range(segments), dtype=np.float64,\n                       columns=['time_to_failure'])","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:32:19.814303Z","iopub.execute_input":"2023-04-27T10:32:19.815252Z","iopub.status.idle":"2023-04-27T10:32:19.824092Z","shell.execute_reply.started":"2023-04-27T10:32:19.815197Z","shell.execute_reply":"2023-04-27T10:32:19.822693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tot_mean = train['acoustic_data'].mean()\ntot_std = train['acoustic_data'].std()\ntot_max = train['acoustic_data'].max()\ntot_min = train['acoustic_data'].min()\ntot_sum = train['acoustic_data'].sum()\ntotal_abs_sum = np.abs(train['acoustic_data']).sum()","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:32:23.209193Z","iopub.execute_input":"2023-04-27T10:32:23.209623Z","iopub.status.idle":"2023-04-27T10:32:32.147823Z","shell.execute_reply.started":"2023-04-27T10:32:23.209588Z","shell.execute_reply":"2023-04-27T10:32:32.146440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#We used the features from https://www.kaggle.com/code/gpreda/lanl-earthquake-eda-and-prediction\ndef create_features(seg_id, seg, X):\n    xc = pd.Series(seg['acoustic_data'].values)\n    zc = np.fft.fft(xc)\n    \n    X.loc[seg_id, 'mean'] = xc.mean()\n    X.loc[seg_id, 'std'] = xc.std()\n    X.loc[seg_id, 'max'] = xc.max()\n    X.loc[seg_id, 'min'] = xc.min()\n    \n    #FFT transform values\n    realFFT = np.real(zc)\n    imagFFT = np.imag(zc)\n    X.loc[seg_id, 'Rmean'] = realFFT.mean()\n    X.loc[seg_id, 'Rstd'] = realFFT.std()\n    X.loc[seg_id, 'Rmax'] = realFFT.max()\n    X.loc[seg_id, 'Rmin'] = realFFT.min()\n    X.loc[seg_id, 'Imean'] = imagFFT.mean()\n    X.loc[seg_id, 'Istd'] = imagFFT.std()\n    X.loc[seg_id, 'Imax'] = imagFFT.max()\n    X.loc[seg_id, 'Imin'] = imagFFT.min()\n    X.loc[seg_id, 'Rmean_last_5000'] = realFFT[-5000:].mean()\n    X.loc[seg_id, 'Rstd__last_5000'] = realFFT[-5000:].std()\n    X.loc[seg_id, 'Rmax_last_5000'] = realFFT[-5000:].max()\n    X.loc[seg_id, 'Rmin_last_5000'] = realFFT[-5000:].min()\n    X.loc[seg_id, 'Rmean_last_15000'] = realFFT[-15000:].mean()\n    X.loc[seg_id, 'Rstd_last_15000'] = realFFT[-15000:].std()\n    X.loc[seg_id, 'Rmax_last_15000'] = realFFT[-15000:].max()\n    X.loc[seg_id, 'Rmin_last_15000'] = realFFT[-15000:].min()\n    \n    X.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(xc))\n#     X.loc[seg_id, 'mean_change_rate'] = np.mean((np.diff(xc) / xc[:-1]))\n    X.loc[seg_id, 'abs_max'] = np.abs(xc).max()\n    X.loc[seg_id, 'abs_min'] = np.abs(xc).min()\n    \n    X.loc[seg_id, 'std_first_50000'] = xc[:50000].std()\n    X.loc[seg_id, 'std_last_50000'] = xc[-50000:].std()\n    X.loc[seg_id, 'std_first_10000'] = xc[:10000].std()\n    X.loc[seg_id, 'std_last_10000'] = xc[-10000:].std()\n    \n    X.loc[seg_id, 'avg_first_50000'] = xc[:50000].mean()\n    X.loc[seg_id, 'avg_last_50000'] = xc[-50000:].mean()\n    X.loc[seg_id, 'avg_first_10000'] = xc[:10000].mean()\n    X.loc[seg_id, 'avg_last_10000'] = xc[-10000:].mean()\n    \n    X.loc[seg_id, 'min_first_50000'] = xc[:50000].min()\n    X.loc[seg_id, 'min_last_50000'] = xc[-50000:].min()\n    X.loc[seg_id, 'min_first_10000'] = xc[:10000].min()\n    X.loc[seg_id, 'min_last_10000'] = xc[-10000:].min()\n    \n    X.loc[seg_id, 'max_first_50000'] = xc[:50000].max()\n    X.loc[seg_id, 'max_last_50000'] = xc[-50000:].max()\n    X.loc[seg_id, 'max_first_10000'] = xc[:10000].max()\n    X.loc[seg_id, 'max_last_10000'] = xc[-10000:].max()\n    \n    X.loc[seg_id, 'max_to_min'] = xc.max() / np.abs(xc.min())\n    X.loc[seg_id, 'max_to_min_diff'] = xc.max() - np.abs(xc.min())\n    X.loc[seg_id, 'count_big'] = len(xc[np.abs(xc) > 500])\n    X.loc[seg_id, 'sum'] = xc.sum()\n    \n#     X.loc[seg_id, 'mean_change_rate_first_50000'] = np.mean((np.diff(xc[:50000]) / xc[:50000][:-1]))\n#     X.loc[seg_id, 'mean_change_rate_last_50000'] = np.mean((np.diff(xc[-50000:]) / xc[-50000:][:-1]))\n#     X.loc[seg_id, 'mean_change_rate_first_10000'] = np.mean((np.diff(xc[:10000]) / xc[:10000][:-1]))\n#     X.loc[seg_id, 'mean_change_rate_last_10000'] = np.mean((np.diff(xc[-10000:]) / xc[-10000:][:-1]))\n    \n    X.loc[seg_id, 'q95'] = np.quantile(xc, 0.95)\n    X.loc[seg_id, 'q99'] = np.quantile(xc, 0.99)\n    X.loc[seg_id, 'q05'] = np.quantile(xc, 0.05)\n    X.loc[seg_id, 'q01'] = np.quantile(xc, 0.01)\n    \n    X.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(xc), 0.95)\n    X.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(xc), 0.99)\n    X.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(xc), 0.05)\n    X.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(xc), 0.01)\n    \n    X.loc[seg_id, 'trend'] = add_trend_feature(xc)\n    X.loc[seg_id, 'abs_trend'] = add_trend_feature(xc, abs_values=True)\n    X.loc[seg_id, 'abs_mean'] = np.abs(xc).mean()\n    X.loc[seg_id, 'abs_std'] = np.abs(xc).std()\n    \n    X.loc[seg_id, 'mad'] = xc.mad()\n    X.loc[seg_id, 'kurt'] = xc.kurtosis()\n    X.loc[seg_id, 'skew'] = xc.skew()\n    X.loc[seg_id, 'med'] = xc.median()\n    \n    X.loc[seg_id, 'Hilbert_mean'] = np.abs(hilbert(xc)).mean()\n    X.loc[seg_id, 'Hann_window_mean'] = (convolve(xc, hann(150), mode='same') / sum(hann(150))).mean()\n    X.loc[seg_id, 'classic_sta_lta1_mean'] = classic_sta_lta(xc, 500, 10000).mean()\n    X.loc[seg_id, 'classic_sta_lta2_mean'] = classic_sta_lta(xc, 5000, 100000).mean()\n    X.loc[seg_id, 'classic_sta_lta3_mean'] = classic_sta_lta(xc, 3333, 6666).mean()\n    X.loc[seg_id, 'classic_sta_lta4_mean'] = classic_sta_lta(xc, 10000, 25000).mean()\n    X.loc[seg_id, 'Moving_average_700_mean'] = xc.rolling(window=700).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_1500_mean'] = xc.rolling(window=1500).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_3000_mean'] = xc.rolling(window=3000).mean().mean(skipna=True)\n    X.loc[seg_id, 'Moving_average_6000_mean'] = xc.rolling(window=6000).mean().mean(skipna=True)\n    ewma = pd.Series.ewm\n    X.loc[seg_id, 'exp_Moving_average_300_mean'] = (ewma(xc, span=300).mean()).mean(skipna=True)\n    X.loc[seg_id, 'exp_Moving_average_3000_mean'] = ewma(xc, span=3000).mean().mean(skipna=True)\n    X.loc[seg_id, 'exp_Moving_average_30000_mean'] = ewma(xc, span=6000).mean().mean(skipna=True)\n    no_of_std = 2\n    X.loc[seg_id, 'MA_700MA_std_mean'] = xc.rolling(window=700).std().mean()\n    X.loc[seg_id,'MA_700MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X.loc[seg_id,'MA_700MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n    X.loc[seg_id, 'MA_400MA_std_mean'] = xc.rolling(window=400).std().mean()\n    X.loc[seg_id,'MA_400MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X.loc[seg_id,'MA_400MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n    X.loc[seg_id, 'MA_1000MA_std_mean'] = xc.rolling(window=1000).std().mean()\n    \n    X.loc[seg_id, 'iqr'] = np.subtract(*np.percentile(xc, [75, 25]))\n    X.loc[seg_id, 'q999'] = np.quantile(xc,0.999)\n    X.loc[seg_id, 'q001'] = np.quantile(xc,0.001)\n    X.loc[seg_id, 'ave10'] = stats.trim_mean(xc, 0.1)\n    \n    for windows in [10, 100, 1000]:\n        x_roll_std = xc.rolling(windows).std().dropna().values\n        x_roll_mean = xc.rolling(windows).mean().dropna().values\n        \n        X.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n#         X.loc[seg_id, 'av_change_rate_roll_std_' + str(windows)] = np.mean((np.diff(x_roll_std) / x_roll_std[:-1]))\n        X.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n#         X.loc[seg_id, 'av_change_rate_roll_mean_' + str(windows)] = np.mean((np.diff(x_roll_mean) / x_roll_mean[:-1]))\n        X.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:32:36.948734Z","iopub.execute_input":"2023-04-27T10:32:36.949183Z","iopub.status.idle":"2023-04-27T10:32:36.997651Z","shell.execute_reply.started":"2023-04-27T10:32:36.949144Z","shell.execute_reply":"2023-04-27T10:32:36.996315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# iterate over all segments\nfor seg_id in tqdm(range(segments)):\n    seg = train.iloc[seg_id*rows:seg_id*rows+rows]\n    create_features(seg_id, seg, X_train)\n    y_train.loc[seg_id, 'time_to_failure'] = seg['time_to_failure'].values[-1]","metadata":{"execution":{"iopub.status.busy":"2023-04-27T10:32:44.086511Z","iopub.execute_input":"2023-04-27T10:32:44.086936Z","iopub.status.idle":"2023-04-27T10:51:48.195859Z","shell.execute_reply.started":"2023-04-27T10:32:44.086902Z","shell.execute_reply":"2023-04-27T10:51:48.193696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train.shape\nX_train.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:29:59.662064Z","iopub.execute_input":"2023-04-27T11:29:59.663759Z","iopub.status.idle":"2023-04-27T11:29:59.738233Z","shell.execute_reply.started":"2023-04-27T11:29:59.663636Z","shell.execute_reply":"2023-04-27T11:29:59.736968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### Normalize the features","metadata":{}},{"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(X_train)\nscaled_X_train = pd.DataFrame(scaler.transform(X_train), columns=X_train.columns)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:03.657528Z","iopub.execute_input":"2023-04-27T11:30:03.658558Z","iopub.status.idle":"2023-04-27T11:30:03.697860Z","shell.execute_reply.started":"2023-04-27T11:30:03.658509Z","shell.execute_reply":"2023-04-27T11:30:03.696666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaled_X_train.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:09.870411Z","iopub.execute_input":"2023-04-27T11:30:09.871341Z","iopub.status.idle":"2023-04-27T11:30:09.899862Z","shell.execute_reply.started":"2023-04-27T11:30:09.871292Z","shell.execute_reply":"2023-04-27T11:30:09.898586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### PCA","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import PCA","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:16.382332Z","iopub.execute_input":"2023-04-27T11:30:16.382806Z","iopub.status.idle":"2023-04-27T11:30:16.492875Z","shell.execute_reply.started":"2023-04-27T11:30:16.382763Z","shell.execute_reply":"2023-04-27T11:30:16.491477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"principal_components = 10\npca = PCA(n_components=principal_components)\npca = pca.fit(scaled_X_train)\nX_pca_skl = pca.transform(scaled_X_train)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:18.319439Z","iopub.execute_input":"2023-04-27T11:30:18.319902Z","iopub.status.idle":"2023-04-27T11:30:18.429203Z","shell.execute_reply.started":"2023-04-27T11:30:18.319860Z","shell.execute_reply":"2023-04-27T11:30:18.427252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#plot \nplt.bar(range(0,principal_components), pca.explained_variance_ratio_, label=\"individual var\");\nplt.step(range(0,principal_components), np.cumsum(pca.explained_variance_ratio_),'r', label=\"cumulative var\");\nplt.xlabel('Principal component index'); plt.ylabel('explained variance ratio %');\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:21.384156Z","iopub.execute_input":"2023-04-27T11:30:21.384587Z","iopub.status.idle":"2023-04-27T11:30:21.673255Z","shell.execute_reply.started":"2023-04-27T11:30:21.384550Z","shell.execute_reply":"2023-04-27T11:30:21.672213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training and Results\n<font color=#6698FF> For our classification we used two different methods. We used one linear model and one non-linear model. For our linear model the support vector regression (SVR) method was used, while for our non-linear model Random Forest was used.</font>\n\n<font color=#6698FF> Before choosing our final linear and non-linear models, various other models were tested. This includes linear regression, K Neighbors and Deicision Tree.<font> \n    \n<font color=#6698FF> To assess the performance of our models, we looked at their score, MSE and MAE. Our Random Forest model had the highest score for our training data, while our SVR model had a higher score for our test data. It appears that our non-linear model overfits the data, which is why the score for training is much higher than its score for test.   ","metadata":{}},{"cell_type":"markdown","source":"##### Train-test split\n<font color=#6698FF> After some testing, we realised that our model performed better without using PCA than with the application of PCA, so in the end in train-test split we used `scaled_X_train` instead of `X_pca_skl`. We used a test size of 30% of the total data and set our random state to 102. Random state was set to 102 so we achieve the same results across different executions.</font>","metadata":{}},{"cell_type":"code","source":"x_train, x_test, Y_train, Y_test = train_test_split(scaled_X_train, y_train, test_size=0.3,shuffle=False, random_state=102)\n#Standard normalization\nscaler = StandardScaler().fit(x_train)\nnorm_x_train = scaler.transform(x_train)\nnorm_x_test = scaler.transform(x_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:31.930881Z","iopub.execute_input":"2023-04-27T11:30:31.931781Z","iopub.status.idle":"2023-04-27T11:30:31.961537Z","shell.execute_reply.started":"2023-04-27T11:30:31.931733Z","shell.execute_reply":"2023-04-27T11:30:31.960206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Linear models","metadata":{}},{"cell_type":"markdown","source":"#### Support Vector Regression (SVR)","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.svm import SVR\nlinmodel = SVR()\nlinmodel.fit(norm_x_train, Y_train.values.flatten())\n\n# Prediction\npredict_Y_train = linmodel.predict(norm_x_train)\npredict_Y_test = linmodel.predict(norm_x_test)\n\n# Accuracy score\ntrain_score = linmodel.score(norm_x_train, Y_train)\ntest_score = linmodel.score(norm_x_test, Y_test)\n\nprint(\"\\nTrain results\\nScore:\",train_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_train, predict_Y_train))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_train, predict_Y_train))\n\nprint(\"\\nTest results\\nScore:\",test_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_test, predict_Y_test))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_test, predict_Y_test))","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:35.438628Z","iopub.execute_input":"2023-04-27T11:30:35.439161Z","iopub.status.idle":"2023-04-27T11:30:42.115764Z","shell.execute_reply.started":"2023-04-27T11:30:35.439094Z","shell.execute_reply":"2023-04-27T11:30:42.114806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Prediction of train (linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the linear model (Train data)\")\nplt.plot(predict_Y_train, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y_train, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:50.818351Z","iopub.execute_input":"2023-04-27T11:30:50.818774Z","iopub.status.idle":"2023-04-27T11:30:51.258237Z","shell.execute_reply.started":"2023-04-27T11:30:50.818739Z","shell.execute_reply":"2023-04-27T11:30:51.256802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Prediction of test (linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\nt = np.linspace(0, 1259, 1259)\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the linear model (Test data)\")\nplt.plot(predict_Y_test, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.07, 0.9))\nax2 = ax1.twinx()\nplt.plot(t, Y_test, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.07, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:30:59.804507Z","iopub.execute_input":"2023-04-27T11:30:59.804952Z","iopub.status.idle":"2023-04-27T11:31:00.252099Z","shell.execute_reply.started":"2023-04-27T11:30:59.804914Z","shell.execute_reply":"2023-04-27T11:31:00.250599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Linear Regression","metadata":{}},{"cell_type":"code","source":"linmodel2 = LinearRegression()\nlinmodel2.fit(norm_x_train, Y_train.values.flatten())\n\n# Prediction\npredict_Y_train = linmodel2.predict(norm_x_train)\npredict_Y_test = linmodel2.predict(norm_x_test)\n\n# Accuracy score\ntrain_score = linmodel2.score(norm_x_train, Y_train)\ntest_score = linmodel2.score(norm_x_test, Y_test)\n\nprint(\"\\nTrain results\\nScore:\",train_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_train, predict_Y_train))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_train, predict_Y_train))\n\nprint(\"\\nTest results\\nScore:\",test_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_test, predict_Y_test))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_test, predict_Y_test))\n\n#Prediction of train (linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the linear model (Train data)\")\nplt.plot(predict_Y_train, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y_train, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)\n\n#Prediction of test (linear)\nfig, ax3 = plt.subplots(figsize=(20, 4))\nt = np.linspace(0, 1259, 1259)\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the linear model (Test data)\")\nplt.plot(predict_Y_test, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.07, 0.9))\nax4 = ax3.twinx()\nplt.plot(t, Y_test, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.07, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:32:04.362346Z","iopub.execute_input":"2023-04-27T11:32:04.363554Z","iopub.status.idle":"2023-04-27T11:32:05.372297Z","shell.execute_reply.started":"2023-04-27T11:32:04.363497Z","shell.execute_reply":"2023-04-27T11:32:05.371373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#6698FF> As can be seen in the figures above, the linear regression model has terrible test results and SVR performs much better, which is why SVR was chosen as our linear model in the end.</font>","metadata":{}},{"cell_type":"markdown","source":"### Nonlinear model","metadata":{}},{"cell_type":"markdown","source":"#### Random Forest","metadata":{}},{"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\nnonlinmodel = RandomForestRegressor(max_depth = 10, random_state=102)\nnonlinmodel.fit(norm_x_train, Y_train.values.flatten())\n\ntrain_pred = nonlinmodel.predict(norm_x_train)\ntest_pred = nonlinmodel.predict(norm_x_test)\n\ntrain_score = nonlinmodel.score(norm_x_train, Y_train)\ntest_score = nonlinmodel.score(norm_x_test, Y_test)\n\nprint(\"\\nTrain results\\nScore:\",train_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_train, train_pred))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_train, train_pred))\n\nprint(\"\\nTest results\\nScore:\",test_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_test, test_pred))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_test, test_pred))","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:33:53.364143Z","iopub.execute_input":"2023-04-27T11:33:53.365514Z","iopub.status.idle":"2023-04-27T11:34:12.290281Z","shell.execute_reply.started":"2023-04-27T11:33:53.365462Z","shell.execute_reply":"2023-04-27T11:34:12.288906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Prediction of train (linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the non-linear model (Train data)\")\nplt.plot(train_pred, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y_train, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:40:07.622024Z","iopub.execute_input":"2023-04-27T11:40:07.622615Z","iopub.status.idle":"2023-04-27T11:40:08.075157Z","shell.execute_reply.started":"2023-04-27T11:40:07.622567Z","shell.execute_reply":"2023-04-27T11:40:08.074155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction of test (non-linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\nt = np.linspace(0, 1259, 1259)\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the non-linear model (Test data)\")\nplt.plot(test_pred, color='b')\nax1.set_ylabel('time_to_failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(t, Y_test, color='r')\nax2.set_ylabel('time_to_failure', color='r')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.8))\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:40:13.965128Z","iopub.execute_input":"2023-04-27T11:40:13.965560Z","iopub.status.idle":"2023-04-27T11:40:14.394835Z","shell.execute_reply.started":"2023-04-27T11:40:13.965517Z","shell.execute_reply":"2023-04-27T11:40:14.393500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### Feature Importance","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(14,20))\nplt.barh(X_train.columns, nonlinmodel.feature_importances_)\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2023-04-27T12:39:58.212992Z","iopub.execute_input":"2023-04-27T12:39:58.213852Z","iopub.status.idle":"2023-04-27T12:40:01.737883Z","shell.execute_reply.started":"2023-04-27T12:39:58.213804Z","shell.execute_reply":"2023-04-27T12:40:01.736413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### K Neighbors","metadata":{}},{"cell_type":"code","source":"from sklearn.neighbors import KNeighborsRegressor\nknnmodel = KNeighborsRegressor()\nknnmodel.fit(norm_x_train, Y_train.values.flatten())\n\ntrain_pred = knnmodel.predict(norm_x_train)\ntest_pred = knnmodel.predict(norm_x_test)\n\ntrain_score = knnmodel.score(norm_x_train, Y_train)\ntest_score = knnmodel.score(norm_x_test, Y_test)\n\nprint(\"\\nTrain results\\nScore:\",train_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_train, train_pred))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_train, train_pred))\n\nprint(\"\\nTest results\\nScore:\",test_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_test, test_pred))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_test, test_pred))\n\n#Prediction of train (linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the non-linear model (Train data)\")\nplt.plot(train_pred, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y_train, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)\n\n# Prediction of test (non-linear)\nfig, ax3 = plt.subplots(figsize=(20, 4))\nt = np.linspace(0, 1259, 1259)\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the non-linear model (Test data)\")\nplt.plot(test_pred, color='b')\nax1.set_ylabel('time_to_failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax4 = ax3.twinx()\nplt.plot(t, Y_test, color='r')\nax2.set_ylabel('time_to_failure', color='r')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.8))\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T12:00:12.036792Z","iopub.execute_input":"2023-04-27T12:00:12.037317Z","iopub.status.idle":"2023-04-27T12:00:14.430167Z","shell.execute_reply.started":"2023-04-27T12:00:12.037276Z","shell.execute_reply":"2023-04-27T12:00:14.428711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Decision Tree","metadata":{}},{"cell_type":"code","source":"from sklearn.tree import DecisionTreeRegressor\ndtmodel = DecisionTreeRegressor(max_leaf_nodes = 30) #max leaf nodes was found empirically\ndtmodel.fit(norm_x_train, Y_train.values.flatten())\n\ntrain_pred = dtmodel.predict(norm_x_train)\ntest_pred = dtmodel.predict(norm_x_test)\n\ntrain_score = dtmodel.score(norm_x_train, Y_train)\ntest_score = dtmodel.score(norm_x_test, Y_test)\n\nprint(\"\\nTrain results\\nScore:\",train_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_train, train_pred))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_train, train_pred))\n\nprint(\"\\nTest results\\nScore:\",test_score)\nprint(\"Mean Squared Error (MSE):\",mean_squared_error(Y_test, test_pred))\nprint(\"Mean Absolute Error (MAE):\",mean_absolute_error(Y_test, test_pred))\n\n#Prediction of train (linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the non-linear model (Train data)\")\nplt.plot(train_pred, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y_train, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)\n\n# Prediction of test (non-linear)\nfig, ax3 = plt.subplots(figsize=(20, 4))\nt = np.linspace(0, 1259, 1259)\nplt.title(\"Comparison of the predicted time to failure vs. the actual time to failure using the non-linear model (Test data)\")\nplt.plot(test_pred, color='b')\nax1.set_ylabel('time_to_failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax4 = ax3.twinx()\nplt.plot(t, Y_test, color='r')\nax2.set_ylabel('time_to_failure', color='r')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.8))\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T12:05:24.516488Z","iopub.execute_input":"2023-04-27T12:05:24.518171Z","iopub.status.idle":"2023-04-27T12:05:25.606868Z","shell.execute_reply.started":"2023-04-27T12:05:24.518082Z","shell.execute_reply":"2023-04-27T12:05:25.605554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#6698FF>As can be seen above, Random Forest performs the best out of the three different non-linear models that we have tried. Overall, it seems that SVR performs the best out of all methods.</font>","metadata":{}},{"cell_type":"markdown","source":"### Process the Test Data","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(columns=X_train.columns, dtype=np.float64, index=submission.index)\nsubmission.shape, X_test.shape\n\n#Perform feature extraction\nfor seg_id in tqdm_notebook(X_test.index):\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv')\n    create_features(seg_id, seg, X_test)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:43:59.307611Z","iopub.execute_input":"2023-04-27T11:43:59.308171Z","iopub.status.idle":"2023-04-27T11:56:37.652485Z","shell.execute_reply.started":"2023-04-27T11:43:59.308130Z","shell.execute_reply":"2023-04-27T11:56:37.651015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_test.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T11:57:22.583865Z","iopub.execute_input":"2023-04-27T11:57:22.584841Z","iopub.status.idle":"2023-04-27T11:57:22.634092Z","shell.execute_reply.started":"2023-04-27T11:57:22.584779Z","shell.execute_reply":"2023-04-27T11:57:22.632810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaler = StandardScaler().fit(X_test)\nnorm_X_test = scaler.transform(X_test)\n# #Perform PCA\n# principal_components = 10\n# pca = PCA(n_components = principal_components)\n# pca = pca.fit(norm_X_test)\n# X_test_pca = pca.transform(X_test)\n# X_test_pca = X_test_pca[:,0:2]","metadata":{"execution":{"iopub.status.busy":"2023-04-27T12:16:46.913023Z","iopub.execute_input":"2023-04-27T12:16:46.914367Z","iopub.status.idle":"2023-04-27T12:16:46.930724Z","shell.execute_reply.started":"2023-04-27T12:16:46.914319Z","shell.execute_reply":"2023-04-27T12:16:46.929445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### Submission","metadata":{}},{"cell_type":"code","source":"# Submission\nsubmission['time_to_failure'] = nonlinmodel.predict(norm_X_test) \nsubmission.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-04-27T13:54:50.501249Z","iopub.execute_input":"2023-04-27T13:54:50.501769Z","iopub.status.idle":"2023-04-27T13:54:50.563222Z","shell.execute_reply.started":"2023-04-27T13:54:50.501731Z","shell.execute_reply":"2023-04-27T13:54:50.562025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T13:54:55.871497Z","iopub.execute_input":"2023-04-27T13:54:55.872289Z","iopub.status.idle":"2023-04-27T13:54:55.885730Z","shell.execute_reply.started":"2023-04-27T13:54:55.872242Z","shell.execute_reply":"2023-04-27T13:54:55.884306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Discussion and Conclusion\n\n<font color=#6698FF>The objective of this competition was to determine when the earthquake will take place. First we processed the data and then extracted features. The features we extracted included the basic statistical measures such as the mean, standard deviation, etc., FFT features and more. We performed this on the first and last 10000 and 50000 samples. After the feature extraction we performed PCA to reduce the dimensionality of our data. Later on, we realised that our models performed worse when PCA was implemented so in the end PCA was not applied.   </font>\n\n<font color=#6698FF> Before running our models, we split our data 70-30 using `train_test_split`. This method was mainly chosen because of its convenience, but in the future it may be better to use KFold cross validation instead as it may lead to higher accuracies.\n    \n<font color=#6698FF> We have tried 5 different models, namely SVR, linear regression, Random Forest, K neighbors, and Decision Tree. Their performance was assessed by looking at their score, MSE and MAE. It seemed that out of the five models random forest and SVR performed the \"best\". Random forest had the highest score for the training data, while SVR had the highest score for the test data. Random forest had a much lower score for its test data than training data which insinuates that the model was overfitted. In the future we could experiment with LightGBM,XGBoost and Catboost which are models used by the highest scoring teams.\n    \n<font color=#6698FF> Once we assessed the models we processed the test data. We submitted both our nonlinear (Random Forest) and linear model (SVR) to the competition as we were curious which model would perform better. Using the non linear model, we achieved a public score of 1.73448 and a private score of 2.45312. Using the linear model, we achieved a public score of 1.74125 and a private score of 2.50103.\n    \n<font color=#6698FF>All in all we are quite satisfied with this score. As said before in the future we may implement other methods such a k fold cross validation, LightGBM, Catboost, and XGBoost.","metadata":{}},{"cell_type":"markdown","source":"### References:\n* https://www.kaggle.com/code/gpreda/lanl-earthquake-eda-and-prediction\n* https://www.kaggle.com/code/inversion/basic-feature-benchmark/notebook\n* https://www.kaggle.com/code/oguzkoroglu/andrews-new-script-genetic-program-and-gplearn\n* https://www.kaggle.com/code/artgor/even-more-features\n* https://www.kaggle.com/code/vipul97/lanl-earthquake-prediction-model\n* https://www.kaggle.com/code/seneralkan/non-linear-regressions-models","metadata":{}}]}