{"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":"code","source":"!pip install autograd","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-28T08:35:26.711418Z","iopub.execute_input":"2023-04-28T08:35:26.712321Z","iopub.status.idle":"2023-04-28T08:35:41.570358Z","shell.execute_reply.started":"2023-04-28T08:35:26.712270Z","shell.execute_reply":"2023-04-28T08:35:41.568688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport autograd.numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom matplotlib import pyplot as plt \n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n for filename in filenames:\n #print(os.path.join(dirname, filename))\n pass\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-28T08:35:41.573137Z","iopub.execute_input":"2023-04-28T08:35:41.573705Z","iopub.status.idle":"2023-04-28T08:35:44.934159Z","shell.execute_reply.started":"2023-04-28T08:35:41.573620Z","shell.execute_reply":"2023-04-28T08:35:44.932927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Introduction\n

NOTE: I talked to Bahareh and got some exceptions for the assignment

\nName: Daniël Paul
\nUsername in kaggle: Daniël Paul
\nScore: 4.04954
\nPublic score: 3.41049
\n","metadata":{}},{"cell_type":"markdown","source":"## 2. Data","metadata":{}},{"cell_type":"markdown","source":"#### 2.1.1 Train-test split\nI set the test_size to 0.2. This way, there is enough test data to test the models on, but still enough training data.","metadata":{}},{"cell_type":"code","source":"#%%time\ntrain = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:35:44.935816Z","iopub.execute_input":"2023-04-28T08:35:44.936489Z","iopub.status.idle":"2023-04-28T08:39:58.149768Z","shell.execute_reply.started":"2023-04-28T08:35:44.936453Z","shell.execute_reply":"2023-04-28T08:39:58.148081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#goal: split training set in snippets of 150k\n\ntrain = train.to_numpy()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:39:58.154046Z","iopub.execute_input":"2023-04-28T08:39:58.154942Z","iopub.status.idle":"2023-04-28T08:40:01.096472Z","shell.execute_reply.started":"2023-04-28T08:39:58.154882Z","shell.execute_reply":"2023-04-28T08:40:01.095116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"l_samples = 150000\nn_samples = int(train.shape[0] / l_samples)\nprint(n_samples)\n\ntrain_3d = train[:n_samples*l_samples,:].reshape(n_samples, l_samples, 2)\nX = train_3d[:,:,0]\nY = train_3d[:,l_samples-1,1]\n\nprint(X.shape) \nprint(Y.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:01.098201Z","iopub.execute_input":"2023-04-28T08:40:01.098784Z","iopub.status.idle":"2023-04-28T08:40:01.107209Z","shell.execute_reply.started":"2023-04-28T08:40:01.098742Z","shell.execute_reply":"2023-04-28T08:40:01.105520Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#split data into training and test\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=42)\n\nprint(X_train.shape)\nprint(X_test.shape)\nprint(Y_train.shape)\nprint(Y_test.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:01.108885Z","iopub.execute_input":"2023-04-28T08:40:01.109322Z","iopub.status.idle":"2023-04-28T08:40:02.984098Z","shell.execute_reply.started":"2023-04-28T08:40:01.109285Z","shell.execute_reply":"2023-04-28T08:40:02.982732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Data Exploration\nTo get a sense of the data, X_train and Y_train are visualised.
\nSome right performance metrics include MSE, MAE and R squared. These metrics are suited for regression, unlike accuracy etc. MSE is more vulnerable to outliers than MAE, and later in this notebook they will be tested. R squared is a metric which measures the correlation of Y_true and Y_pred. With perfect regression, it will output '1'. This is a good way to measure the quality of the model, because the goal is to get the best regression. R squared will therefore be used to judge the quality.","metadata":{}},{"cell_type":"code","source":"print(X_train[0:5].shape)\nplt.plot(X_train[0:5].T)\nplt.title('Samples of X_train')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:02.985783Z","iopub.execute_input":"2023-04-28T08:40:02.986145Z","iopub.status.idle":"2023-04-28T08:40:03.796580Z","shell.execute_reply.started":"2023-04-28T08:40:02.986112Z","shell.execute_reply":"2023-04-28T08:40:03.795003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(Y_train, bins=200)\nplt.xlabel('Y_train')\nplt.ylabel('Frequency')\nplt.title('Distribution of Y_train')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:03.798336Z","iopub.execute_input":"2023-04-28T08:40:03.798728Z","iopub.status.idle":"2023-04-28T08:40:04.375336Z","shell.execute_reply.started":"2023-04-28T08:40:03.798691Z","shell.execute_reply":"2023-04-28T08:40:04.373932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"See 2.3 for data exploration after feature extraction","metadata":{}},{"cell_type":"markdown","source":"### 2.3 Data Preparation","metadata":{}},{"cell_type":"markdown","source":"#### 2.3.1 Feature extraction\nfor feature extraction, I tried a number of features to select the best from. I divided them in basic features (with features like RMS, standard deviation etc), percentile features, features using scipy stats, features I used earlier in a different project and some MFCC features. I got inspriation for a number of features from this notebook: https://www.kaggle.com/code/abhishek/quite-a-few-features-1-51/notebook","metadata":{}},{"cell_type":"code","source":"#feature extraction\nimport librosa\nimport scipy as sp\nimport time\n\n\n#root mean square\ndef RMS(x):\n return np.apply_along_axis(lambda x: librosa.feature.rms(y=x).mean(), axis=1, arr=x)\n\n\n#autocorrelation\ndef acl(x):\n ac = librosa.autocorrelate(y=x)\n return (ac/ac[0]).mean()\n \n#autocorrelation\ndef autocorrelation(x):\n ac = np.apply_along_axis(acl, axis=1, arr=x)\n return ac\n#std\ndef std(x):\n return np.std(x, axis=1)\n#min\ndef minim(x):\n return np.min(x, axis=1)\n#max\ndef maxim(x):\n return np.max(x, axis=1)\n#min time\ndef min_time(x):\n return np.argmin(x, axis=1) \n#max time\ndef max_time(x):\n return np.argmax(x, axis=1)\n\n#from https://www.kaggle.com/code/abhishek/quite-a-few-features-1-51/notebook\n\n#percentile 10\ndef perc10(x):\n return np.percentile(x, 10, axis=1)\n#percentile 30\ndef perc30(x):\n return np.percentile(x, 30, axis=1)\n#percentile 50\ndef perc50(x):\n return np.percentile(x, 50, axis=1)\n#percentile 70\ndef perc70(x):\n return np.percentile(x, 70, axis=1)\n#percentile 90\ndef perc90(x):\n return np.percentile(x, 90, axis=1)\n\n\n#skew\ndef skew(x):\n return sp.stats.skew(x, axis=1)\ndef kurt(x):\n return sp.stats.kurtosis(x, axis=1)\n'''def kstat1(x):\n return sp.stats.kstat(x, 1, axis=1, )\ndef kstat2(x):\n return sp.stats.kstat(x, 2, axis=1) \ndef kstat3(x):\n return sp.stats.kstat(x, 3, axis=1)\ndef kstat4(x):\n return sp.stats.kstat(x, 4, axis=1)\ndef moment1(x):\n return sp.stats.moment(x, 1, axis=1)\ndef moment2(x):\n return sp.stats.moment(x, 2, axis=1)\ndef moment3(x):\n return sp.stats.moment(x, 3, axis=1)\ndef moment4(x):\n return sp.stats.moment(x, 4, axis=1)'''\n\n#''', kstat1(x),kstat2(x),kstat3(x),kstat4(x), moment1(x),moment2(x),moment3(x),moment4(x)'''\n\n#some features from week 4\n\n#spectral centroid\ndef spectral_centroid(x):\n return librosa.feature.spectral_centroid(y=x).mean()\n#spectral flatness\ndef spectral_flatness(x):\n return librosa.feature.spectral_flatness(y=x).mean()\n#spectral rolloff\ndef spectral_rolloff(x):\n return librosa.feature.spectral_rolloff(y=x).mean()\n#spectral bandwidth\ndef spectral_bandwidth(x):\n return librosa.feature.spectral_bandwidth(y=x).mean()\n#MFCC\ndef MFCC(x, n=3):\n n_mfcc = n\n return librosa.feature.mfcc(y=x, n_mfcc=n_mfcc).mean(axis=1)\n\n\n#all fetures\ndef feature_extraction(x):\n t0 = time.time()\n #features_basic = np.array([RMS(x), autocorrelation(x), std(x), minim(x), maxim(x), min_time(x),max_time(x)]).T\n features_basic = np.array([std(x), maxim(x), min_time(x),max_time(x)]).T\n t1 = time.time()\n print(f'features_basic: done after {t1-t0} s')\n \n #features_percentiles = np.array([perc10(x),perc30(x),perc50(x),perc70(x),perc90(x)]).T\n features_percentiles = np.array([perc10(x),perc70(x)]).T\n t2 = time.time()\n print(f'features_percentiles: done after {t2-t1} s')\n \n #features_sp = np.array([skew(x), kurt(x)]).T\n features_sp = np.array([]).T\n t3 = time.time()\n print(f'features_sp: done after {t3-t2} s')\n \n #features_week4 = np.array([np.apply_along_axis(spectral_centroid, 1, x), np.apply_along_axis(spectral_flatness, 1, x), np.apply_along_axis(spectral_rolloff, 1, x), np.apply_along_axis(spectral_bandwidth, 1, x)]).T\n features_week4 = np.array([]).T\n \n t4 = time.time()\n print(f'features_week4: done after {t4-t3} s')\n \n features_MFCC = np.apply_along_axis(MFCC, 1, x)\n t5 = time.time()\n print(f'features_MFCC: done after {t5-t4} s')\n \n #features = np.concatenate((features_basic, features_percentiles, features_sp, features_week4, features_MFCC), axis=1)\n features = np.concatenate((features_basic, features_percentiles, features_MFCC), axis=1) \n \n return features","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:04.377353Z","iopub.execute_input":"2023-04-28T08:40:04.378339Z","iopub.status.idle":"2023-04-28T08:40:04.410063Z","shell.execute_reply.started":"2023-04-28T08:40:04.378281Z","shell.execute_reply":"2023-04-28T08:40:04.408215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#1,4,9,10,12,13,14,15,16,17,18,20,22,25,27 dont do much","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:04.414554Z","iopub.execute_input":"2023-04-28T08:40:04.415180Z","iopub.status.idle":"2023-04-28T08:40:04.421937Z","shell.execute_reply.started":"2023-04-28T08:40:04.415137Z","shell.execute_reply":"2023-04-28T08:40:04.420179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(X_train.shape)\nprint(RMS(X_train[0:20]))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:04.423919Z","iopub.execute_input":"2023-04-28T08:40:04.424325Z","iopub.status.idle":"2023-04-28T08:40:18.607208Z","shell.execute_reply.started":"2023-04-28T08:40:04.424288Z","shell.execute_reply":"2023-04-28T08:40:18.605329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#make features\n\nX_train_features = feature_extraction(X_train)\nX_test_features = feature_extraction(X_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:40:18.609467Z","iopub.execute_input":"2023-04-28T08:40:18.610235Z","iopub.status.idle":"2023-04-28T08:43:20.413282Z","shell.execute_reply.started":"2023-04-28T08:40:18.610192Z","shell.execute_reply":"2023-04-28T08:43:20.409899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(X_train_features.shape)\nprint(X_test_features.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:20.418618Z","iopub.execute_input":"2023-04-28T08:43:20.423267Z","iopub.status.idle":"2023-04-28T08:43:20.436326Z","shell.execute_reply.started":"2023-04-28T08:43:20.423160Z","shell.execute_reply":"2023-04-28T08:43:20.434627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 2.3.2 Feature normalization\nIn this section, the features get normalized. The normalization is defined by only X_train, and performed on X_train and X_test.","metadata":{}},{"cell_type":"code","source":"#normalize input\n\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(X_train_features)\nX_train_normalized = scaler.transform(X_train_features)\nX_test_normalized = scaler.transform(X_test_features)\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:20.438686Z","iopub.execute_input":"2023-04-28T08:43:20.441054Z","iopub.status.idle":"2023-04-28T08:43:20.462374Z","shell.execute_reply.started":"2023-04-28T08:43:20.440974Z","shell.execute_reply":"2023-04-28T08:43:20.460021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 2.3.3 exploring the features\nTo discover what features were important and what features could better be left behind, it is useful to plot the features. I tried different combinations of features to plot, and listed my observations. I used these observations to remove some features in section 2.3.1 to simplify the model.","metadata":{}},{"cell_type":"code","source":"# Plotting X_train_normalized\n\nprint(X_train_normalized.shape)\nfig = plt.figure(figsize=(10,10))\nax1 = fig.add_subplot(2,2,1, projection='3d')\nax1.scatter(X_train_normalized[:, 0], X_train_normalized[:, 1], X_train_normalized[:, 2], c=Y_train, cmap='coolwarm')\nax1.set_xlabel('std')\n#ax1.set_xlim([-1, 1])\nax1.set_ylabel('maxim')\n#ax1.set_ylim([-1, 1])\nax1.set_zlabel('min_time')\n\nax2 = fig.add_subplot(2,2,2, projection='3d')\nax2.scatter(X_train_normalized[:, 3], X_train_normalized[:, 4], X_train_normalized[:, 5], c=Y_train, cmap='coolwarm')\nax2.set_xlabel('max_time')\nax2.set_ylabel('perc10')\nax2.set_zlabel('perc70')\n\nax3 = fig.add_subplot(2,2,3,projection='3d')\nax3.scatter(X_train_normalized[:, 6], X_train_normalized[:, 7], X_train_normalized[:, 8], c=Y_train, cmap='coolwarm')\nax3.set_xlabel('MFCC 1')\nax3.set_ylabel('MFCC 2')\nax3.set_zlabel('MFCC 3')\n\n#ax4 = fig.add_subplot(2,2,4,projection='3d')\n#ax4.scatter(X_train_normalized[:, 9], X_train_normalized[:, 10], X_train_normalized[:, 11], c=Y_train, cmap='coolwarm')\n#ax4.set_xlabel('MFCC 4')\n#ax4.set_ylabel('MFCC 5')\n#ax4.set_zlabel('MFCC 6')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:20.466298Z","iopub.execute_input":"2023-04-28T08:43:20.467702Z","iopub.status.idle":"2023-04-28T08:43:21.481196Z","shell.execute_reply.started":"2023-04-28T08:43:20.467609Z","shell.execute_reply":"2023-04-28T08:43:21.480118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"observations:\n1. maxtime and mintime are very correlated\n2. maxim and std barely have variance\n3. perc70 is either -1 or 1\n4. MFCC's are often correlated\n5. perc10 is a good indicator\n6. std and maxim are correlated","metadata":{}},{"cell_type":"markdown","source":"#### 2.3.4 PCA\nTo better understand the data, PCA is performed to see if the features are very correlated. The PCA result could also be used as input to the model.","metadata":{}},{"cell_type":"code","source":"#PCA\n\nfrom sklearn.decomposition import PCA\npca = PCA(n_components=3)\npca = pca.fit(X_train_normalized)\n\nX_train_pca = pca.transform(X_train_normalized)\nX_test_pca = pca.transform(X_test_normalized)\n\nsingular_values = pca.singular_values_\nprint(f'singular values: {singular_values}')\n\n# Variance\nplt.bar(range(0,3), pca.explained_variance_ratio_[:3], label=\"individual var\");\nplt.step(range(0,3), np.cumsum(pca.explained_variance_ratio_[:3]),'r', label=\"cumulative var\");\nplt.xlabel('Principal component index'); plt.ylabel('explained variance ratio %');\nplt.legend()\nplt.show()\n\n# Plotting PCA\nfig = plt.figure()\nax = fig.add_subplot(projection='3d')\nax.scatter(X_train_pca[:, 0], X_train_pca[:, 1], X_train_pca[:, 2], c=Y_train, cmap='coolwarm')\nax.set_xlabel('PC1')\nax.set_ylabel('PC2')\nax.set_zlabel('PC3')\nplt.show()\n\nplt.plot(X_train_pca[:, 0], X_train_pca[:, 1], 'o')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:21.482411Z","iopub.execute_input":"2023-04-28T08:43:21.482775Z","iopub.status.idle":"2023-04-28T08:43:22.527946Z","shell.execute_reply.started":"2023-04-28T08:43:21.482740Z","shell.execute_reply":"2023-04-28T08:43:22.526310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Training and Results","metadata":{}},{"cell_type":"markdown","source":"#### 3.1 linear model\nA linear regression model is chosen as a linear model for its simplicity.","metadata":{}},{"cell_type":"code","source":"X_train_input = X_train_normalized\n#X_train_input = X_train_pca\nX_test_input = X_test_normalized\n#X_test_input = X_test_pca","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:22.529773Z","iopub.execute_input":"2023-04-28T08:43:22.530175Z","iopub.status.idle":"2023-04-28T08:43:22.537209Z","shell.execute_reply.started":"2023-04-28T08:43:22.530137Z","shell.execute_reply":"2023-04-28T08:43:22.535417Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#define linear model\n\ndef model(x,w):\n a = w[0] + x @ w[1:]\n return (a.T)\n\nn_features = X_train_input.shape[1]\nw = np.random.randn(n_features+1)\nprint(w)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:22.539113Z","iopub.execute_input":"2023-04-28T08:43:22.539631Z","iopub.status.idle":"2023-04-28T08:43:22.550619Z","shell.execute_reply.started":"2023-04-28T08:43:22.539591Z","shell.execute_reply":"2023-04-28T08:43:22.549691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#define cost function\nimport sklearn\n\nprint(model(X_train_input,w).shape)\n\ndef huber_loss(y, y_pred):\n residual = np.abs(y - y_pred)\n is_outlier = residual > 1\n quadratic = np.minimum(residual, 1)\n linear = residual - quadratic\n return np.mean(0.5 * quadratic**2 + 1 * linear)\n\ndef cost_function(y,y_pred):\n #MSE loss\n return np.mean((y-y_pred)**2)\n \n #MAE loss\n #return np.mean(np.abs(y-y_pred))\n \n #huber loss\n #return huber_loss(y, y_pred)\n\nY_train_pred = model(X_train_input,w)\n\nprint(cost_function(Y_train, Y_train_pred))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T09:53:49.229671Z","iopub.execute_input":"2023-04-28T09:53:49.230086Z","iopub.status.idle":"2023-04-28T09:53:49.242482Z","shell.execute_reply.started":"2023-04-28T09:53:49.230052Z","shell.execute_reply":"2023-04-28T09:53:49.241198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#define gradient descent\nfrom autograd import grad\n\ndef dummy(w):\n return cost_function(Y_train, model(X_train_input, w))\n\ndef gradient_descent(cost_function=cost_function, alpha=0.01, n_epochs=10, w=w):\n gradient = grad(dummy)\n \n Y_train_pred = model(X_train_input,w)\n Y_test_pred = model(X_test_input,w)\n weight_history = [w]\n cost_history = [cost_function(Y_train, Y_train_pred)]\n cost_history_test = [cost_function(Y_test, Y_test_pred)]\n \n for epoch in range(n_epochs):\n print(f'epoch {epoch}', end='\\r')\n Y_train_pred = model(X_train_input,w)\n Y_test_pred = model(X_test_input,w)\n \n grad_eval = gradient(w)\n \n w = alpha * -grad_eval + w\n \n weight_history.append(w)\n cost_history.append(cost_function(Y_train, Y_train_pred))\n cost_history_test.append(cost_function(Y_test, Y_test_pred))\n print(f'MSE: {np.mean((Y_test-Y_test_pred)**2)}')\n print(f'MAE: {np.mean(np.abs(Y_test-Y_test_pred))}')\n print(f'huber: {huber_loss(Y_test, Y_test_pred)}')\n plt.plot(cost_history)\n plt.plot(cost_history_test)\n plt.yscale('log')\n plt.show()\n \n \n \n\ngradient_descent(n_epochs=2000)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T09:53:49.730696Z","iopub.execute_input":"2023-04-28T09:53:49.731504Z","iopub.status.idle":"2023-04-28T09:53:54.078749Z","shell.execute_reply.started":"2023-04-28T09:53:49.731463Z","shell.execute_reply":"2023-04-28T09:53:54.077113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### comparison of cost functions\nafter 2000 epochs, alpha = 0.01
\n#### cost function: MSE\nMSE: 7.55133497637975
\nMAE: 2.2116551821267816
\nhuber: 1.747647870461226
\n\n#### cost function: MAE\nMSE: 7.794476006891011
\nMAE: 2.200750433336176
\nhuber: 1.7409732308642503
\n\n\n#### cost function: Huber\nMSE: 7.794240121126477
\nMAE: 2.203431468955108
\nhuber: 1.7430145533145445
\n
\nIt is difficult to compare the effectiveness of cost functions, because the error is defined by the cost function itself. However, if we calculate each cost function value for each cost function, we might get a better idea. It seems that MSE works the best for MSE (obviously), but scores better on MSA then MSA scores on MSE. For this reason, MSE is chosen.","metadata":{}},{"cell_type":"code","source":"#take a look at predicted vs true values of time before earthquake\nY_train_pred = model(X_train_input,w)\nfor i in range(10):\n print(f'Y value: {Y_train[i]}')\n print(f'Y predicted: {Y_train_pred[i]}')\n print()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T09:54:23.547957Z","iopub.execute_input":"2023-04-28T09:54:23.548444Z","iopub.status.idle":"2023-04-28T09:54:23.560727Z","shell.execute_reply.started":"2023-04-28T09:54:23.548405Z","shell.execute_reply":"2023-04-28T09:54:23.558168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get importance\nimportance = w\n\n#plot feature importance\nplt.bar([x for x in range(len(importance))], importance)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T09:54:35.377132Z","iopub.execute_input":"2023-04-28T09:54:35.378258Z","iopub.status.idle":"2023-04-28T09:54:35.599488Z","shell.execute_reply.started":"2023-04-28T09:54:35.378205Z","shell.execute_reply":"2023-04-28T09:54:35.597875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3.2: non-linear model\nA random forest regression model is chosen.","metadata":{}},{"cell_type":"code","source":"#define non-linear model\n\nfrom sklearn.ensemble import RandomForestRegressor\nmodel_RF = RandomForestRegressor()\nmodel_RF.fit(X_train_input, Y_train)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:24.578245Z","iopub.execute_input":"2023-04-28T08:43:24.578567Z","iopub.status.idle":"2023-04-28T08:43:27.021548Z","shell.execute_reply.started":"2023-04-28T08:43:24.578536Z","shell.execute_reply":"2023-04-28T08:43:27.020031Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#take a look at predicted vs true values of time before earthquake\nY_train_pred = model_RF.predict(X_train_input)\nfor i in range(10):\n print(f'Y value: {Y_train[i]}')\n print(f'Y predicted: {Y_train_pred[i]}')\n print()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:27.023707Z","iopub.execute_input":"2023-04-28T08:43:27.024470Z","iopub.status.idle":"2023-04-28T08:43:27.125293Z","shell.execute_reply.started":"2023-04-28T08:43:27.024417Z","shell.execute_reply":"2023-04-28T08:43:27.123545Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3.3 visualisation\nTo look at the results of both models and to compare them, the true and predicted y-values are plotted for both.","metadata":{}},{"cell_type":"code","source":"Y_pred_train_LR = model(X_train_input, w)\nY_pred_train_RF = model_RF.predict(X_train_input)\n\nY_pred_test_LR = model(X_test_input, w)\nY_pred_test_RF = model_RF.predict(X_test_input)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:27.127061Z","iopub.execute_input":"2023-04-28T08:43:27.127577Z","iopub.status.idle":"2023-04-28T08:43:27.327035Z","shell.execute_reply.started":"2023-04-28T08:43:27.127542Z","shell.execute_reply":"2023-04-28T08:43:27.325695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import r2_score\n\nr2_LR = r2_score(Y_test, Y_pred_test_LR)\nr2_RF = r2_score(Y_test, Y_pred_test_RF)\nprint(\"R-squared score linear regression model:\", r2_LR)\nprint(\"R-squared score random forest model:\", r2_RF)\n\nplt.scatter(Y_pred_test_LR, Y_test)\nplt.plot(np.linspace(0,10,100),np.linspace(0,10,100), c='red')\nplt.xlabel('Y_pred_test_LR')\nplt.ylabel('Y_test')\nplt.title('Linear regression, R2 = '+ str(r2_LR))\nplt.show()\n\nplt.scatter(Y_pred_test_RF, Y_test)\nplt.plot(np.linspace(0,10,100),np.linspace(0,10,100), c='red')\nplt.xlabel('Y_pred_test_RF')\nplt.ylabel('Y_test')\nplt.title('Random Forest, R2 = '+ str(r2_RF))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T11:18:32.029105Z","iopub.execute_input":"2023-04-28T11:18:32.029675Z","iopub.status.idle":"2023-04-28T11:18:32.551159Z","shell.execute_reply.started":"2023-04-28T11:18:32.029615Z","shell.execute_reply":"2023-04-28T11:18:32.549630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from mpl_toolkits import mplot3d\n\nfig = plt.figure(figsize=(15,15))\nfig.suptitle('Linear regression model')\nax1 = fig.add_subplot(2,4,2, projection='3d')\nscatter1 = ax1.scatter(X_train_input[:, 0], X_train_input[:, 1], X_train_input[:, 2], c=Y_pred_train_LR, cmap='coolwarm', vmin=0, vmax=15)\nax1.set_xlabel('feature 0')\n#ax1.set_xlim([-1, 1])\nax1.set_ylabel('feature 1')\n#ax1.set_ylim([-1, 1])\nax1.set_zlabel('feature 2')\nfig.colorbar(scatter1, fraction=0.046, pad=0.04)\nax1.set_title('Y_train_pred')\n\nax2 = fig.add_subplot(2,4,1, projection='3d')\nscatter2 = ax2.scatter(X_train_input[:, 0], X_train_input[:, 1], X_train_input[:, 2], c=Y_train, cmap='coolwarm', vmin=0, vmax=15)\nax2.set_xlabel('feature 0')\n#ax2.set_xlim([-1, 1])\nax2.set_ylabel('feature 1')\n#ax2.set_ylim([-1, 1])\nax2.set_zlabel('feature 2')\nfig.colorbar(scatter2, fraction=0.046, pad=0.04)\nax2.set_title('Y_train')\n\nax3 = fig.add_subplot(2,4,6, projection='3d')\nscatter3 = ax3.scatter(X_train_input[:, 7], X_train_input[:, 8], X_train_input[:, 6], c=Y_pred_train_LR, cmap='coolwarm', vmin=0, vmax=15)\nax3.set_xlabel('feature 3')\nax3.set_xlim([-1, 1])\nax3.set_ylabel('feature 4')\nax3.set_ylim([-1, 1])\nax3.set_zlabel('feature 5')\nfig.colorbar(scatter3, fraction=0.046, pad=0.04)\nax3.set_title('Y_train_pred')\n\nax4 = fig.add_subplot(2,4,5, projection='3d')\nscatter4 = ax4.scatter(X_train_input[:, 7], X_train_input[:, 8], X_train_input[:, 6], c=Y_train, cmap='coolwarm', vmin=0, vmax=15)\nax4.set_xlabel('feature 3')\nax4.set_xlim([-1, 1])\nax4.set_ylabel('feature 4')\nax4.set_ylim([-1, 1])\nax4.set_zlabel('feature 5')\nfig.colorbar(scatter4, fraction=0.046, pad=0.04)\nax4.set_title('Y_train')\n\nax5 = fig.add_subplot(2,4,4, projection='3d')\nscatter5 = ax5.scatter(X_test_input[:, 0], X_test_input[:, 1], X_test_input[:, 2], c=Y_pred_test_LR, cmap='coolwarm', vmin=0, vmax=15)\nax5.set_xlabel('feature 0')\n#ax2.set_xlim([-1, 1])\nax5.set_ylabel('feature 1')\n#ax2.set_ylim([-1, 1])\nax5.set_zlabel('feature 2')\nfig.colorbar(scatter5, fraction=0.046, pad=0.04)\nax5.set_title('Y_test_pred')\n\nax6 = fig.add_subplot(2,4,8, projection='3d')\nscatter6 = ax6.scatter(X_test_input[:, 7], X_test_input[:, 8], X_test_input[:, 6], c=Y_pred_test_LR, cmap='coolwarm', vmin=0, vmax=15)\nax6.set_xlabel('feature 7')\n#ax2.set_xlim([-1, 1])\nax6.set_ylabel('feature 8')\n#ax2.set_ylim([-1, 1])\nax6.set_zlabel('feature 6')\nfig.colorbar(scatter6, fraction=0.046, pad=0.04)\nax6.set_title('Y_test_pred')\n\nax7 = fig.add_subplot(2,4,3, projection='3d')\nscatter7 = ax7.scatter(X_test_input[:, 0], X_test_input[:, 1], X_test_input[:, 2], c=Y_test, cmap='coolwarm', vmin=0, vmax=15)\nax7.set_xlabel('feature 0')\n#ax2.set_xlim([-1, 1])\nax7.set_ylabel('feature 1')\n#ax2.set_ylim([-1, 1])\nax7.set_zlabel('feature 2')\nfig.colorbar(scatter7, fraction=0.046, pad=0.04)\nax7.set_title('Y_test')\n\nax8 = fig.add_subplot(2,4,7, projection='3d')\nscatter8 = ax8.scatter(X_test_input[:, 7], X_test_input[:, 8], X_test_input[:, 6], c=Y_test, cmap='coolwarm', vmin=0, vmax=15)\nax8.set_xlabel('feature 7')\n#ax2.set_xlim([-1, 1])\nax8.set_ylabel('feature 8')\n#ax2.set_ylim([-1, 1])\nax8.set_zlabel('feature 6')\nfig.colorbar(scatter8, fraction=0.046, pad=0.04)\nax8.set_title('Y_test')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:27.328606Z","iopub.execute_input":"2023-04-28T08:43:27.328946Z","iopub.status.idle":"2023-04-28T08:43:30.844628Z","shell.execute_reply.started":"2023-04-28T08:43:27.328916Z","shell.execute_reply":"2023-04-28T08:43:30.843554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(15,15))\nfig.suptitle('Random Forest regression model')\nax1 = fig.add_subplot(2,4,2, projection='3d')\nscatter1 = ax1.scatter(X_train_input[:, 0], X_train_input[:, 1], X_train_input[:, 2], c=Y_pred_train_RF, cmap='coolwarm', vmin=0, vmax=15)\nax1.set_xlabel('feature 0')\nax1.set_xlim([-0.5, 2])\nax1.set_ylabel('feature 1')\nax1.set_ylim([-0.5, 2])\nax1.set_zlabel('feature 2')\nfig.colorbar(scatter1, fraction=0.046, pad=0.04)\nax1.set_title('Y_train_pred')\n\nax2 = fig.add_subplot(2,4,1, projection='3d')\nscatter2 = ax2.scatter(X_train_input[:, 0], X_train_input[:, 1], X_train_input[:, 2], c=Y_train, cmap='coolwarm', vmin=0, vmax=15)\nax2.set_xlabel('feature 0')\nax2.set_xlim([-0.5, 2])\nax2.set_ylabel('feature 1')\nax2.set_ylim([-0.5, 2])\nax2.set_zlabel('feature 2')\nfig.colorbar(scatter2, fraction=0.046, pad=0.04)\nax2.set_title('Y_train')\n\nax3 = fig.add_subplot(2,4,6, projection='3d')\nscatter3 = ax3.scatter(X_train_input[:, 7], X_train_input[:, 8], X_train_input[:, 6], c=Y_pred_train_RF, cmap='coolwarm', vmin=0, vmax=15)\nax3.set_xlabel('feature 3')\nax3.set_xlim([-1, 1])\nax3.set_ylabel('feature 4')\nax3.set_ylim([-1, 1])\nax3.set_zlabel('feature 5')\nfig.colorbar(scatter3, fraction=0.046, pad=0.04)\nax3.set_title('Y_train_pred')\n\nax4 = fig.add_subplot(2,4,5, projection='3d')\nscatter4 = ax4.scatter(X_train_input[:, 7], X_train_input[:, 8], X_train_input[:, 6], c=Y_train, cmap='coolwarm', vmin=0, vmax=15)\nax4.set_xlabel('feature 3')\nax4.set_xlim([-1, 1])\nax4.set_ylabel('feature 4')\nax4.set_ylim([-1, 1])\nax4.set_zlabel('feature 5')\nfig.colorbar(scatter4, fraction=0.046, pad=0.04)\nax4.set_title('Y_train')\n\nax5 = fig.add_subplot(2,4,4, projection='3d')\nscatter5 = ax5.scatter(X_test_input[:, 0], X_test_input[:, 1], X_test_input[:, 2], c=Y_pred_test_RF, cmap='coolwarm', vmin=0, vmax=15)\nax5.set_xlabel('feature 0')\nax5.set_xlim([-0.5, 2])\nax5.set_ylabel('feature 1')\nax5.set_ylim([-0.5, 2])\nax5.set_zlabel('feature 2')\nfig.colorbar(scatter5, fraction=0.046, pad=0.04)\nax5.set_title('Y_test_pred')\n\nax6 = fig.add_subplot(2,4,8, projection='3d')\nscatter6 = ax6.scatter(X_test_input[:, 7], X_test_input[:, 8], X_test_input[:, 6], c=Y_pred_test_RF, cmap='coolwarm', vmin=0, vmax=15)\nax6.set_xlabel('feature 7')\n#ax2.set_xlim([-1, 1])\nax6.set_ylabel('feature 8')\n#ax2.set_ylim([-1, 1])\nax6.set_zlabel('feature 6')\nfig.colorbar(scatter6, fraction=0.046, pad=0.04)\nax6.set_title('Y_test_pred')\n\nax7 = fig.add_subplot(2,4,3, projection='3d')\nscatter7 = ax7.scatter(X_test_input[:, 0], X_test_input[:, 1], X_test_input[:, 2], c=Y_test, cmap='coolwarm', vmin=0, vmax=15)\nax7.set_xlabel('feature 0')\nax7.set_xlim([-0.5, 2])\nax7.set_ylabel('feature 1')\nax7.set_ylim([-0.5, 2])\nax7.set_zlabel('feature 2')\nfig.colorbar(scatter7, fraction=0.046, pad=0.04)\nax7.set_title('Y_test')\n\nax8 = fig.add_subplot(2,4,7, projection='3d')\nscatter8 = ax8.scatter(X_test_input[:, 7], X_test_input[:, 8], X_test_input[:, 6], c=Y_test, cmap='coolwarm', vmin=0, vmax=15)\nax8.set_xlabel('feature 7')\n#ax2.set_xlim([-1, 1])\nax8.set_ylabel('feature 8')\n#ax2.set_ylim([-1, 1])\nax8.set_zlabel('feature 6')\nfig.colorbar(scatter8, fraction=0.046, pad=0.04)\nax8.set_title('Y_test')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:30.846255Z","iopub.execute_input":"2023-04-28T08:43:30.847816Z","iopub.status.idle":"2023-04-28T08:43:34.599077Z","shell.execute_reply.started":"2023-04-28T08:43:30.847763Z","shell.execute_reply":"2023-04-28T08:43:34.598074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is clearly visible that the Random Forest works much better. Not only from the R2 values, but also the plots.
\nThe linear model seems to play very save and to always output a value close to the average. This is probably because it is not complex enough, and because of outliers. It is, however, very computationally efficient.
\nThe random forest has much less difficulty with outliers, and can capture non-linear relationships.","metadata":{}},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion\n\nTo get the best earthquake predictions, I made the following steps:\n1. Data exploration\n* I plotted some samples of X_train to see how the data looks like. With this information, I was able to think of some features that could describe them.\n* I plotted the distribution of Y_train to get an idea of the results that the models were supposed to give.\n2. Data preprocessing\n* I cut the training data in pieces of 150,000 elements, because they resemble the test data this way. I could have used another approach, in which I did not cut up the data, but choose n random elements to start a 150,000 element long snippet with. This way, there is overlapping data between snippets, but you are able to extract a lot more snippets from the data. I chose not to do this, because the models I wanted to use didn't need more data than generated in the process first described.\n* I split the data in testing and training data, with a 1/4 ratio.\n* I tried many different features, and looked at which of those were useful. I then selected the best features to continue with.\n* I normalized the remaining features\n* I performed PCA to check how much the features were still correlated.\n3. Training models\n* I made a Linear Regression model and a Random Forest model and trained them on the data using different cost functions. I chose the best cost function to continue with.\n* I visualised the predicted outcomes and compared them to the true outcomes. With this information, I experimented to optimise the Linear Regression model.\n
\n
\n\nStrong points of my approach include:\n* testing many features\n* testing different cost functions\n* investigating features and results by plotting them with Y values as color\n
\n\nImprovements could be:\n* using k-fold validation\n* using different models, maybe neural networks\n* diving deeper into earthquake prediction to learn new important features\n\n","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### make submission","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id')\nprint(submission.head)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:34.604768Z","iopub.execute_input":"2023-04-28T08:43:34.606288Z","iopub.status.idle":"2023-04-28T08:43:34.647267Z","shell.execute_reply.started":"2023-04-28T08:43:34.606239Z","shell.execute_reply":"2023-04-28T08:43:34.645952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"folder_path = '/kaggle/input/LANL-Earthquake-Prediction/test/'\n\n# List all CSV files in the folder\ncsv_files = [os.path.join(folder_path, f) for f in os.listdir(folder_path) if f.endswith('.csv')]\n\n# Read each CSV file into a pandas DataFrame and extract the numerical data\ndata_list = []\nfor csv_file in csv_files:\n df = pd.read_csv(csv_file)\n data = df.values # extract numerical data from DataFrame\n data_list.append(data)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:43:34.648861Z","iopub.execute_input":"2023-04-28T08:43:34.649235Z","iopub.status.idle":"2023-04-28T08:44:42.420679Z","shell.execute_reply.started":"2023-04-28T08:43:34.649199Z","shell.execute_reply":"2023-04-28T08:44:42.419116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Combine all data arrays into a single numpy array\ndata_array = np.concatenate(data_list, axis=1).T","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:44:42.423599Z","iopub.execute_input":"2023-04-28T08:44:42.424162Z","iopub.status.idle":"2023-04-28T08:44:54.483584Z","shell.execute_reply.started":"2023-04-28T08:44:42.424108Z","shell.execute_reply":"2023-04-28T08:44:54.482259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_array = data_array.astype(float)\nprint(data_array.dtype)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:44:54.485406Z","iopub.execute_input":"2023-04-28T08:44:54.485988Z","iopub.status.idle":"2023-04-28T08:44:56.532342Z","shell.execute_reply.started":"2023-04-28T08:44:54.485937Z","shell.execute_reply":"2023-04-28T08:44:56.530889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_submission_features = feature_extraction(data_array)\nX_submission_normalized = scaler.transform(X_submission_features)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:44:56.534277Z","iopub.execute_input":"2023-04-28T08:44:56.534827Z","iopub.status.idle":"2023-04-28T08:48:28.893260Z","shell.execute_reply.started":"2023-04-28T08:44:56.534775Z","shell.execute_reply":"2023-04-28T08:48:28.891260Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# use our model to make predictions\nY_submission_pred = model_RF.predict(X_submission_normalized)\n\nprint(Y_submission_pred[:10])","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:48:28.896000Z","iopub.execute_input":"2023-04-28T08:48:28.904319Z","iopub.status.idle":"2023-04-28T08:48:29.028586Z","shell.execute_reply.started":"2023-04-28T08:48:28.904193Z","shell.execute_reply":"2023-04-28T08:48:29.027145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission['time_to_failure'] = Y_submission_pred\nprint(submission)\nsubmission.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T08:48:29.030798Z","iopub.execute_input":"2023-04-28T08:48:29.031752Z","iopub.status.idle":"2023-04-28T08:48:29.064275Z","shell.execute_reply.started":"2023-04-28T08:48:29.031698Z","shell.execute_reply":"2023-04-28T08:48:29.062835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}