{"cells":[{"metadata":{"_uuid":"382b29be994a1e723b0bfb037f3ca89e2630750c"},"cell_type":"markdown","source":"# Earthquake Prediction: Linear Regression on 180 features."},{"metadata":{"_uuid":"82c29fc5ed1fabcd3351949d7e75637b08598c2f"},"cell_type":"markdown","source":"In this notebook, we will perform linear regression on 180 features generated from the raw data. We won't go into detail on the feature extraction itself, as there are enough other kernels that do this already. In stead, this notebook tries to provide a simple scaffolding for future more advanced notebooks where the Linear Regression can be swapped for more advanced models."},{"metadata":{"_uuid":"c22d211583c33e5fc3f52c9c563ed57ce0d3d63c"},"cell_type":"markdown","source":"## Imports"},{"metadata":{"trusted":false,"_uuid":"d0ea67d29cd1c7345888ccf311744fcd0c04dc7a"},"cell_type":"code","source":"# numerical computation\nimport numpy as np\n\n# dataframes\nimport pandas as pd\n\n# visualization\n# we could do this with matplotlib, but I wanted to try\n# something new... Do not fear: only two visualizations ;)\n# altair is a very nice plotting library by the way!\nimport altair as alt  # plots\nalt.renderers.enable(\"kaggle\")\nfrom IPython.display import display  # pretty display of e.g. dataframes\n\n# progress bars\nfrom tqdm import tqdm_notebook as tqdm\n\n# simple models and preprocessing from scikit learn\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\n\n# some constants:\nnum_lines = 629_145_480  # total number of lines in the CSV file\nnum_lines_per_segment = 150_000  # number of lines in each test segment","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6408f23130f28a13b1e255a314fec4a336936ab3"},"cell_type":"markdown","source":"## Feature extraction"},{"metadata":{"_uuid":"953385f2afa514222cf83552cc9fe5d265dafcae"},"cell_type":"markdown","source":"Given a segment of data (usually 150000 elements), the following function will construct a single vector of 180 features. This vector is stored as a DataFrame with 180 named columns and one single index (given as an argument to the function).\n\nOther kernels with more information about feature extraction are for example:\n * [artgor/earthquakes-fe-more-features-and-samples](https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples)\n * [andrekos/basic-feature-benchmark-with-quantiles](https://www.kaggle.com/andrekos/basic-feature-benchmark-with-quantiles)\n * [wimwim/rolling-quantiles](https://www.kaggle.com/wimwim/rolling-quantiles)"},{"metadata":{"trusted":false,"_uuid":"d7bf91461a79a6ebb89967f2cb86254acca09a7d"},"cell_type":"code","source":"def create_features_for_segment(idx, raw_data):\n    \"\"\" create features for a segment of raw data \n    \n    Args:\n        idx: index to save the features under\n        raw_data: the raw data segment to calculate the features for\n        \n    Returns:\n        features: a pandas DataFrame with 180 feature columns and a single\n            index specified by `idx`.\n    \n    \"\"\"\n    data = pd.DataFrame(index=[idx])\n    data.loc[idx, \"data\"] = raw_data.iloc[-1].item()\n    data.loc[idx, \"mean\"] = raw_data.mean().item()\n    data.loc[idx, \"std\"] = raw_data.std().item()\n    data.loc[idx, \"max\"] = raw_data.max().item()\n    data.loc[idx, \"min\"] = raw_data.min().item()\n    data.loc[idx, \"mad\"] = raw_data.mad().item()\n    data.loc[idx, \"kurt\"] = raw_data.kurtosis().item()\n    data.loc[idx, \"skew\"] = raw_data.skew().item()\n    data.loc[idx, \"median\"] = raw_data.median().item()\n    data.loc[idx, \"q01\"] = np.quantile(raw_data, 0.01)\n    data.loc[idx, \"q05\"] = np.quantile(raw_data, 0.05)\n    data.loc[idx, \"q95\"] = np.quantile(raw_data, 0.95)\n    data.loc[idx, \"q99\"] = np.quantile(raw_data, 0.99)\n    data.loc[idx, \"iqr\"] = np.subtract(*np.percentile(raw_data, [75, 25]))\n    data.loc[idx, \"abs_mean\"] = raw_data.abs().mean().item()\n    data.loc[idx, \"abs_std\"] = raw_data.abs().std().item()\n    data.loc[idx, \"abs_max\"] = raw_data.abs().max().item()\n    data.loc[idx, \"abs_min\"] = raw_data.abs().min().item()\n    data.loc[idx, \"abs_mad\"] = raw_data.abs().mad().item()\n    data.loc[idx, \"abs_kurt\"] = raw_data.abs().kurtosis().item()\n    data.loc[idx, \"abs_skew\"] = raw_data.abs().skew().item()\n    data.loc[idx, \"abs_median\"] = raw_data.abs().median().item()\n    data.loc[idx, \"abs_q01\"] = np.quantile(raw_data.abs(), 0.01)\n    data.loc[idx, \"abs_q05\"] = np.quantile(raw_data.abs(), 0.05)\n    data.loc[idx, \"abs_q95\"] = np.quantile(raw_data.abs(), 0.95)\n    data.loc[idx, \"abs_q99\"] = np.quantile(raw_data.abs(), 0.99)\n    data.loc[idx, \"abs_iqr\"] = np.subtract(*np.percentile(raw_data.abs(), [75, 25]))\n\n    for window in [10, 100, 1000]:\n\n        data_roll_mean = raw_data.rolling(window).mean().dropna()\n        data.loc[idx, f\"mean_mean_{window}\"] = data_roll_mean.mean().item()\n        data.loc[idx, f\"std_mean_{window}\"] = data_roll_mean.std().item()\n        data.loc[idx, f\"max_mean_{window}\"] = data_roll_mean.max().item()\n        data.loc[idx, f\"min_mean_{window}\"] = data_roll_mean.min().item()\n        data.loc[idx, f\"mad_mean_{window}\"] = data_roll_mean.mad().item()\n        data.loc[idx, f\"kurt_mean_{window}\"] = data_roll_mean.kurtosis().item()\n        data.loc[idx, f\"skew_mean_{window}\"] = data_roll_mean.skew().item()\n        data.loc[idx, f\"median_mean_{window}\"] = data_roll_mean.median().item()\n        data.loc[idx, f\"q01_mean_{window}\"] = np.quantile(data_roll_mean, 0.01)\n        data.loc[idx, f\"q05_mean_{window}\"] = np.quantile(data_roll_mean, 0.05)\n        data.loc[idx, f\"q95_mean_{window}\"] = np.quantile(data_roll_mean, 0.95)\n        data.loc[idx, f\"q99_mean_{window}\"] = np.quantile(data_roll_mean, 0.99)\n        data.loc[idx, f\"iqr_mean_{window}\"] = np.subtract(\n            *np.percentile(data_roll_mean, [75, 25])\n        )\n        data.loc[idx, f\"abs_mean_mean_{window}\"] = data_roll_mean.abs().mean().item()\n        data.loc[idx, f\"abs_std_mean_{window}\"] = data_roll_mean.abs().std().item()\n        data.loc[idx, f\"abs_max_mean_{window}\"] = data_roll_mean.abs().max().item()\n        data.loc[idx, f\"abs_min_mean_{window}\"] = data_roll_mean.abs().min().item()\n        data.loc[idx, f\"abs_mad_mean_{window}\"] = data_roll_mean.abs().mad().item()\n        data.loc[idx, f\"abs_kurt_mean_{window}\"] = (\n            data_roll_mean.abs().kurtosis().item()\n        )\n        data.loc[idx, f\"abs_skew_mean_{window}\"] = data_roll_mean.abs().skew().item()\n        data.loc[idx, f\"abs_median_mean_{window}\"] = (\n            data_roll_mean.abs().median().item()\n        )\n        data.loc[idx, f\"abs_q01_mean_{window}\"] = np.quantile(\n            data_roll_mean.abs(), 0.01\n        )\n        data.loc[idx, f\"abs_q05_mean_{window}\"] = np.quantile(\n            data_roll_mean.abs(), 0.05\n        )\n        data.loc[idx, f\"abs_q95_mean_{window}\"] = np.quantile(\n            data_roll_mean.abs(), 0.95\n        )\n        data.loc[idx, f\"abs_q99_mean_{window}\"] = np.quantile(\n            data_roll_mean.abs(), 0.99\n        )\n        data.loc[idx, f\"abs_iqr_mean_{window}\"] = np.subtract(\n            *np.percentile(data_roll_mean.abs(), [75, 25])\n        )\n\n        data_roll_std = raw_data.rolling(window).std().dropna()\n        data.loc[idx, f\"mean_std_{window}\"] = data_roll_std.mean().item()\n        data.loc[idx, f\"std_std_{window}\"] = data_roll_std.std().item()\n        data.loc[idx, f\"max_std_{window}\"] = data_roll_std.max().item()\n        data.loc[idx, f\"min_std_{window}\"] = data_roll_std.min().item()\n        data.loc[idx, f\"mad_std_{window}\"] = data_roll_std.mad().item()\n        data.loc[idx, f\"kurt_std_{window}\"] = data_roll_std.kurtosis().item()\n        data.loc[idx, f\"skew_std_{window}\"] = data_roll_std.skew().item()\n        data.loc[idx, f\"median_std_{window}\"] = data_roll_std.median().item()\n        data.loc[idx, f\"q01_std_{window}\"] = np.quantile(data_roll_mean, 0.01)\n        data.loc[idx, f\"q05_std_{window}\"] = np.quantile(data_roll_mean, 0.05)\n        data.loc[idx, f\"q95_std_{window}\"] = np.quantile(data_roll_mean, 0.95)\n        data.loc[idx, f\"q99_std_{window}\"] = np.quantile(data_roll_mean, 0.99)\n        data.loc[idx, f\"iqr_std_{window}\"] = np.subtract(\n            *np.percentile(data_roll_std, [75, 25])\n        )\n        data.loc[idx, f\"abs_mean_std_{window}\"] = data_roll_std.abs().mean().item()\n        data.loc[idx, f\"abs_std_std_{window}\"] = data_roll_std.abs().std().item()\n        data.loc[idx, f\"abs_max_std_{window}\"] = data_roll_std.abs().max().item()\n        data.loc[idx, f\"abs_min_std_{window}\"] = data_roll_std.abs().min().item()\n        data.loc[idx, f\"abs_mad_std_{window}\"] = data_roll_std.abs().mad().item()\n        data.loc[idx, f\"abs_kurt_std_{window}\"] = data_roll_std.abs().kurtosis().item()\n        data.loc[idx, f\"abs_skew_std_{window}\"] = data_roll_std.abs().skew().item()\n        data.loc[idx, f\"abs_median_std_{window}\"] = data_roll_std.abs().median().item()\n        data.loc[idx, f\"abs_q01_std_{window}\"] = np.quantile(data_roll_std.abs(), 0.01)\n        data.loc[idx, f\"abs_q05_std_{window}\"] = np.quantile(data_roll_std.abs(), 0.05)\n        data.loc[idx, f\"abs_q95_std_{window}\"] = np.quantile(data_roll_std.abs(), 0.95)\n        data.loc[idx, f\"abs_q99_std_{window}\"] = np.quantile(data_roll_std.abs(), 0.99)\n        data.loc[idx, f\"iqr_std_{window}\"] = np.subtract(\n            *np.percentile(data_roll_std, [75, 25])\n        )\n\n    return data","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fa6fb4220c1b9d46c1e404d57d51fefac4f03c6f"},"cell_type":"markdown","source":"## Data Loading"},{"metadata":{"_uuid":"f192c7080629434a181f2f7752245f9267bfe0ec"},"cell_type":"markdown","source":"The following function will load the raw data from `train.csv`:"},{"metadata":{"trusted":false,"_uuid":"fcd83f0f367b61cf804d4bc97f29932dbb9ee539"},"cell_type":"code","source":"def load_raw_train_data():\n    \"\"\" load raw train data as a dataframe \"\"\"\n    train_data = pd.read_csv(\n        filepath_or_buffer=\"../input/train.csv\",\n        dtype={\"acoustic_data\": np.int16, \"time_to_failure\": np.float32},\n    )\n    return train_data","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3f4b66a9331c8ee07da3f2c00f2629af4e0f9f93"},"cell_type":"markdown","source":"Given the raw data, the following function will calculate the 180 features for a certain number of segments and returns them as a single features dataframe. It also returns the target for each of these segments in a second dataframe. The target is defined as the last `time_to_failure` for each segment of raw data."},{"metadata":{"trusted":false,"_uuid":"8a6bafa0b60b06a9520096cc31bc26357e4e11d8"},"cell_type":"code","source":"def load_train_features_and_target():\n    \"\"\" load raw train data and transform to two (feature and target) dataframes \"\"\"\n    print(\"loading raw train data... [this takes about 2 min]\")\n    raw_data = load_raw_train_data()\n    ram_mb = raw_data.memory_usage(deep=True).sum().item() / 1024 ** 2\n    print(f\"raw train data loaded. RAM Usage: {ram_mb:.2f} MB\")\n    idxs = np.arange(num_lines_per_segment, num_lines, num_lines_per_segment // 2)\n    target_values = np.zeros((idxs.shape[0], 1))\n    feature_values = np.zeros((idxs.shape[0], 180))\n    print(\"transforming raw data into feature dataframe. This takes about 30 min...\")\n    for i, idx in enumerate(tqdm(idxs)):\n        segment = raw_data.iloc[idx - num_lines_per_segment + 1 : idx + 1]\n        target_values[i] = segment.time_to_failure.values[-1:]\n        segment = segment[[\"acoustic_data\"]]\n        feature_row = create_features_for_segment(idx, segment)\n        feature_values[i, :] = feature_row.values\n    features = pd.DataFrame(\n        index=idxs, data=feature_values, columns=feature_row.columns\n    )\n    print(\"train feature dataframe created\")\n    target = pd.DataFrame(index=idxs, data=target_values, columns=[\"time_to_failure\"])\n    print(\"train target dataframe created\")\n    return features, target","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8440af55214a659bbe4d7164eb700afc6230a60c"},"cell_type":"markdown","source":"The following function will load all the raw test data segments and perform the feature extraction on them. It returns a single dataframe with the test features."},{"metadata":{"trusted":false,"_uuid":"16e01337d9ffb3665171024cbb5ff8eb2936d750"},"cell_type":"code","source":"def load_test_features():\n    \"\"\" load raw test data and transform to a feature dataframe \"\"\"\n    print(\"loading train segment ids...\")\n    seg_ids = pd.read_csv(\"../input/sample_submission.csv\", index_col=\"seg_id\").index\n    feature_values = np.zeros((seg_ids.shape[0], 180))\n    print(\"converting test segments into feature dataframes. \"\n          \"This takes a about 30 min...\")\n    for i, seg_id in enumerate(tqdm(seg_ids)):\n        segment = pd.read_csv(\n            filepath_or_buffer=f\"../input/test/{seg_id}.csv\",\n            dtype={\"acoustic_data\": np.int16, \"time_to_failure\": np.float32},\n        )\n        feature_row = create_features_for_segment(seg_id, segment)\n        feature_values[i, :] = feature_row\n    features = pd.DataFrame(\n        index=seg_ids, data=feature_values, columns=feature_row.columns\n    )\n    print(\"test feature dataframe created\")\n    return features","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0a0526920d8eb5d6c38ab52bfe8ac5938728fe18"},"cell_type":"markdown","source":"## Load Data"},{"metadata":{"_uuid":"006bb2efbf050751f5570fac33c5a4c00ab139a8"},"cell_type":"markdown","source":"With this out of the way, we can start with loading the data. We load the features and the target *[this takes about 30 min!]* and scale the former with a standard scaler from scikit-learn. The transformation of this scaler is nothing more than subtracting the mean and dividing by the standard deviation for each column."},{"metadata":{"trusted":false,"_uuid":"af288cdde16cc100f25408634247e23e0f3296e8"},"cell_type":"code","source":"# load data\nfeatures, target = load_train_features_and_target()\n\n# scale features inplace\nfeature_scaler = StandardScaler(copy=False)\nfeature_scaler.fit_transform(features)\n\nprint(\"\\n\\n\\ndata:\")\nprint(features.shape)\ndisplay(features.head())\n\nprint(\"\\n\\n\\ntarget:\")\nprint(target.shape)\ndisplay(target.head())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4aa24dc34ffb426e5639499b32e20a4d03308a0b"},"cell_type":"markdown","source":"## Train"},{"metadata":{"_uuid":"be2cf154b00d652d59d6e3e97a5292fe693d348e"},"cell_type":"markdown","source":"Now that all the data is loaded, we can go on to the training of our model. In this case it is nothing more than defining a `LinearRegression` instance of scikit-learn and fitting it to the data. However, more complex models can be easily inserted here."},{"metadata":{"trusted":false,"_uuid":"442d02253aac2dcd68fbc236256508f1378eba84"},"cell_type":"code","source":"model = LinearRegression()\nmodel.fit(features, target)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8a9e887bc2efa0e5d617979d458b2214630f4d81"},"cell_type":"markdown","source":"## Evaluation"},{"metadata":{"_uuid":"c6be8305e307c337cacfffce236ae8b526c28e57"},"cell_type":"markdown","source":"The next step after training is evaluating how well our model performs. Normally you should evaluate the model on a seperate validation dataset (usually a split of the training set). However, since we used the complete training set, we will report the training error."},{"metadata":{"_uuid":"76a7f6d68fd64fc8a0d947dacf6a581285f221ed"},"cell_type":"markdown","source":"We first define a custom function `make_prediction`, which is nothing more than a wrapper around `model.predict` returning a dataframe in stead of a numpy array:"},{"metadata":{"trusted":false,"_uuid":"c9fc2eecccf1b8adaca2c639cfe268eb62daede0"},"cell_type":"code","source":"def make_prediction(features, column_name=\"prediction\"):\n    \"\"\" custom prediction function that returns a dataframe in stead of a numpy array\"\"\"\n    prediction = pd.DataFrame(\n        index=features.index, data=model.predict(features.values), columns=[column_name]\n    )\n    return prediction","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ab631c53c54d1b0ceca179ee2a965f09e57550b8"},"cell_type":"markdown","source":"We can use this custom function to make a prediction on our data:"},{"metadata":{"trusted":false,"_uuid":"b863e906589d879ec09701cd3b06bcbe05828391"},"cell_type":"code","source":"prediction = make_prediction(features)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1705aea8b5ee6c329362eb41f805d0d03efd16ea"},"cell_type":"markdown","source":"This prediction can then be used to calculate the train error:"},{"metadata":{"trusted":false,"_uuid":"1463cf6197f697baf24d82ab5e037d278d26486b"},"cell_type":"code","source":"np.mean(np.abs(prediction.values - target.values))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c11d9d50eaf36beb18439b485f7cb3f005213912"},"cell_type":"markdown","source":"## Visualize"},{"metadata":{"_uuid":"5bc61142d87189b5204e2749296cfb01b847e190"},"cell_type":"markdown","source":"Altair stores *all* the data it recieves internally in the notebook as json. It is therefore usefull to define a slightly smaller dataframe with just the features we need for visualization to avoid bloating the size of our notebook:"},{"metadata":{"trusted":false,"_uuid":"8160ecdc4e91988bc89461d4a1274d92f73f8031"},"cell_type":"code","source":"chart_data = pd.concat([prediction, target], 1).iloc[::10].reset_index()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"52e8ff96f1b02b89165906e3f629cc49855d7c4e"},"cell_type":"markdown","source":"Using this chart data, we see a clear correlation between the predicted values and target values, which is of course good:"},{"metadata":{"trusted":false,"_uuid":"bc9ecdc5ba656beafc65bf404b3ac9721303c297"},"cell_type":"code","source":"alt.Chart(chart_data).mark_point().encode(x=\"prediction\", y=\"time_to_failure\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8ba98f5648d7dd17b6995c8196070e570aca6bf1"},"cell_type":"markdown","source":"We also see that our simple linear regression model does a fair job estimating the `time_to_failure`"},{"metadata":{"trusted":false,"_uuid":"9ef4250fcb0f4e74343b1c4c4c147a48849efbd0"},"cell_type":"code","source":"chart1 = alt.Chart(chart_data).mark_line().encode(\n    x = \"index\",\n    y = \"time_to_failure\",\n    \n)\nchart2 = alt.Chart(chart_data).mark_line().encode(\n    x = \"index\",\n    y = \"prediction\",\n    color=alt.value(\"red\")\n)\nchart1 + chart2","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fc97c9886906bde49bc06cd2d767995df729f043"},"cell_type":"markdown","source":"However, we also see that our linear regression predictor sometimes makes negative predictions for `time_to_failure`. this is ofcourse impossible. To make our final submission, we will set those values to zero."},{"metadata":{"_uuid":"bcd940c3e91fe792b74344d88f8a0d8840a9e5a2"},"cell_type":"markdown","source":"## Submission"},{"metadata":{"_uuid":"82b0834ef71e00174d886a84bb24bb113bcad713"},"cell_type":"markdown","source":"Finally, we get to the point that we can make a submission. We do this by first loading the test features *[this takes about 30min!]* after which a prediction is made, which is then stored as our `submission.csv`."},{"metadata":{"trusted":false,"_uuid":"be13accc5058eccb708f0748cf525f5fa506f81c"},"cell_type":"code","source":"test_features = load_test_features()\nfeature_scaler.transform(test_features)\nprediction = make_prediction(test_features, column_name=\"time_to_failure\")\nprediction.time_to_failure[prediction.time_to_failure < 0] = 0\nprediction.index.name = \"seg_id\"\nprediction.to_csv(\"submission.csv\")\nprediction.head()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.2"}},"nbformat":4,"nbformat_minor":1}