{"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","metadata":{"id":"3NnAgB-y0knj"}},{"cell_type":"markdown","source":"<font color=#6698FF> Name : Robbert Bithray and Bert van lange \n    \n<font color=#6698FF>Username: robbertbithray\n\n<font color=#6698FF>score :\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":"code","source":"import numpy as np # linear algebra\nimport pandas as pd \nimport matplotlib.pyplot as plt # For plots\n\n\nimport os\nprint(os.listdir(\"../input/LANL-Earthquake-Prediction\"))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:45:05.534688Z","iopub.execute_input":"2023-04-28T21:45:05.535836Z","iopub.status.idle":"2023-04-28T21:45:05.544764Z","shell.execute_reply.started":"2023-04-28T21:45:05.535782Z","shell.execute_reply":"2023-04-28T21:45:05.543374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64}) # read the csv data file\nprint(train.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T19:58:46.262267Z","iopub.execute_input":"2023-04-28T19:58:46.262717Z","iopub.status.idle":"2023-04-28T20:01:16.685144Z","shell.execute_reply.started":"2023-04-28T19:58:46.262674Z","shell.execute_reply":"2023-04-28T20:01:16.683802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This was used during development with a pre generated lest of feature values based on `extract_features()` and the provided data set\n```python\nfeatures = np.loadtxt('../input/lanl-earthquake-features/features.csv', delimiter = ',')\ny_train = np.loadtxt('../input/lanl-earthquake-features/Y_train.csv', delimiter = ',')[:,1]\nprint(features.shape)\nprint(y_train.shape)\n```","metadata":{"execution":{"iopub.status.busy":"2023-04-28T20:38:10.475467Z","iopub.execute_input":"2023-04-28T20:38:10.476692Z","iopub.status.idle":"2023-04-28T20:38:10.950635Z","shell.execute_reply.started":"2023-04-28T20:38:10.476642Z","shell.execute_reply":"2023-04-28T20:38:10.949370Z"}}},{"cell_type":"markdown","source":"Plot the training set, to get a fealing of the input.\n\n","metadata":{}},{"cell_type":"code","source":"fig,axis = plt.subplots()\naxis.plot(train['acoustic_data'].values[::100], label =\"axoustic_data\",color=\"blue\")\naxis.set_xlabel(\"samples/100\")\naxis.set_ylabel(\"audio amplitude\")\naxis2= axis.twinx()\naxis2.set_ylabel(\"time to failure\")\naxis2.plot(train['time_to_failure'].values[::100], label = \"time to failure\",color = \"green\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T13:39:55.412541Z","iopub.execute_input":"2023-04-28T13:39:55.412974Z","iopub.status.idle":"2023-04-28T13:39:58.374974Z","shell.execute_reply.started":"2023-04-28T13:39:55.412938Z","shell.execute_reply":"2023-04-28T13:39:58.373681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from statistics import  mode\nmode(np.diff(train['time_to_failure'].values))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T14:13:29.610964Z","iopub.execute_input":"2023-04-28T14:13:29.611368Z","iopub.status.idle":"2023-04-28T14:15:22.679240Z","shell.execute_reply.started":"2023-04-28T14:13:29.611334Z","shell.execute_reply":"2023-04-28T14:15:22.677883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"to get a feeling of wat the ml algaridem neads to work with 3 test files are plotted","metadata":{}},{"cell_type":"code","source":"test_file = [\"seg_00030f\",\"seg_0012b5\",\"seg_00184e\",\"seg_003339\"]\n#len of file 15000\n\n\nfig, axs = plt.subplots(len(test_file),sharex=True)\nfig.suptitle(\"Test data\")\nfor x in range(len(test_file)):\n    audio_data = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/'+test_file[x]+'.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\n    axs[x].plot(audio_data)\n    axs[x].set_ylim(-300,300)\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T14:27:14.082032Z","iopub.execute_input":"2023-04-28T14:27:14.082452Z","iopub.status.idle":"2023-04-28T14:27:14.837682Z","shell.execute_reply.started":"2023-04-28T14:27:14.082400Z","shell.execute_reply":"2023-04-28T14:27:14.836438Z"},"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. Set a value for the `test_size` yourself. Argue why the test value can not be too small or too large. You can also use k-fold cross validation.\n","metadata":{"id":"udqAFg1C3JHe"}},{"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":"markdown","source":"","metadata":{"execution":{"iopub.status.busy":"2023-04-28T14:34:46.890932Z","iopub.execute_input":"2023-04-28T14:34:46.891425Z","iopub.status.idle":"2023-04-28T14:34:46.908865Z","shell.execute_reply.started":"2023-04-28T14:34:46.891361Z","shell.execute_reply":"2023-04-28T14:34:46.906889Z"}}},{"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":"code","source":"# Code based on Basic Feature Benchmark by INVERSION, https://www.kaggle.com/code/inversion/basic-feature-benchmark\n# Code based on Earthquakes FE. More features and samples by ANDREW LUKYANENKO, https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples\n# Code based on Rolling + Quantiles by PATRICK YAM, https://www.kaggle.com/code/wimwim/rolling-quantiles/notebook\nfrom scipy import stats\nfrom tqdm import tqdm\n\ndef feature_extraction(x):\n    features = []\n\n    # Basic statistical measures\n    features.append(x.mean())\n    features.append(x.std())\n    features.append(x.max())\n    features.append(x.min())\n    features.append(np.mean(np.diff(x)))\n    features.append(np.abs(x).max())\n    features.append(np.abs(x).min())\n\n    # Basic statistical measures over intervals\n    features.append(x[:50000].std())\n    features.append(x[-50000:].std())\n    features.append(x[:10000].std())\n    features.append(x[-10000:].std())\n\n    features.append(x[:50000].mean())\n    features.append(x[-50000:].mean())\n    features.append(x[:10000].mean())\n    features.append(x[-10000:].mean())\n\n    features.append(x[:50000].min())\n    features.append(x[-50000:].min())\n    features.append(x[:10000].min())\n    features.append(x[-10000:].min())\n\n    features.append(x[:50000].max())\n    features.append(x[-50000:].max())\n    features.append(x[:10000].max())\n    features.append(x[-10000:].max())\n\n    features.append(x.max() / np.abs(x.min()))\n    features.append(x.max() - np.abs(x.min()))\n    features.append(len(x[np.abs(x) > 500]))\n    features.append(x.sum())\n\n    # Quantiles\n    features.append(np.quantile(x, 0.95))\n    features.append(np.quantile(x, 0.99))\n    features.append(np.quantile(x, 0.05))\n    features.append(np.quantile(x, 0.01))\n\n    features.append(np.quantile(np.abs(x), 0.95))\n    features.append(np.quantile(np.abs(x), 0.99))\n    features.append(np.quantile(np.abs(x), 0.05))\n    features.append(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        features.append(x_roll_std.mean())\n        features.append(x_roll_std.std())\n        features.append(x_roll_std.max())\n        features.append(x_roll_std.min())\n        features.append(np.quantile(x_roll_std, 0.01))\n        features.append(np.quantile(x_roll_std, 0.05))\n        features.append(np.quantile(x_roll_std, 0.95))\n        features.append(np.quantile(x_roll_std, 0.99))\n        features.append(np.mean(np.diff(x_roll_std)))\n        features.append(np.abs(x_roll_std).max())\n\n        features.append(x_roll_mean.mean())\n        features.append(x_roll_mean.std())\n        features.append(x_roll_mean.max())\n        features.append(x_roll_mean.min())\n        features.append(np.quantile(x_roll_mean, 0.01))\n        features.append(np.quantile(x_roll_mean, 0.05))\n        features.append(np.quantile(x_roll_mean, 0.95))\n        features.append(np.quantile(x_roll_mean, 0.99))\n        features.append(np.mean(np.diff(x_roll_mean)))\n        features.append(np.abs(x_roll_mean).max())\n\n    features = np.asarray(features, dtype = np.float64)\n\n    return features.T","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:08:39.325356Z","iopub.execute_input":"2023-04-28T21:08:39.325765Z","iopub.status.idle":"2023-04-28T21:08:39.351653Z","shell.execute_reply.started":"2023-04-28T21:08:39.325731Z","shell.execute_reply":"2023-04-28T21:08:39.350698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_train = pd.DataFrame(index=range(int(len(train['acoustic_data'].values)/150000)), dtype=np.float64, columns=['time_to_failure'])\nfeatures = np.empty([int(len(train['acoustic_data'].values)/150000),95])\nfor x in tqdm(range(int(len(train['acoustic_data'].values)/150000))):\n    acoustic_data =  pd.Series(train['acoustic_data'].values[x*150000:x*150000+150000])\n    y_train[x] = train['time_to_failure'].values[x*150000+150000]\n    features[x] = feature_extraction(acoustic_data)\n    ","metadata":{"execution":{"iopub.status.busy":"2023-04-28T20:13:10.279286Z","iopub.execute_input":"2023-04-28T20:13:10.279715Z","iopub.status.idle":"2023-04-28T20:18:45.218204Z","shell.execute_reply.started":"2023-04-28T20:13:10.279676Z","shell.execute_reply":"2023-04-28T20:18:45.217012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The features computed above","metadata":{}},{"cell_type":"code","source":"# Before normalizing the data and applaying the pca\n# Split the data into train and test data\nfrom sklearn.model_selection import train_test_split\nX_train, X_test, y_train_split, y_test = train_test_split(features, y_train, test_size=0.25,shuffle=False, random_state=162)\n\n# Standard normalize the data \n# Use the trainding data to creat a scaler and scale both sets of data\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(X_train)\nX_train = scaler.transform(X_train)\nX_test = scaler.transform(X_test)\n\n\n# Train the pca with the training data \nfrom sklearn.decomposition import PCA\npca = PCA(n_components=20) # doing pca and keeping only 10 components\npca = pca.fit(X_train)\n\n# Preform the pca on the features  (training and test data separate)\nX_train_pca=pca.transform(X_train)\nX_test_pca=pca.transform(X_test)    \n\n# Plotting the training data to see the seperation of the categories\nprint(\"Shape of the data after the pca: \" + str(np.shape(X_train_pca)))\n\n\n# Plot how important each feature is for the classification\nplt.bar(range(len(pca.explained_variance_ratio_)), pca.explained_variance_ratio_);\nplt.step(range(len(pca.explained_variance_ratio_)), np.cumsum(pca.explained_variance_ratio_),'r');\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T20:44:42.083657Z","iopub.execute_input":"2023-04-28T20:44:42.084100Z","iopub.status.idle":"2023-04-28T20:44:42.417909Z","shell.execute_reply.started":"2023-04-28T20:44:42.084062Z","shell.execute_reply":"2023-04-28T20:44:42.416673Z"},"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<\\font>\n","metadata":{"id":"IUcc2ue7yrOv"}},{"cell_type":"code","source":"#includig libary needed for the models\nfrom sklearn.svm import SVR\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.ensemble import RandomForestRegressor","metadata":{"execution":{"iopub.status.busy":"2023-04-28T20:46:02.662262Z","iopub.execute_input":"2023-04-28T20:46:02.662704Z","iopub.status.idle":"2023-04-28T20:46:02.668927Z","shell.execute_reply.started":"2023-04-28T20:46:02.662666Z","shell.execute_reply":"2023-04-28T20:46:02.667567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### training an testing the linear moddel","metadata":{}},{"cell_type":"code","source":"#creating a linear moddel\nlr_model = SVR()\nlr_model.fit(X_train_pca, y_train_split)\n\n\n\n# Evalueating the moddel\ny_pred_lin = lr_model.predict(X_test_pca)\nmse = mean_squared_error(y_test, y_pred_lin)\nprint(\"Mean Squared Error of linear moddel\", mse)\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T20:46:11.502546Z","iopub.execute_input":"2023-04-28T20:46:11.503624Z","iopub.status.idle":"2023-04-28T20:46:12.297968Z","shell.execute_reply.started":"2023-04-28T20:46:11.503579Z","shell.execute_reply":"2023-04-28T20:46:12.296504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create and train the random forest regression model\nrf_model = RandomForestRegressor(n_estimators=100, random_state=0)\nrf_model.fit(X_train_pca, y_train_split)\n\n# Predict the target values for the test set\ny_pred_non_lin = rf_model.predict(X_test_pca)\n\n# Evaluate the performance of the model\nmse = mean_squared_error(y_test, y_pred_non_lin)\nprint(\"Mean Squared Error:\", mse)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T20:47:54.184778Z","iopub.execute_input":"2023-04-28T20:47:54.185247Z","iopub.status.idle":"2023-04-28T20:47:58.828858Z","shell.execute_reply.started":"2023-04-28T20:47:54.185210Z","shell.execute_reply":"2023-04-28T20:47:58.827538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### plotting preformance","metadata":{}},{"cell_type":"code","source":"fig,axis = plt.subplots()\naxis.plot(y_test, label =\"axoustic_data\",color=\"blue\")\naxis.set_xlabel(\"samples/100\")\naxis.set_ylabel(\"time to failure\")\naxis2= axis.twinx()\naxis2.plot(y_pred_lin, label = \"linear moddel\",color = \"green\")\naxis3= axis.twinx()\naxis3.plot(y_pred_non_lin, label = \"non linear moddel\",color = \"red\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T20:47:58.831128Z","iopub.execute_input":"2023-04-28T20:47:58.832081Z","iopub.status.idle":"2023-04-28T20:47:59.192038Z","shell.execute_reply.started":"2023-04-28T20:47:58.832029Z","shell.execute_reply":"2023-04-28T20:47:59.190851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### calculating the time of failure of the test files","metadata":{}},{"cell_type":"code","source":"#creating the result from linear model\nsample_submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv' ,dtype={'seg_id':str, 'time_to_failure': np.float64}) # read the csv data file\nlin_moddel = sample_submission\nnon_lin_moddel =   sample_submission\n\n\nfor index in tdqm(range(len(sample_submission['seg_id']))):\n    test_file =  sample_submission['seg_id'][index]\n    test_data = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/'+test_file+'.csv' ,dtype={'acoustic_data': np.int16})\n    test_data = feature_extraction(pd.DataFrame(test_data))\n    test_data[np.isnan(test_data)] = 0\n    test_data = scaler.transform(test_data.reshape(1,-1))\n    test_data = pca.transform(test_data)\n    \n    lin_moddel.loc['time_to_failure', index]     = lr_model.predict(test_data)\n    non_lin_moddel.loc['time_to_failure', index] = lr_model.predict(test_data)\n\nlin_moddel.to_csv('lin_model_submission.csv')\nnon_lin_moddel.to_csv('non_lin_moddel_submission.csv')\n\n    \n    \n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:45:14.669472Z","iopub.execute_input":"2023-04-28T21:45:14.669903Z","iopub.status.idle":"2023-04-28T21:45:14.711487Z","shell.execute_reply.started":"2023-04-28T21:45:14.669840Z","shell.execute_reply":"2023-04-28T21:45:14.709764Z"},"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 our submission for the Kaggle LANL earthquake prediction contest, we utilized a wide range of features based on the results found in other submissions. And segmented the input data into segents of a length equal to that of the test files. However, there were several potential improvements that we did not implement. These include incorporating overlapping segments during feature extraction, utilizing wavelet feature extraction, implementing K-fold cross-validation, and exploring alternative machine learning algorithms. These enhancements could have further improved the performance and accuracy of our model in predicting earthquakes. But were ommitted due to the added complexity this would have added to the already lengthy project, given we had a reasonably adequate result.","metadata":{}},{"cell_type":"markdown","source":"## 5. References","metadata":{"id":"wn7PBEmWESIi"}},{"cell_type":"markdown","source":"Code from Basic Feature Benchmark by INVERSION, https://www.kaggle.com/code/inversion/basic-feature-benchmark<br>\nCode from Earthquakes FE. More features and samples by ANDREW LUKYANENKO, https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples<br>\nCode from Rolling + Quantiles by PATRICK YAM, https://www.kaggle.com/code/wimwim/rolling-quantiles/notebook<br>","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:25:19.205199Z","iopub.execute_input":"2023-04-28T21:25:19.205988Z","iopub.status.idle":"2023-04-28T21:25:19.215206Z","shell.execute_reply.started":"2023-04-28T21:25:19.205940Z","shell.execute_reply":"2023-04-28T21:25:19.213745Z"}}}]}