{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":11000,"databundleVersionId":875412,"sourceType":"competition"}],"dockerImageVersionId":30458,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":">## By Mounzir Baroud and Dennis Landman\n>## Kaggle Username: Mounzir and Dennis Landman\n\n","metadata":{"id":"3NnAgB-y0knj"}},{"cell_type":"markdown","source":"## 1. Introduction","metadata":{}},{"cell_type":"markdown","source":"The LANL earthquake competition was an event that challenged participants to develop innovative solutions for predicting seismic activity. The competition was open to participants from around the world, offering a unique opportunity for students to apply their knowledge of machine learning and data analysis to a real-world problem.\nEventhough the competition has ended, it is still possible to submit solutions to see how they would perform.\n\nOur best effort would have landed us around the 1800 mark with a score of 2.59383.\n","metadata":{}},{"cell_type":"markdown","source":"## 2. Data Analysis & Feature Extraction\n","metadata":{"id":"PB7FLdQn-dmK"}},{"cell_type":"markdown","source":"### Load packages","metadata":{}},{"cell_type":"code","source":"from tqdm.notebook import tqdm\nimport numpy as np\nimport pandas as pd\npd.options.display.precision = 15\nimport os\nprint(os.listdir(\"../input\"))\nimport matplotlib.pyplot as plt\nfrom sklearn.linear_model import LinearRegression\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.svm import NuSVR, SVR\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.ensemble import RandomForestRegressor as RFR\nfrom sklearn.model_selection import train_test_split, KFold\nimport lightgbm ","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:49:20.813873Z","iopub.execute_input":"2024-07-06T13:49:20.814322Z","iopub.status.idle":"2024-07-06T13:49:20.823489Z","shell.execute_reply.started":"2024-07-06T13:49:20.814279Z","shell.execute_reply":"2024-07-06T13:49:20.822297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.1.1 Data exploration","metadata":{}},{"cell_type":"code","source":"# Load training data\ntrain = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:49:20.825962Z","iopub.execute_input":"2024-07-06T13:49:20.826326Z","iopub.status.idle":"2024-07-06T13:52:15.162972Z","shell.execute_reply.started":"2024-07-06T13:49:20.826289Z","shell.execute_reply":"2024-07-06T13:52:15.160529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the figure it is clear that all the spikes in the time_to_failure (green line) correspond to the spikes in the acoustic_data (blue line).","metadata":{}},{"cell_type":"code","source":"# Plot the whole training data to get a visual insight\n\ntrain_acoustic_data_small = train['acoustic_data'].values[::50]\ntrain_time_to_failure_small = train['time_to_failure'].values[::50]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel train_acoustic_data_small\ndel train_time_to_failure_small","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:52:15.166018Z","iopub.execute_input":"2024-07-06T13:52:15.167241Z","iopub.status.idle":"2024-07-06T13:52:29.828335Z","shell.execute_reply.started":"2024-07-06T13:52:15.167181Z","shell.execute_reply":"2024-07-06T13:52:29.826956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A closer in look at the first earth quake sees a delay between the spike of the acoustic_data and the actual moment of failure.","metadata":{}},{"cell_type":"code","source":"# Plotting a small portion of the full data to get a more zoomed in look at the first quake\ntrain_acoustic_data_small = train['acoustic_data'].values[:8091455]\ntrain_time_to_failure_small = train['time_to_failure'].values[:8091455]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure.\")\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel train_acoustic_data_small\ndel train_time_to_failure_small","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:52:29.830300Z","iopub.execute_input":"2024-07-06T13:52:29.831110Z","iopub.status.idle":"2024-07-06T13:52:38.775973Z","shell.execute_reply.started":"2024-07-06T13:52:29.831057Z","shell.execute_reply":"2024-07-06T13:52:38.774783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Checking for NaNs\ntrain.isna().sum()","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:52:38.779082Z","iopub.execute_input":"2024-07-06T13:52:38.779577Z","iopub.status.idle":"2024-07-06T13:52:40.966319Z","shell.execute_reply.started":"2024-07-06T13:52:38.779537Z","shell.execute_reply":"2024-07-06T13:52:40.965072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.1.2 Feature Extraction & Data Processing","metadata":{}},{"cell_type":"code","source":"# Splitting the data in smaller sections\nrows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:52:40.967808Z","iopub.execute_input":"2024-07-06T13:52:40.968238Z","iopub.status.idle":"2024-07-06T13:52:40.973983Z","shell.execute_reply.started":"2024-07-06T13:52:40.968198Z","shell.execute_reply":"2024-07-06T13:52:40.972621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extracting features function\ndef extract_features_from_segment(segment, seg, X):\n    x = pd.Series(seg['acoustic_data'].values)\n\n    X.loc[segment, 'mean'] = x.mean()\n    X.loc[segment, 'std'] = x.std()\n    X.loc[segment, 'max'] = x.max()\n    X.loc[segment, 'min'] = x.min()\n    X.loc[segment, 'kurt'] = x.kurtosis()\n    X.loc[segment, 'skew'] = x.skew()\n    X.loc[segment, 'q99'] = np.quantile(x, 0.99)\n    X.loc[segment, 'q95'] = np.quantile(x, 0.95)\n    X.loc[segment, 'q05'] = np.quantile(x, 0.05)\n    X.loc[segment, 'q01'] = np.quantile(x, 0.01)\n    \n    \n#Rolling window of size 50\n    x_roll_std = x.rolling(window = 50).std().dropna().values\n    x_roll_mean = x.rolling(window= 50).mean().dropna().values\n    \n    #Features from rolling mean\n    X.loc[segment, 'rolling_mean'] =x_roll_mean.mean()\n    X.loc[segment, 'rolling_mean'] =x_roll_mean.std()\n    X.loc[segment, 'rolling_mean'] =x_roll_mean.min()\n    X.loc[segment, 'rolling_mean'] =x_roll_mean.max()\n    X.loc[segment, 'rolling_mean_q99'] = np.quantile(x_roll_mean, 0.99)\n    X.loc[segment, 'rolling_mean_q01'] = np.quantile(x_roll_mean, 0.01)\n    \n    #Features from rolling std\n    X.loc[segment, 'rolling_std'] = x_roll_std.mean()\n    X.loc[segment, 'rolling_std'] = x_roll_std.std()\n    X.loc[segment, 'rolling_std'] = x_roll_std.min()\n    X.loc[segment, 'rolling_std'] = x_roll_std.max()\n    X.loc[segment, 'rolling_std_q99'] = np.quantile(x_roll_std, 0.99)\n    X.loc[segment, 'rolling_std_q01'] = np.quantile(x_roll_std, 0.01)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-07-06T13:52:40.976062Z","iopub.execute_input":"2024-07-06T13:52:40.976766Z","iopub.status.idle":"2024-07-06T13:52:41.219159Z","shell.execute_reply.started":"2024-07-06T13:52:40.976710Z","shell.execute_reply":"2024-07-06T13:52:41.217630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Extracting features on training data\n\nX_train = pd.DataFrame(index=range(segments), dtype=np.float64)\ny_train = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\nfor segment in tqdm(range(segments)):\n    seg = train.iloc[segment*rows:segment*rows+rows]\n    extract_features_from_segment(segment, seg, X_train)\n    y_train.loc[segment, 'time_to_failure'] = seg['time_to_failure'].values[-1]","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:52:41.220841Z","iopub.execute_input":"2024-07-06T13:52:41.221475Z","iopub.status.idle":"2024-07-06T13:54:54.747298Z","shell.execute_reply.started":"2024-07-06T13:52:41.221412Z","shell.execute_reply":"2024-07-06T13:54:54.745715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Small peak at the features that are extracted\nX_train.head(10)","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:54:54.749271Z","iopub.execute_input":"2024-07-06T13:54:54.749784Z","iopub.status.idle":"2024-07-06T13:54:54.811103Z","shell.execute_reply.started":"2024-07-06T13:54:54.749731Z","shell.execute_reply":"2024-07-06T13:54:54.809631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Normalizing data\nscaler = StandardScaler()\nscaler.fit(X_train)\nX_train_scaled = scaler.transform(X_train)","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:54:54.812815Z","iopub.execute_input":"2024-07-06T13:54:54.813286Z","iopub.status.idle":"2024-07-06T13:54:54.828996Z","shell.execute_reply.started":"2024-07-06T13:54:54.813240Z","shell.execute_reply":"2024-07-06T13:54:54.827548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Models","metadata":{}},{"cell_type":"markdown","source":"The assignment required the use of both a linear and a non-linear model. For the linear model we chose to use LinearRegression due to its simplicity. As for the non-linear model we experimented with different options like Support Vector Regression (SVR), Random Forest Regression (RFR) and Light Gradient Boosted Machine Regression (LGBMR) with Cross Validation. These models offer more flexibility to conform to the real world dataset, but can be prone to overfitting.","metadata":{}},{"cell_type":"markdown","source":"> ## 3.1 Linear","metadata":{}},{"cell_type":"markdown","source":"### 3.1.1 LinearRegression","metadata":{}},{"cell_type":"code","source":"# Model #1 - LinearRegression\n\n# define the model\nmodel_lr = LinearRegression()\n# fit the model\nmodel_lr.fit(X_train_scaled, y_train)\ny_pred_lin = model_lr.predict(X_train_scaled)\n\nscore = mean_absolute_error(y_train.values.flatten(), y_pred_lin)\nprint(f'Score: {score:0.3f}')","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:54:54.830994Z","iopub.execute_input":"2024-07-06T13:54:54.831814Z","iopub.status.idle":"2024-07-06T13:54:54.856539Z","shell.execute_reply.started":"2024-07-06T13:54:54.831761Z","shell.execute_reply":"2024-07-06T13:54:54.855196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> ## 3.2 non-Linear","metadata":{}},{"cell_type":"markdown","source":"### 3.2.1 Nu Support Vector Machine","metadata":{}},{"cell_type":"code","source":"#Model #2 - Nu Support Vector Machine \n\nmodel_svm = NuSVR()\nmodel_svm.fit(X_train_scaled, y_train.values.flatten())\ny_pred = model_svm.predict(X_train_scaled)\n\nscore = mean_absolute_error(y_train.values.flatten(), y_pred)\nprint(f'Score: {score:0.3f}')","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:54:54.858250Z","iopub.execute_input":"2024-07-06T13:54:54.858678Z","iopub.status.idle":"2024-07-06T13:54:56.166761Z","shell.execute_reply.started":"2024-07-06T13:54:54.858635Z","shell.execute_reply":"2024-07-06T13:54:56.165399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.2.2 Random Forest Regression","metadata":{}},{"cell_type":"code","source":"# Model #3 - Random Forest Regression\n\nmodel_rfr = RFR(n_estimators = 100, random_state = 42)\nmodel_rfr.fit(X_train_scaled, y_train.values.flatten())\ny_pred2 = model_rfr.predict(X_train_scaled)\n\nscore = mean_absolute_error(y_train.values.flatten(), y_pred2)\nprint(f'Score: {score:0.3f}')","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:54:56.168553Z","iopub.execute_input":"2024-07-06T13:54:56.168942Z","iopub.status.idle":"2024-07-06T13:55:00.515449Z","shell.execute_reply.started":"2024-07-06T13:54:56.168904Z","shell.execute_reply":"2024-07-06T13:55:00.514217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.2.3 Light Gradient Boosting","metadata":{}},{"cell_type":"code","source":"# Model #4 - Light Gradient Boosting using Cross Validation\n\ncv = KFold(n_splits=5, random_state=1, shuffle=True)\n\nparam = {\n    'num_leaves': 37,\n    'objective': 'regression',\n    'max_depth':9,\n    'learning_rate': 0.01,\n    'boosting': 'gbdt',\n    'feature_fraction': 0.311473528266664,\n    'bagging_freq': 14,\n    'bagging_fraction': 0.4190138426177889,\n    'bagging_seed': 42,\n    'metric': 'mae',\n    'lambda_l1': 0.0003661365327691201,\n    'lambda_l2': 0.2855479018074398,\n    'verbosity': -1,\n    'nthread': -1,\n    'random_state': 0,\n    'min_data_in_leaf': 40\n}\nfor fold, (train_index, test_index) in enumerate(cv.split(X_train)):\n    \n    X_traincv, X_validcv = X_train.iloc[train_index], X_train.iloc[test_index]\n    y_traincv, y_validcv = y_train.iloc[train_index], y_train.iloc[test_index]\n    \n    train_data = lightgbm.Dataset(X_traincv, label=y_traincv)\n    valid_data = lightgbm.Dataset(X_validcv, label=y_validcv)\n    \n    model_lgbm = lightgbm.train(param,\n                            train_data,\n                            valid_sets=valid_data,\n                            num_boost_round=5000,\n                            early_stopping_rounds=50)","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2024-07-06T13:55:00.520183Z","iopub.execute_input":"2024-07-06T13:55:00.520782Z","iopub.status.idle":"2024-07-06T13:55:10.206699Z","shell.execute_reply.started":"2024-07-06T13:55:00.520744Z","shell.execute_reply":"2024-07-06T13:55:10.205388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Results for train\ny_train_pred = model_lgbm.predict(X_traincv)\ny_valid_pred = model_lgbm.predict(X_validcv)\nscore = mean_absolute_error(y_traincv.values.flatten(), y_train_pred)\nprint(f'Score: {score:0.3f}')","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:55:10.208400Z","iopub.execute_input":"2024-07-06T13:55:10.208817Z","iopub.status.idle":"2024-07-06T13:55:10.293010Z","shell.execute_reply.started":"2024-07-06T13:55:10.208777Z","shell.execute_reply":"2024-07-06T13:55:10.291626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Results for validation\nscorevalid = mean_absolute_error(y_validcv.values.flatten(), y_valid_pred)\nprint(f'Scorevalid: {scorevalid:0.3f}')\n","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:55:10.294469Z","iopub.execute_input":"2024-07-06T13:55:10.294812Z","iopub.status.idle":"2024-07-06T13:55:10.302085Z","shell.execute_reply.started":"2024-07-06T13:55:10.294778Z","shell.execute_reply":"2024-07-06T13:55:10.300930Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n","metadata":{}},{"cell_type":"markdown","source":"## 4. Test Data Feature Extraction & Submission","metadata":{}},{"cell_type":"markdown","source":"> ### 4.1 Test Data Feature Extraction","metadata":{}},{"cell_type":"code","source":"#Getting features from test data\n\nsubmission = pd.read_csv('/kaggle/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)\nfor seg_id in tqdm(X_test.index):\n    seg = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv')\n    extract_features_from_segment(seg_id, seg, X_test)","metadata":{"execution":{"iopub.status.busy":"2024-07-06T13:55:36.890037Z","iopub.execute_input":"2024-07-06T13:55:36.890509Z","iopub.status.idle":"2024-07-06T13:57:43.480958Z","shell.execute_reply.started":"2024-07-06T13:55:36.890466Z","shell.execute_reply":"2024-07-06T13:57:43.479629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> ### 4.2 Prediction & Submission","metadata":{}},{"cell_type":"code","source":"#Scaling test data and saving to the submission file using desired model\n\nmodels = {\"lr\": model_lr, \"svm\": model_svm, \"rfr\": model_rfr, \"lgbm\": model_lgbm}\nX_test_scaled = scaler.transform(X_test)\n\nfor label, model in models.items():\n    submission['time_to_failure'] = model.predict(X_test_scaled)\n    submission.to_csv(f'{label}_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2024-07-06T14:18:19.080988Z","iopub.execute_input":"2024-07-06T14:18:19.082271Z","iopub.status.idle":"2024-07-06T14:18:19.603967Z","shell.execute_reply.started":"2024-07-06T14:18:19.082206Z","shell.execute_reply":"2024-07-06T14:18:19.602814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5. Results","metadata":{}},{"cell_type":"markdown","source":"\n\n| Model | Internal Score | Public Score | Private Score |\n|-------|----------------|--------------|---------------|\n| Linear Regression | 2.167 | 1.81016 | 2.59383 |\n| SVM | 2.051 | 1.55187 | 2.59515 |\n| Random Forest Regression | 0.797 | 1.67659 | 2.59729 |\n| Light Gradient Boost | 1.852 | 5.42115 | 4.14112 |\n","metadata":{}},{"cell_type":"markdown","source":"## 6. Discussion and Conclusion","metadata":{}},{"cell_type":"markdown","source":"We first started with a LinearRegression model and a handfull of features to get a baseline. With this we already got reasonable scores, however we were not satisfied. We increased the complexities of the models used and added more features. In some instances this significantly improved our results, however most did not produce a noteworthy boost in score.\n\nAfter submission we observe a considerable decrease in all scores. Particularly the amazing score achieved by the Random Forest Regression model dramatically reduced. Similarly the Light Gradient Boost Machine Regression model performance degraded to the point of its score becoming insignificant compared to the other models. We suspect this to be a product of overfitting on the training data.\n\nTo further improve our rank we should consider fine-tuning hyperparameter to avoid overfitting in the more complex non-linear models. Furthermore adding more features and subsequently filtering for importance should further improve the performance of all our models.","metadata":{}},{"cell_type":"markdown","source":"## 7. References","metadata":{}},{"cell_type":"markdown","source":"[1] LANL Earthquake EDA and Prediction https://www.kaggle.com/code/gpreda/lanl-earthquake-eda-and-prediction#Prepare-the-data-analysis\n\n[2] Shaking Earth https://www.kaggle.com/code/allunia/shaking-earth#Rolling-features\n\n[3] Earthquakes FE. More features and samples https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples#Building-models\n\n[4] 3rd place memo - Discussion https://www.kaggle.com/competitions/LANL-Earthquake-Prediction/discussion/94459#549323\n\n[5] Baseline with multiple models https://www.kaggle.com/code/jsaguiar/baseline-with-multiple-models#1.-About-this-notebook","metadata":{}}]}