{"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":"# Final Project\n\nAs in all machine learning problems you should complete the following steps:\n1. Load and explore the data (using plots and histograms )\n2. Clean/preprocess/transform the data if necessary (first performed on training data and next on test data)\n3. Train the machine learning model (in this assignment two models, one linear and one non-linear)\n4. Evaluate and optimise the model\n\nHowever, while performing the steps you should consider our main goals and research questions for doing this project and try to adress them. This can be either in each step or afterwards in discussion and conclusion section (your call!). Below are the research question we are interested in:\n\n* What are the necessary step to clean and prepare the dataset as we have categorial features and missing data.\n* What model/classifier provides the best result for this application.\n* What is the best choice for our cost function and performance metrics for this problem.\n\nPlease also note the following points during the assignment:\n\n* Use functions from open source libraries like sci-kit learn and keras and avoid using your own hand-written functions from previous assignments. You can also use our [cheatsheet](https://colab.research.google.com/drive/12h-QBlsaWXkjGIRoJXfF4yqi1elnX9qn?usp=sharing).\n* Feel free to contact us on Teams if you need more description.\n* This notebook is structured like a scientific paper. The text should provide a high-level overview of your approach. Please don't include any details about your code in the text but add them as comments in the code itself. Your code should be cleane and readable with enough comments.\n* There are some instructions and questions in each section, remove the highlighted text in blue and replace it with your explanations and answers.\n\nYou can delete this section before submission.","metadata":{"id":"GF4WsCxg-lL9"}},{"cell_type":"markdown","source":"## 1. Introduction\n","metadata":{"id":"3NnAgB-y0knj"}},{"cell_type":"markdown","source":"<font color=#6698FF> Names: Sabita Kalila Putri Kurnia, Gordon Myburg\n    \n<font color=#6698FF> Kaggle Username: SabitaKalila, Gordon_Myburg\n    \n<font color=#6698FF> Score: \n    \n<font color=#6698FF> Leaderboard Rank:","metadata":{"id":"DPJV5e8sDVnT"}},{"cell_type":"markdown","source":"## 2. Data\n","metadata":{"id":"PB7FLdQn-dmK"}},{"cell_type":"markdown","source":"<font color=#6698FF> The needed libraries","metadata":{}},{"cell_type":"code","source":"from sklearn.linear_model import LinearRegression\nimport numpy as np\nimport pandas as pd\nimport gc\nimport os\nimport time\nimport logging\nimport datetime\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport xgboost as xgb\nimport lightgbm as lgb\nfrom scipy import stats\nfrom scipy.signal import hann\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\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.preprocessing import LabelEncoder\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\nwarnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2024-04-21T06:47:39.869267Z","iopub.execute_input":"2024-04-21T06:47:39.870149Z","iopub.status.idle":"2024-04-21T06:47:42.863256Z","shell.execute_reply.started":"2024-04-21T06:47:39.870100Z","shell.execute_reply":"2024-04-21T06:47:42.862292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.1 Dataset\n<font color=#6698FF>The data set is read and loaded here.","metadata":{"id":"vfGikOxUAxJB"}},{"cell_type":"code","source":"train_data = pd.read_csv(\"../input/LANL-Earthquake-Prediction/train.csv\",dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\ntrain_data.head(10)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T06:47:46.188562Z","iopub.execute_input":"2024-04-21T06:47:46.188980Z","iopub.status.idle":"2024-04-21T06:52:30.529068Z","shell.execute_reply.started":"2024-04-21T06:47:46.188941Z","shell.execute_reply":"2024-04-21T06:52:30.527466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"###### 2.1.1 Train-test split\n<font color=#6698FF>In the below, we split the train data into a test and a train set.","metadata":{"id":"udqAFg1C3JHe"}},{"cell_type":"code","source":"#Reference: jessteijn's and artgor's Kaggle Notebooks\n\n#artgor: https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples\n\ndef calc_change_rate(x):\n    change = (np.diff(x) / x[:-1]).values\n    change = change[np.nonzero(change)[0]]\n    change = change[~np.isnan(change)]\n    change = change[change != -np.inf]\n    change = change[change != np.inf]\n    return np.mean(change)\n\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\n#jesseijn = https://www.kaggle.com/code/jessteijn/ml23-week8\n\ndef x_features(x_data, X, segment):\n    X.loc[segment, 'ave'] = x_data.mean()\n    X.loc[segment, 'std'] = x_data.std()\n    X.loc[segment, 'max'] = x_data.max()\n    X.loc[segment, 'min'] = x_data.min()\n    \n    X.loc[segment, 'q01'] = np.quantile(x_data,0.01)\n    X.loc[segment, 'q05'] = np.quantile(x_data,0.05)\n    X.loc[segment, 'q95'] = np.quantile(x_data,0.95)\n    X.loc[segment, 'q99'] = np.quantile(x_data,0.99)\n    \n    X.loc[segment, 'mean_change_abs'] = np.mean(np.diff(x_data))\n    X.loc[segment, 'mean_change_rate'] = calc_change_rate(x_data)\n    X.loc[segment, 'abs_max'] = np.abs(x_data).max()\n    X.loc[segment, 'abs_min'] = np.abs(x_data).min()\n\n   \n    X.loc[segment, 'max_to_min'] = x_data.max() / np.abs(x_data.min())\n    X.loc[segment, 'max_to_min_diff'] = x_data.max() - np.abs(x_data.min())\n    X.loc[segment, 'count_big'] = len(x_data[np.abs(x_data) > 500])\n    X.loc[segment, 'sum'] = x_data.sum()\n\n    X.loc[segment, 'abs_trend'] = add_trend_feature(x_data, abs_values=True)\n    X.loc[segment, 'abs_mean'] = np.abs(x_data).mean()\n    X.loc[segment, 'abs_std'] = np.abs(x_data).std()\n    X.loc[segment, 'abs_median'] = np.median(np.abs(x_data))\n\n    X.loc[segment, 'trend'] = add_trend_feature(x_data)\n    X.loc[segment, 'mad'] = x_data.mad()\n    X.loc[segment, 'kurt'] = x_data.kurtosis()\n    X.loc[segment, 'skew'] = x_data.skew()\n    X.loc[segment, 'med'] = x_data.median()\n    \n    X.loc[segment, 'abs_q95'] = np.quantile(np.abs(x_data),0.95)\n    X.loc[segment, 'abs_q99'] = np.quantile(np.abs(x_data),0.99)\n    X.loc[segment, 'abs_q05'] = np.quantile(np.abs(x_data), 0.05)\n    X.loc[segment, 'abs_q01'] = np.quantile(np.abs(x_data), 0.01)\n    X.loc[segment, 'F_test'], X.loc[segment, 'p_test'] = stats.f_oneway(x_data[:40000],x_data[40000:80000],x_data[80000:100000],x_data[100000:130000],x_data[130000:])\n    X.loc[segment, 'av_change_abs'] = np.mean(np.diff(x_data))\n        \n    return X\n\ndef extract_features(train_data,segments):\n    X = pd.DataFrame(index=range(segments), dtype=np.float64)\n    y = pd.DataFrame(index=range(segments), dtype=np.float64,columns=['time_to_failure'])\n    for segment in tqdm(range(segments)):\n        seg = train_data.iloc[segment*rows:segment*rows+rows]\n        x_data = pd.Series(seg['acoustic_data'].values)\n        y_data = seg['time_to_failure'].values[-1]\n\n        y.loc[segment,'time_to_failure'] = y_data\n        \n        X = x_features(x_data, X, segment)\n        \n    return X, y","metadata":{"execution":{"iopub.status.busy":"2024-04-21T06:52:33.902377Z","iopub.execute_input":"2024-04-21T06:52:33.902810Z","iopub.status.idle":"2024-04-21T06:52:33.932828Z","shell.execute_reply.started":"2024-04-21T06:52:33.902765Z","shell.execute_reply":"2024-04-21T06:52:33.931240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rows = 150_000\nsegments = int(np.floor(train_data.shape[0] / rows))\n\nX,y = extract_features(train_data, segments)\n\ny=y.values.flatten()\nprint(X.shape)\nprint(y.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T06:52:37.668200Z","iopub.execute_input":"2024-04-21T06:52:37.669032Z","iopub.status.idle":"2024-04-21T07:00:58.172982Z","shell.execute_reply.started":"2024-04-21T06:52:37.668985Z","shell.execute_reply":"2024-04-21T07:00:58.171657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Splitting the data into test and train set\n\nfrom sklearn.model_selection import train_test_split\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,shuffle=True, random_state=100)\n\n#checking the shape of the training and test data\n\nprint(X_train.shape)\nprint(X_test.shape)\nprint(y_train.shape)\nprint(y_test.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:02.037864Z","iopub.execute_input":"2024-04-21T07:01:02.038341Z","iopub.status.idle":"2024-04-21T07:01:02.062921Z","shell.execute_reply.started":"2024-04-21T07:01:02.038299Z","shell.execute_reply":"2024-04-21T07:01:02.061300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Data Exploration\n\n<font color=#6698FF> Here, the data that were uploaded are explored by plotting them and see how the data behaves.","metadata":{"id":"QN5jS2_p_Y9N"}},{"cell_type":"code","source":"#Plotting the 1% of the data to have a better overview. Here we sample every 50 points of the data.\n\n#Reference: gpreda's Kaggle\n#https://www.kaggle.com/code/gpreda/lanl-earthquake-eda-and-prediction#Prepare-the-data-analysis\n\nimport matplotlib.pyplot as plt\n\ntrain_ad_sample_df = train_data['acoustic_data'].values[::50] \ntrain_ttf_sample_df = train_data['time_to_failure'].values[::50]\n\ndef plot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% sampled data\"):\n    fig, ax1 = plt.subplots(figsize=(12, 8))\n    plt.title(title)\n    plt.plot(train_ad_sample_df, color='r')\n    ax1.set_ylabel('acoustic data', color='r')\n    plt.legend(['acoustic data'], loc=(0.01, 0.95))\n    ax2 = ax1.twinx()\n    plt.plot(train_ttf_sample_df, color='b')\n    ax2.set_ylabel('time to failure', color='b')\n    plt.legend(['time to failure'], loc=(0.01, 0.9))\n    plt.grid(True)\n\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df)\ndel train_ad_sample_df\ndel train_ttf_sample_df","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:04.739973Z","iopub.execute_input":"2024-04-21T07:01:04.741173Z","iopub.status.idle":"2024-04-21T07:01:10.379726Z","shell.execute_reply.started":"2024-04-21T07:01:04.741118Z","shell.execute_reply":"2024-04-21T07:01:10.378200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#6698FF> From the graph we can see that acoustic data peaks before time to failure. Also the largest amplitude of the acoustic data is when time to failure is zero.","metadata":{}},{"cell_type":"code","source":"#Zoomed in portion of the first quake\n\n# Reference: gpreda's Kaggle\n#https://www.kaggle.com/code/gpreda/lanl-earthquake-eda-and-prediction#Prepare-the-data-analysis\n\n\ntrain_ad_sample_df = train_data['acoustic_data'].values[:6291455]\ntrain_ttf_sample_df = train_data['time_to_failure'].values[:6291455]\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% of data\")\ndel train_ad_sample_df\ndel train_ttf_sample_df","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:13.957576Z","iopub.execute_input":"2024-04-21T07:01:13.961144Z","iopub.status.idle":"2024-04-21T07:01:19.182308Z","shell.execute_reply.started":"2024-04-21T07:01:13.960847Z","shell.execute_reply":"2024-04-21T07:01:19.180990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#6698FF> From the zoomed in plot, we can see that there are small oscillations before the big one in the acoustic data. After the largest oscillation, the time to failure occurs.","metadata":{}},{"cell_type":"markdown","source":"### 2.3 Data Preparation\n\n<font color=#6698FF> Standard normalize the data using Scikit preprocessing.","metadata":{"id":"4A5GsVzgATRu"}},{"cell_type":"code","source":"#standard normalize the data\n\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(X_train)\nX_train = scaler.transform(X_train)\nX_test = scaler.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:22.262012Z","iopub.execute_input":"2024-04-21T07:01:22.263049Z","iopub.status.idle":"2024-04-21T07:01:22.280407Z","shell.execute_reply.started":"2024-04-21T07:01:22.262996Z","shell.execute_reply":"2024-04-21T07:01:22.278839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 3. Training and Results\n\n<font color=#6698FF> For this project both linear and non-linear model has to be used. For the linear model, LinearRegression is used because it is simple. Since we are working with a large data set, it is handy to use a simple model for the linear model part. A drawback of using LinearRegression is that it is sensitive for outliers and noise. \n\n<font color=#6698FF> For the non-linear model several options can be explored. For this project, the following non-linear models are explored: NuSVR and Random Forest Regression. An advantage of NuSVR is it provides a better control over model complexity. A disadvantage is that it can still be sensitive to outliers. An advantage of RandomForestRegression is that it is less prone to overfitting. A disadvantage is it offers limited insight of the relation between the features and target variables.\n    \n<font color=#6698FF> Each models, both linear and non-linear, are tested on the test and training data. The test and train score are calculated and printed. This is also done with the MSE. Based on the performances of the different models, the one with the lowest MSE score will be used since this indicates that the data points are distributed closely to the mean.","metadata":{"id":"IUcc2ue7yrOv"}},{"cell_type":"code","source":"#Linear model: LinearRegression\n\nfrom sklearn.linear_model import LinearRegression\n#define the model \nlinear = LinearRegression()\n#fit the model\nlinear.fit(X_train, y_train)\n\n#checking the score for both the train and test data\n\nprint('Train score:', linear.score(X_train, y_train))\nprint('Test score:', linear.score(X_test, y_test))\n\nmse= mean_squared_error(y_test, linear.predict(X_test))\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:25.364192Z","iopub.execute_input":"2024-04-21T07:01:25.365471Z","iopub.status.idle":"2024-04-21T07:01:25.398398Z","shell.execute_reply.started":"2024-04-21T07:01:25.365380Z","shell.execute_reply":"2024-04-21T07:01:25.391954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Non-linear model: NuSVR\nfrom sklearn.svm import NuSVR, SVR\n\nsvm = NuSVR()\nsvm.fit(X_train, y_train)\n\n#check the model score:\nscore = svm.score(X_test, y_test)\nprint('Test score:', score)\nprint('Train score:',svm.score(X_train, y_train))\n\nmse = mean_squared_error(y_test, svm.predict(X_test))\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:28.090935Z","iopub.execute_input":"2024-04-21T07:01:28.091927Z","iopub.status.idle":"2024-04-21T07:01:30.727043Z","shell.execute_reply.started":"2024-04-21T07:01:28.091880Z","shell.execute_reply":"2024-04-21T07:01:30.725482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Non-linear model: Random Forest Regression\n\nfrom sklearn.ensemble import RandomForestRegressor\n\n#fitting the random forest regression to the dataset\n#regressor = RandomForestRegressor(n_estimators = 200, random_state = 40, oob_state = True)\nregressor = RandomForestRegressor(n_estimators = 200, random_state = 40)\n#fitting the regressor to the x and y data\nregressor.fit(X_train, y_train)\n\n#printing the oob score\n#print('out-of-bag score:', regressor.oob_score_)\n\n#checking the score of the model to see whether it is a good model for both the testing and training data\nprint('Test score:',regressor.score(X_test, y_test))\nprint('Train score:', regressor.score(X_train, y_train))\n\nmse = mean_squared_error(y_test, regressor.predict(X_test))\nprint(mse)\n","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:33.648331Z","iopub.execute_input":"2024-04-21T07:01:33.649262Z","iopub.status.idle":"2024-04-21T07:01:44.904681Z","shell.execute_reply.started":"2024-04-21T07:01:33.649215Z","shell.execute_reply":"2024-04-21T07:01:44.903333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Plot the predicted time_to_failure of the entire dataset with the actual predicted time_to_failure\nfig, ax1 = plt.subplots(figsize=(20, 4))\nplt.title(\"NuSVR model with the data\")\nplt.plot(svm.predict(scaler.transform(X)), color='b')\nax1.set_ylabel('predicted time_to_failure')\nplt.legend(['predicted time_to_failure'], loc=(0.7, 0.9))\nax2 = ax1.twinx()\nplt.plot(y, color='r')\nax2.set_ylabel('time_to_failure')\nplt.legend(['time_to_failure'], loc=(0.8, 0.9))","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:01:48.752319Z","iopub.execute_input":"2024-04-21T07:01:48.752798Z","iopub.status.idle":"2024-04-21T07:01:49.798955Z","shell.execute_reply.started":"2024-04-21T07:01:48.752756Z","shell.execute_reply":"2024-04-21T07:01:49.797556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion\n\n<font color=#6698FF> After trying out different models and comparing them with each other, the result as shown above are the best that could be achieved with the tools and knowledge available. Out of the three models that were used, the NuSVR performs the best. This can be seen by the low MSE score and the score on the test data. On the train data, the NuSVR falls short compared to the RandomForestRegression but in the end the NuSVR is chosen. \n    \n<font color=#6698FF> The biggest difference between the score on the train and test data can be observed when using the RandomForestRegression. This might be because the model is overfitted. To further improve the model, exploring other features might help. Finetuning the hyperparameters to avoid overfitting might also help improve the performances of the models.","metadata":{"id":"yekgnxu7y0gE"}},{"cell_type":"code","source":"#Saving the submission document into the desirable format\n\nsubmission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(dtype=np.float64, index=submission.index)\n\nfor seg_id in tqdm(X_test.index):\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv')\n    \n    x = pd.Series(seg['acoustic_data'].values)\n    X_test = x_features(x, X_test, seg_id)\n    \n#checking the shape of the to be submitted data\n\nprint(x.shape)\nprint(X_test.shape)\n\n#scale the test data\n\nX_test = scaler.transform(X_test)\nsubmission['time_to_failure'] = svm.predict(X_test)\nsubmission.to_csv('submission.csv')\nsubmission.head(10)\n\nprint(submission)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T07:15:38.905460Z","iopub.execute_input":"2024-04-21T07:15:38.906309Z","iopub.status.idle":"2024-04-21T07:21:20.257943Z","shell.execute_reply.started":"2024-04-21T07:15:38.906264Z","shell.execute_reply":"2024-04-21T07:21:20.256714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 5. References","metadata":{"id":"wn7PBEmWESIi"}},{"cell_type":"markdown","source":"* LANL Earthquake EDA and Prediction (gpreda) https://www.kaggle.com/code/gpreda/lanl-earthquake-eda-and-prediction#Prepare-the-data-analysis\n* ML23_week8 (jessteijn) https://www.kaggle.com/code/jessteijn/ml23-week8\n* Noura_Hannah_earthquake (nouradeklerk) https://www.kaggle.com/code/nouradeklerk/noura-hannah-earthquake\n* Earthquakes FE. More features and samples (artgor) https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples\n* Discussion: https://www.kaggle.com/competitions/LANL-Earthquake-Prediction/discussion/84601","metadata":{}}]}