{"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":"# 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> Please write your name, username in Kaggle, score and Leaderboard rank here.</font>\n\n    Names: Coen Buitink\n           Jeroen Lohman\n    \n    Usernames Kaggle: coenbuitink\n                      JeroenLoL\n    Leaderboard rank:\n","metadata":{"id":"DPJV5e8sDVnT"}},{"cell_type":"markdown","source":"## 2. Data\n","metadata":{"id":"PB7FLdQn-dmK"}},{"cell_type":"markdown","source":"### 2.1 Dataset\n<font color=#6698FF>In this section, we load and explore the dataset.","metadata":{"id":"vfGikOxUAxJB"}},{"cell_type":"markdown","source":"**Libraries**","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport tensorflow as tf\nimport tensorflow_decision_forests as tfdf\nfrom tqdm import tqdm_notebook\nfrom scipy import stats","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:00:25.366827Z","iopub.execute_input":"2023-04-28T19:00:25.367219Z","iopub.status.idle":"2023-04-28T19:00:35.605169Z","shell.execute_reply.started":"2023-04-28T19:00:25.367185Z","shell.execute_reply":"2023-04-28T19:00:35.603736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"TensorFlow v\" + tf.__version__)\nprint(\"TensorFlow Decision Forests v\" + tfdf.__version__)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:00:35.607714Z","iopub.execute_input":"2023-04-28T19:00:35.608372Z","iopub.status.idle":"2023-04-28T19:00:35.614671Z","shell.execute_reply.started":"2023-04-28T19:00:35.608331Z","shell.execute_reply":"2023-04-28T19:00:35.613559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load the Data**","metadata":{}},{"cell_type":"code","source":"path_training_data = \"/kaggle/input/LANL-Earthquake-Prediction/train.csv\"\npath_test_data = \"/kaggle/input/LANL-Earthquake-Prediction/test\"\ntraining_data = pd.read_csv(path_training_data, dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\n\nprint(training_data.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:00:35.616038Z","iopub.execute_input":"2023-04-28T19:00:35.616584Z","iopub.status.idle":"2023-04-28T19:04:42.922121Z","shell.execute_reply.started":"2023-04-28T19:00:35.616551Z","shell.execute_reply":"2023-04-28T19:04:42.920935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Data Exploration\n\n<font color=#6698FF> Explore the features and target variables of the dataset. Think about making some scatter plots, box plots, histograms or printing the data, but feel free to choose any method that suits you.\nWhat do you think is the right performance\nmetric to use for this dataset? Clearly explain which performance metric you\nchoose and why.\nAlgorithmic bias can be a real problem in Machine Learning.  Explain what you believe.</font>","metadata":{"id":"QN5jS2_p_Y9N"}},{"cell_type":"code","source":"training_data.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:04:42.923830Z","iopub.execute_input":"2023-04-28T19:04:42.924561Z","iopub.status.idle":"2023-04-28T19:04:42.957528Z","shell.execute_reply.started":"2023-04-28T19:04:42.924519Z","shell.execute_reply":"2023-04-28T19:04:42.956058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot the raw data**\nIn this section, plotting the data coded using this source, only now plotting the entire dataset to see the periodicity:\nhttps://www.kaggle.com/code/ahemateja19bec1025/data-visualization-earthquake-review3-19bec1025\n","metadata":{}},{"cell_type":"code","source":"#visualize 1% of samples data, first 100 datapoints\ntrain_acoustic_data_sample_df = training_data['acoustic_data'].values[::100]\ntrain_ttf_sample_df = training_data['time_to_failure'].values[::100]\n\n#function for plotting based on both features\ndef plot_acc_ttf_data(train_acoustic_data_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure over all the data\"):\n    fig, ax1 = plt.subplots(figsize=(12, 8))\n    plt.title(title)\n    plt.plot(train_acoustic_data_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_acoustic_data_sample_df, train_ttf_sample_df)\ndel train_acoustic_data_sample_df\ndel train_ttf_sample_df","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:04:42.961153Z","iopub.execute_input":"2023-04-28T19:04:42.961593Z","iopub.status.idle":"2023-04-28T19:04:46.600513Z","shell.execute_reply.started":"2023-04-28T19:04:42.961552Z","shell.execute_reply":"2023-04-28T19:04:46.598922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Compare this with test data**","metadata":{}},{"cell_type":"code","source":"# Load in the first three data files\ntest_1 = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/seg_00030f.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\ntest_2 = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/seg_0012b5.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\ntest_3 = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/seg_00184e.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\n\n# Initialize the subplots\nfig, ax = plt.subplots(3, 1, figsize=(16,6))\n\n# Plot the test files\nax[0].plot(test_1, color='g')\nax[0].legend(['acoustic_data'])\nax[0].set_ylabel('acoustic_data')\n\nax[1].plot(test_2, color='g')\nax[1].legend(['acoustic_data'])\nax[1].set_ylabel('acoustic_data')\n\nax[2].plot(test_3, color='g')\nax[2].legend(['acoustic_data'])\nax[2].set_ylabel('acoustic_data')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:04:46.603023Z","iopub.execute_input":"2023-04-28T19:04:46.604180Z","iopub.status.idle":"2023-04-28T19:04:48.248024Z","shell.execute_reply.started":"2023-04-28T19:04:46.604117Z","shell.execute_reply":"2023-04-28T19:04:48.246469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(test_1), len(test_2), len(test_3))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:04:48.249844Z","iopub.execute_input":"2023-04-28T19:04:48.250299Z","iopub.status.idle":"2023-04-28T19:04:48.256797Z","shell.execute_reply.started":"2023-04-28T19:04:48.250259Z","shell.execute_reply":"2023-04-28T19:04:48.255500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since the test data has files of 150000 samples long, the training data should be analysed in similar sized 'rows'.","metadata":{}},{"cell_type":"markdown","source":"### 2.3 Data Preparation\n\n<font color=#6698FF>  This dataset hasn’t been cleaned yet. Meaning that some attributes (features) are in numerical format and some are in categorial format. Moreover, there are missing values as well. However, all Scikit-learn’s implementations of these algorithms expect numerical features. Check for all features if they are in categorial and use a method to transform them to numerical values. For the numerical data, handle the missing data and normalize the data. \nNote that you are only allowed to use training data for preprocessing but you then need to perform similar changes on test data too.\nYou can use [pipelining](https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html) to help with the preprocessing.</font>","metadata":{"id":"4A5GsVzgATRu"}},{"cell_type":"markdown","source":"**Feature Extraction**\nhttps://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples","metadata":{}},{"cell_type":"code","source":"from tqdm import tqdm\nfrom scipy import stats\n\n#Make 150000 rows, since the submission sample data has files of 150000 datapoints\nrows = 150000\nsegments = int(np.floor(training_data.shape[0] / rows))\nX = pd.DataFrame(index=range(segments), dtype=np.float64)\nY = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\n\nfor segment in tqdm(range(segments)):\n    seg = training_data.iloc[segment*rows:segment*rows+rows]\n    x = pd.Series(seg['acoustic_data'].values)\n    y = seg['time_to_failure'].values[-1]\n    \n    Y.loc[segment, 'time_to_failure'] = y\n    \n    # Basic statistical measures \n    X.loc[segment, 'ave'] = x.mean()\n    X.loc[segment, 'std'] = x.std()\n    X.loc[segment, 'max'] = x.max()\n    X.loc[segment, 'min'] = x.min()\n    \n    X.loc[segment, 'mean_change_abs'] = np.mean(np.diff(x))\n    X.loc[segment, 'abs_max'] = np.abs(x).max()\n    X.loc[segment, 'abs_min'] = np.abs(x).min()\n    \n    # Basic statistical measures with intervals\n    X.loc[segment, 'std_first_50000'] = x[:50000].std()\n    X.loc[segment, 'std_last_50000'] = x[-50000:].std()\n    X.loc[segment, 'std_first_10000'] = x[:10000].std()\n    X.loc[segment, 'std_last_10000'] = x[-10000:].std()\n    \n    X.loc[segment, 'avg_first_50000'] = x[:50000].mean()\n    X.loc[segment, 'avg_last_50000'] = x[-50000:].mean()\n    X.loc[segment, 'avg_first_10000'] = x[:10000].mean()\n    X.loc[segment, 'avg_last_10000'] = x[-10000:].mean()\n    \n    X.loc[segment, 'min_first_50000'] = x[:50000].min()\n    X.loc[segment, 'min_last_50000'] = x[-50000:].min()\n    X.loc[segment, 'min_first_10000'] = x[:10000].min()\n    X.loc[segment, 'min_last_10000'] = x[-10000:].min()\n    \n    X.loc[segment, 'max_first_50000'] = x[:50000].max()\n    X.loc[segment, 'max_last_50000'] = x[-50000:].max()\n    X.loc[segment, 'max_first_10000'] = x[:10000].max()\n    X.loc[segment, 'max_last_10000'] = x[-10000:].max()\n    \n    X.loc[segment, 'max_to_min'] = x.max() / np.abs(x.min())\n    X.loc[segment, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n    X.loc[segment, 'count_big'] = len(x[np.abs(x) > 500])\n    X.loc[segment, 'sum'] = x.sum()\n    \n    # Quantiles\n    X.loc[segment, 'q95'] = np.quantile(x, 0.95)\n    X.loc[segment, 'q99'] = np.quantile(x, 0.99)\n    X.loc[segment, 'q05'] = np.quantile(x, 0.05)\n    X.loc[segment, 'q01'] = np.quantile(x, 0.01)\n    \n    X.loc[segment, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n    X.loc[segment, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n    X.loc[segment, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n    X.loc[segment, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    # Rolling statistical measures\n    for windows in [10, 100, 1000]:\n        x_roll_std = x.rolling(windows).std().dropna().values\n        x_roll_mean = x.rolling(windows).mean().dropna().values\n        \n        X.loc[segment, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X.loc[segment, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X.loc[segment, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X.loc[segment, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X.loc[segment, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X.loc[segment, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X.loc[segment, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X.loc[segment, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X.loc[segment, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        X.loc[segment, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X.loc[segment, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X.loc[segment, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X.loc[segment, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X.loc[segment, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X.loc[segment, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X.loc[segment, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X.loc[segment, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X.loc[segment, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X.loc[segment, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        X.loc[segment, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:04:48.258117Z","iopub.execute_input":"2023-04-28T19:04:48.258490Z","iopub.status.idle":"2023-04-28T19:13:04.917536Z","shell.execute_reply.started":"2023-04-28T19:04:48.258455Z","shell.execute_reply":"2023-04-28T19:13:04.915608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Check the shapes of the created data\nprint(X.shape)\nprint(Y.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:04.920047Z","iopub.execute_input":"2023-04-28T19:13:04.921171Z","iopub.status.idle":"2023-04-28T19:13:04.930394Z","shell.execute_reply.started":"2023-04-28T19:13:04.921099Z","shell.execute_reply":"2023-04-28T19:13:04.929094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:04.932264Z","iopub.execute_input":"2023-04-28T19:13:04.932688Z","iopub.status.idle":"2023-04-28T19:13:04.982016Z","shell.execute_reply.started":"2023-04-28T19:13:04.932655Z","shell.execute_reply":"2023-04-28T19:13:04.980812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The train test split is performed here, because now the features could be extracted first and then be splitted. The majority (70% or 80%) should be allocated for training data, since the given data doesn't have too many distinctive features. Therefore, it is better to use enough data for training to make a model.","metadata":{}},{"cell_type":"code","source":"# Train test split\nfrom sklearn.model_selection import train_test_split\nX_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.3,shuffle=False, random_state=102)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:04.983927Z","iopub.execute_input":"2023-04-28T19:13:04.984266Z","iopub.status.idle":"2023-04-28T19:13:05.200295Z","shell.execute_reply.started":"2023-04-28T19:13:04.984234Z","shell.execute_reply":"2023-04-28T19:13:05.199229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(Y_train.shape)\nprint(X_train.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:05.201564Z","iopub.execute_input":"2023-04-28T19:13:05.201946Z","iopub.status.idle":"2023-04-28T19:13:05.209358Z","shell.execute_reply.started":"2023-04-28T19:13:05.201912Z","shell.execute_reply":"2023-04-28T19:13:05.208364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Standard normalize the data\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":"2023-04-28T19:13:05.211128Z","iopub.execute_input":"2023-04-28T19:13:05.211486Z","iopub.status.idle":"2023-04-28T19:13:05.231361Z","shell.execute_reply.started":"2023-04-28T19:13:05.211438Z","shell.execute_reply":"2023-04-28T19:13:05.230146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 3. Training and Results\n<font color=#6698FF> Briefly introduce algorithms you choose. \nPresent your balanced accuracies for both test and training data for all classifiers. Analyse the performance on test and training in terms of bias and variance. Give one advantage and one drawback of the method you use.\n","metadata":{"id":"IUcc2ue7yrOv"}},{"cell_type":"markdown","source":"**Linear Model**\n<br>\nFirst make plain linear regression model","metadata":{}},{"cell_type":"code","source":"from sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import mean_absolute_error, mean_squared_error\nmodel_LR= LinearRegression()\n\n# fit the model to data\nmodel_LR.fit(X_train,Y_train)\n\npredictions = model_LR.predict(X_test)\n#print(predictions)\nscore = model_LR.score(X_test, Y_test)\n#print(\"\")\nprint(\"Accuracy: \" + str(score))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:05.236167Z","iopub.execute_input":"2023-04-28T19:13:05.236517Z","iopub.status.idle":"2023-04-28T19:13:05.374741Z","shell.execute_reply.started":"2023-04-28T19:13:05.236482Z","shell.execute_reply":"2023-04-28T19:13:05.370668Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"MSE: \", mean_squared_error(Y_test, predictions))\nprint(\"MAE: \", mean_absolute_error(Y_test, predictions))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:05.376546Z","iopub.execute_input":"2023-04-28T19:13:05.377007Z","iopub.status.idle":"2023-04-28T19:13:05.402553Z","shell.execute_reply.started":"2023-04-28T19:13:05.376961Z","shell.execute_reply":"2023-04-28T19:13:05.400604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax1 = plt.subplots(figsize=(20, 4))\nt = np.linspace(0, 1259, 1259)\nplt.title(\"Linear Regression with test splitted data prediction\")\nplt.plot(predictions, color='g')\nax1.set_ylabel('Time to failure', color='g')\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(['Real time to failure'], loc=(0.07, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:05.409847Z","iopub.execute_input":"2023-04-28T19:13:05.411134Z","iopub.status.idle":"2023-04-28T19:13:05.866006Z","shell.execute_reply.started":"2023-04-28T19:13:05.411069Z","shell.execute_reply":"2023-04-28T19:13:05.864554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**PCA with linear regression to reduce feature dimensions**","metadata":{}},{"cell_type":"code","source":"# Train pca\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.decomposition import PCA\n\n\n# Selecting 24 features from the total number of features (right amount to get a >90% cumulative variance)\nn_components=24\npca = PCA(n_components) \npca = pca.fit(X_train, Y_train)\n\n# Perform pca on features\nX_train_pca=pca.transform(X_train);\nX_test_pca=pca.transform(X_test)\n\n# Training and testing on model\n\nmodel = LinearRegression()\n# Fit the model\nmodel.fit(X_train_pca, Y_train)\n\npredictions_pca_train = model.predict(X_train_pca)\npredictions_pca_test = model.predict(X_test_pca)\n\nscore_train = model.score(X_train_pca, Y_train)\n\nprint(\"Training results:\")\nprint(\"Accuracy score with PCA: \" + str(score_train))\n\nprint(\"MSE with PCA: \", mean_squared_error(Y_train, predictions_pca_train))\nprint(\"MAE with PCA: \",mean_absolute_error(Y_train, predictions_pca_train))\n\nprint(\"\")\n\nscore_test = model.score(X_test_pca, Y_test)\n\nprint(\"Test results:\")\nprint(\"Accuracy score with PCA: \" + str(score_test))\n\nprint(\"MSE with PCA: \", mean_squared_error(Y_test, predictions_pca_test))\nprint(\"MAE with PCA: \",mean_absolute_error(Y_test, predictions_pca_test))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:05.867873Z","iopub.execute_input":"2023-04-28T19:13:05.868287Z","iopub.status.idle":"2023-04-28T19:13:06.017917Z","shell.execute_reply.started":"2023-04-28T19:13:05.868253Z","shell.execute_reply":"2023-04-28T19:13:06.013885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Prediction of train (linear)\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Training Data\")\nplt.plot(predictions_pca_train, color='g')\nax1.set_ylabel('Time to failure', color='g')\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(['Real 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(\"Test Data\")\nplt.plot(predictions_pca_test, color='g')\nax1.set_ylabel('Time to failure', color='g')\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(['Real time to failure'], loc=(0.07, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:06.020007Z","iopub.execute_input":"2023-04-28T19:13:06.020435Z","iopub.status.idle":"2023-04-28T19:13:06.869653Z","shell.execute_reply.started":"2023-04-28T19:13:06.020390Z","shell.execute_reply":"2023-04-28T19:13:06.868797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.bar(range(0,n_components), pca.explained_variance_ratio_, label=\"individual var\");\nplt.step(range(0,n_components), np.cumsum(pca.explained_variance_ratio_),'r', label=\"cumulative var\");\nplt.xlabel('Principal component index'); plt.ylabel('explained variance ratio %');\nplt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:06.871094Z","iopub.execute_input":"2023-04-28T19:13:06.872041Z","iopub.status.idle":"2023-04-28T19:13:07.190942Z","shell.execute_reply.started":"2023-04-28T19:13:06.872005Z","shell.execute_reply":"2023-04-28T19:13:07.189664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"With the amount of statistical features selected, it was not expected of this linear model to perform rather impressively. Moreover, it is most probable that the features of the data have a non-linear correlation, meaning that fitting a linear model with great accuracy is hard. Nontheless, the model has been optimized with PCA to improve the accuracy and MSE.","metadata":{}},{"cell_type":"markdown","source":"**Non-Linear Model**","metadata":{}},{"cell_type":"code","source":"#Random forest\nfrom sklearn.ensemble import RandomForestRegressor\nrf = RandomForestRegressor(max_depth = 4, random_state = 42)\nrf.fit(X_train, np.ravel(Y_train))\npredictions_rf = rf.predict(X_test)\nscore_rf = rf.score(X_test, Y_test)\n\nprint(\"Accuracy score Random Forest: \" + str(score_rf))\nMSE_rf = np.sum((predictions_rf.reshape(len(Y_test),1) - Y_test)**2)/Y_test.size\nprint(\"MSE with RF: \" + str(float(MSE_rf)))\nprint(\"MAE with RF: \",mean_absolute_error(Y_test, predictions_rf))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:07.193825Z","iopub.execute_input":"2023-04-28T19:13:07.196071Z","iopub.status.idle":"2023-04-28T19:13:12.854084Z","shell.execute_reply.started":"2023-04-28T19:13:07.196016Z","shell.execute_reply":"2023-04-28T19:13:12.852824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax1 = plt.subplots(figsize=(20, 4))\nt = np.linspace(0, 1259, 1259)\nplt.title(\"Linear Regression with test splitted data prediction\")\nplt.plot(predictions_rf, color='g')\nax1.set_ylabel('Time to failure', color='g')\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(['Real time to failure'], loc=(0.07, 0.85))\nplt.grid(True)\nplt.xlim(0, 1100)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:12.857052Z","iopub.execute_input":"2023-04-28T19:13:12.857840Z","iopub.status.idle":"2023-04-28T19:13:13.295328Z","shell.execute_reply.started":"2023-04-28T19:13:12.857791Z","shell.execute_reply":"2023-04-28T19:13:13.293682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Submission**","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id') # submission csv\nx_test = pd.DataFrame(dtype=np.float64, index=submission.index)\n    \nfor seg_id in x_test.index:\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv') # read every test data and extract features\n    \n    x = pd.Series(seg['acoustic_data'].values)\n    \n    x_test.loc[seg_id, 'ave'] = x.mean()\n    x_test.loc[seg_id, 'std'] = x.std()\n    x_test.loc[seg_id, 'max'] = x.max()\n    x_test.loc[seg_id, 'min'] = x.min()\n    \n    x_test.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(x))\n    x_test.loc[seg_id, 'abs_max'] = np.abs(x).max()\n    x_test.loc[seg_id, 'abs_min'] = np.abs(x).min()\n    \n    x_test.loc[seg_id, 'std_first_50000'] = x[:50000].std()\n    x_test.loc[seg_id, 'std_last_50000'] = x[-50000:].std()\n    x_test.loc[seg_id, 'std_first_10000'] = x[:10000].std()\n    x_test.loc[seg_id, 'std_last_10000'] = x[-10000:].std()\n    \n    x_test.loc[seg_id, 'avg_first_50000'] = x[:50000].mean()\n    x_test.loc[seg_id, 'avg_last_50000'] = x[-50000:].mean()\n    x_test.loc[seg_id, 'avg_first_10000'] = x[:10000].mean()\n    x_test.loc[seg_id, 'avg_last_10000'] = x[-10000:].mean()\n\n    x_test.loc[seg_id, 'min_first_50000'] = x[:50000].min()\n    x_test.loc[seg_id, 'min_last_50000'] = x[-50000:].min()\n    x_test.loc[seg_id, 'min_first_10000'] = x[:10000].min()\n    x_test.loc[seg_id, 'min_last_10000'] = x[-10000:].min()\n    \n    x_test.loc[seg_id, 'max_first_50000'] = x[:50000].max()\n    x_test.loc[seg_id, 'max_last_50000'] = x[-50000:].max()\n    x_test.loc[seg_id, 'max_first_10000'] = x[:10000].max()\n    x_test.loc[seg_id, 'max_last_10000'] = x[-10000:].max()\n    \n    x_test.loc[seg_id, 'max_to_min'] = x.max() / np.abs(x.min())\n    x_test.loc[seg_id, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n    x_test.loc[seg_id, 'count_big'] = len(x[np.abs(x) > 500])\n    x_test.loc[seg_id, 'sum'] = x.sum()\n\n    x_test.loc[seg_id, 'q95'] = np.quantile(x, 0.95)\n    x_test.loc[seg_id, 'q99'] = np.quantile(x, 0.99)\n    x_test.loc[seg_id, 'q05'] = np.quantile(x, 0.05)\n    x_test.loc[seg_id, 'q01'] = np.quantile(x, 0.01)\n    \n    x_test.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n    x_test.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n    x_test.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n    x_test.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    for windows in [10, 100, 1000]:\n        x_roll_std = x.rolling(windows).std().dropna().values\n        x_roll_mean = x.rolling(windows).mean().dropna().values\n        \n        x_test.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        x_test.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        x_test.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        x_test.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        x_test.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        x_test.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        x_test.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        x_test.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        x_test.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        x_test.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        x_test.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        x_test.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        x_test.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        x_test.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        x_test.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        x_test.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        x_test.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        x_test.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        x_test.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        x_test.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:13:13.297480Z","iopub.execute_input":"2023-04-28T19:13:13.298256Z","iopub.status.idle":"2023-04-28T19:19:08.146119Z","shell.execute_reply.started":"2023-04-28T19:13:13.298209Z","shell.execute_reply":"2023-04-28T19:19:08.145051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_test.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:19:08.147630Z","iopub.execute_input":"2023-04-28T19:19:08.148100Z","iopub.status.idle":"2023-04-28T19:19:08.185616Z","shell.execute_reply.started":"2023-04-28T19:19:08.148067Z","shell.execute_reply":"2023-04-28T19:19:08.184413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaler = preprocessing.StandardScaler().fit(x_test)\nx_test = scaler.transform(x_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:19:08.186988Z","iopub.execute_input":"2023-04-28T19:19:08.187419Z","iopub.status.idle":"2023-04-28T19:19:08.204650Z","shell.execute_reply.started":"2023-04-28T19:19:08.187369Z","shell.execute_reply":"2023-04-28T19:19:08.203509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission['time_to_failure'] = rf.predict(x_test) \nsubmission.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:19:08.206011Z","iopub.execute_input":"2023-04-28T19:19:08.206374Z","iopub.status.idle":"2023-04-28T19:19:08.247878Z","shell.execute_reply.started":"2023-04-28T19:19:08.206331Z","shell.execute_reply":"2023-04-28T19:19:08.246847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:19:08.249216Z","iopub.execute_input":"2023-04-28T19:19:08.249590Z","iopub.status.idle":"2023-04-28T19:19:08.260430Z","shell.execute_reply.started":"2023-04-28T19:19:08.249556Z","shell.execute_reply":"2023-04-28T19:19:08.259164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion\n\n<font color=#6698FF>Discuss all the choices you made during the process and your final conclusions.  Highlight the strong points of your approach, discuss its shortcomings  and suggest some future approaches that may improve it. Please be self critical here. The assignment is not about achieving a state of the art performance, but about showing what you have learned the concepts during the course.</font>","metadata":{"id":"yekgnxu7y0gE"}},{"cell_type":"markdown","source":"In this assignment, all the obtained knowlegde about feature extraction methods and regression models was needed to come to a functional prediction model. The used models were: linear regression, PCA + linear regression and as for non linear the random forest regressor. To start with the data exploration, the entire training dataset was plotted to get a view of what data was given for this assignment. Here it could be seen how before every rise in the 'time to failure' was preceded by a spike in the 'acoustic data'. From these repeated spikes was derived that the given data consisted of segments. Furthermore, a few submission sample were plotted to see how the final model should be used to analyse the performance. From these plots was derived as was stated before that each test file was the same amount of samples. Next, with the help of other notebooks like from user 'ANDREW LUKYANENKO' the feature detection algorithm was implemented. \n<br>\n<br>\nAfter the feature extraction, a train test split was performed. Here, a critical mistake was made here that congested progression of this assignment for a long time. Namely that you shouldn't randomly shuffle this train test split, since else it is impossible to plot your prediction performance later on. Nevertheless, the first model implemented was a plain linear regression model with all the 95 features using SKLearn. Then the amount of features was reduced for the linear regression model using PCA and keeping 24 features. This reduction in features was carefully chosen to have a sufficient cumulative variance, as well as a good accuracy for the training and test data (from the split). \nFinally, for the non linear model a random forest regressor was implemented. Here the consideration was made between a good accuracy against a overfitted model, which is very likely to occur with a random forest model. A larger 'max_depth' means in general a better accuracy, but also overfits the model to the training data. \n<br>\n<br>\nFor the measurement of the models, the MSE as well as the MAE were used to evaluate the performance. Overall, the random forest regressor had the lowest MSE and very similar MAE and accuracy score compared to the other tested models. Therefore the random forest model was used for predicting the time to failure of the test data. ","metadata":{}},{"cell_type":"markdown","source":"## 5. References","metadata":{"id":"wn7PBEmWESIi"}},{"cell_type":"markdown","source":"Data visualisation: https://www.kaggle.com/code/ahemateja19bec1025/data-visualization-earthquake-review3-19bec1025\nFeature extraction and submission code: https://www.kaggle.com/code/wimwim/rolling-quantiles/notebook and https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples\nGeneral perception of the assignmnt: https://www.kaggle.com/code/allunia/shaking-earth","metadata":{}}]}