{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "- With this notebook we continue our analysis of `bikeshare` dataset.\n", "- In this part we will do machine learning steps in order to predict the bike rentals for the given Kaggle test set and submit it to the Kaggle\n", "\n", "![time](/images/timecycle.jpg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- In the previous post, we prepared the dataset partially \n", "- We created new columns and dropped the outliers regarding to the `\"count\"` column \n", "- Now we start by loading the dataset and continue to prepare it before applying machine learning algorithms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Outline:\n", " \n", "- About Time Data\n", "- Trigonometric functions for cyclic time data transformation\n", "- Data Split\n", "- TimeSeriesSplit for the cross validation \n", "- Sklearn Pipeline\n", "- Interpratation of RSMLE metric\n", "- Creating a custom scoring function \n", "- Kaggle submission " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import the necessary modules\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import numpy as np\n", "\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.model_selection import TimeSeriesSplit\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.model_selection import cross_val_score\n", "from sklearn.metrics import r2_score\n", "from sklearn.metrics import mean_squared_error\n", "from sklearn.metrics import make_scorer\n", "from sklearn.pipeline import Pipeline\n", "\n", "import warnings \n", "warnings.filterwarnings(\"ignore\")\n", "\n", "%matplotlib inline\n", "sns.set()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Let's remind our dataset " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Kaggle terms for the dataset:\n", "\n", "- The training set is comprised of the first `19 days` of each month, while the test set comprised of the days from `20th` to the end of the month for each month. \n", "\n", "- Predict the `total count of bikes rented during each hour` covered by the test set, using only `information available prior to the rental period`.\n", "\n", "Here is the [Kaggle link for the dataset](https://www.kaggle.com/c/bike-sharing-demand/data)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
seasonholidayworkingdayweathertempatemphumiditywindspeedcasualregisteredcountmonthweekdayhour
datetime
2011-01-01 00:00:0010019.8414.395810.031316150
2011-01-01 01:00:0010019.0213.635800.083240151
2011-01-01 02:00:0010019.0213.635800.052732152
\n", "
" ], "text/plain": [ " season holiday workingday weather temp atemp \\\n", "datetime \n", "2011-01-01 00:00:00 1 0 0 1 9.84 14.395 \n", "2011-01-01 01:00:00 1 0 0 1 9.02 13.635 \n", "2011-01-01 02:00:00 1 0 0 1 9.02 13.635 \n", "\n", " humidity windspeed casual registered count month \\\n", "datetime \n", "2011-01-01 00:00:00 81 0.0 3 13 16 1 \n", "2011-01-01 01:00:00 80 0.0 8 32 40 1 \n", "2011-01-01 02:00:00 80 0.0 5 27 32 1 \n", "\n", " weekday hour \n", "datetime \n", "2011-01-01 00:00:00 5 0 \n", "2011-01-01 01:00:00 5 1 \n", "2011-01-01 02:00:00 5 2 " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read the bike_data_inliers\n", "bike_data=pd.read_csv(\"bike_data_inliers.csv\", parse_dates=[\"datetime\"], index_col=\"datetime\")\n", "bike_data.head(3)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
holidayworkingdaytemphumiditywindspeedcountmonthweekdayhour
datetime
2011-01-01 00:00:00009.84810.016150
2011-01-01 01:00:00009.02800.040151
2011-01-01 02:00:00009.02800.032152
\n", "
" ], "text/plain": [ " holiday workingday temp humidity windspeed count \\\n", "datetime \n", "2011-01-01 00:00:00 0 0 9.84 81 0.0 16 \n", "2011-01-01 01:00:00 0 0 9.02 80 0.0 40 \n", "2011-01-01 02:00:00 0 0 9.02 80 0.0 32 \n", "\n", " month weekday hour \n", "datetime \n", "2011-01-01 00:00:00 1 5 0 \n", "2011-01-01 01:00:00 1 5 1 \n", "2011-01-01 02:00:00 1 5 2 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Filter out the casual and registered columns and \n", "## high correlated columns \"season\" (correlated with \"month\") and \"atemp\" (correlated with \"temp\") should be omitted\n", "## Also in the previous post we noticed that the data record in the \"weather\" column is not relaible\n", "## We will also drop it\n", "bike_data= bike_data.drop([\"casual\", \"registered\", \"season\", \"atemp\", \"weather\"], axis=1)\n", "bike_data.head(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## About Time Data\n", "- Since our data is a timeseries, i.e data recorded with a time dimension, we need to take into account the timeseries practices.\n", "\n", "\n", "- Firstly, we need to be aware of the `cyclic` nature of our time data\n", "\n", "\n", "- As we walked through the time frame charts in the previous post, we saw that ` month` and `hour` variables are very predictive over the number of bike rentals. \n", "\n", "\n", "- This is very resaonable as we might expect that the bike usage and rentals during `January` is less than ` May` due to the weather conditions and rentals at `3am` in the morning is less than at `14h`, during the day. \n", "\n", "\n", "- So we should use the `\"month\"` and `\"hour\"` data as features for our models. However, we have a problem of how to utilize these features. \n", "\n", "- We need to keep them **numeric** for the machine learning algorithms but if we use them without any transformation, it would not be so smart practice. \n", "\n", "To make it clearer here is a question:\n", "- Which one would you expect is more similar to the bike rentals in `January`: rentals in `December` or rentals in `May`? \n", "\n", "\n", "- Regarding the seasonal affects we can tell that number of bike rentals in `December` is more similar to rentals in `January`, however we represent `December` with number `12` and `May` with number `5`. \n", "\n", "- This is a wrong representation of the time features especially for the algorithms that use distance or the algorithms like linear regression.\n", "\n", "> As an example let's assume our model is $$ bikerentals= 15+ 2*hour $$ \n", "At 0h the number of estimated rentals = 15
\n", "At 23h the number of estimated rentals= 61
\n", "Even though there is only 1 hour difference between 0h-23h with this represantation the difference is the largest \n", "\n", "\n", "- Because of the cyclic nature of the months and hours we need to find a way of representing these features such that \n", " - the `12th` month should be closer to `1st` month than `5th` month or \n", " - `23th` hour should be closer to `1st` hour than `7th hour`. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Representation of Cyclic Time Features\n", "- We map each cyclical variable onto a circle such that the `lowest value` for that variable appears right next to the `largest value`. \n", "\n", "- We represent each variable with two components: x axis value and y axis value\n", "\n", "- We compute the `x` and `y` components of that point using `sin` and `cos` functions. \n", "\n", "![hours](/images/hours.jpg)\n", "\n", "- When we perform this transformation for the `\"month\"` variable, we also shift the values down by one such that it extends from `0 to 11`, for convenience\n", "\n", "**NOTE:** Instead of represanting the time features with trigonometric functions we can also use dummy variables. Here we will go with the trigonometric functions" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Create the sin and cos components of hour, month and weekday columns\n", "bike_data['hour_sin'] = np.sin(bike_data.index.hour*(2.*np.pi/24))\n", "bike_data['hour_cos'] = np.cos(bike_data.index.hour*(2.*np.pi/24))\n", "bike_data['month_sin'] = np.sin((bike_data.index.month-1)*(2.*np.pi/12))\n", "bike_data['month_cos'] = np.cos((bike_data.index.month-1)*(2.*np.pi/12))\n", "bike_data[\"weekday_cos\"]= np.cos((bike_data.index.weekday)*(2.*np.pi/7))\n", "bike_data['weekday_sin'] = np.sin((bike_data.index.weekday)*(2.*np.pi/7))\n", "\n", "# Now we can drop the columns \"month\", \"weekday\" and \"hour\"\n", "bike_data=bike_data.drop([\"month\", \"weekday\", \"hour\"], axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Split the data\n", "- Even though there is a given test set (without target variable) for Kaggle submission we will analyse the data as an independent project and follow our data splitting principles\n", "\n", "- At the end we will use the Kaggle test set for submission\n", "\n", "- So, after instantiating a model we will use **cross validation** for evaluating model performance and hyperparameter tuning but we still we need to keep an **hold-out set** for our final evaluation. \n", "\n", "- Since this is a timeseries dataset we must respect to the temporal order of the data. Thus, we must use only the past data to predict the future data. \n", "- So we can take the **last %5** as our hold-out data. (In the Kaggle competition the test set is the last 10 days of the months)\n", "\n", "- We need to split the data before doing transformation (incase we do) because test data points represent real-world data. \n", "- Transformation of the variables use information from the dataset like **mean** and **standard deviation**\n", "- If we transform the data before splitting we take information like mean and variance of the whole dataset thus we will be **introducing future information to the training data**. \n", "\n", "- Therefore, we should perform transformation over the training data. Then perform the same transformation on testing data as well, but this time using the parameters (like mean and variance) of the training data. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The shape of the hold-out dataset: (522, 13)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
datetimeholidayworkingdaytemphumiditywindspeedcounthour_sinhour_cosmonth_sinmonth_cosweekday_cosweekday_sin
99142012-11-16 06:00:000112.30616.00321301.0000006.123234e-17-0.8660250.5-0.900969-0.433884
99152012-11-16 07:00:000112.306111.00143670.965926-2.588190e-01-0.8660250.5-0.900969-0.433884
99162012-11-16 09:00:000113.94537.00153510.707107-7.071068e-01-0.8660250.5-0.900969-0.433884
\n", "
" ], "text/plain": [ " datetime holiday workingday temp humidity windspeed \\\n", "9914 2012-11-16 06:00:00 0 1 12.30 61 6.0032 \n", "9915 2012-11-16 07:00:00 0 1 12.30 61 11.0014 \n", "9916 2012-11-16 09:00:00 0 1 13.94 53 7.0015 \n", "\n", " count hour_sin hour_cos month_sin month_cos weekday_cos \\\n", "9914 130 1.000000 6.123234e-17 -0.866025 0.5 -0.900969 \n", "9915 367 0.965926 -2.588190e-01 -0.866025 0.5 -0.900969 \n", "9916 351 0.707107 -7.071068e-01 -0.866025 0.5 -0.900969 \n", "\n", " weekday_sin \n", "9914 -0.433884 \n", "9915 -0.433884 \n", "9916 -0.433884 " ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Find the starting indice of the last five percent\n", "last_five_percent_ind= int(len(bike_data)* 0.95)\n", "last_five_percent_ind\n", "\n", "# Create the hold-out dataset\n", "hold_out_df=bike_data.reset_index().iloc[last_five_percent_ind: ,:]\n", "\n", "print(\"The shape of the hold-out dataset:\", hold_out_df.shape)\n", "hold_out_df.head(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will do training and cross validation on the rest of the data. Lets create it" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
datetimeholidayworkingdaytemphumiditywindspeedcounthour_sinhour_cosmonth_sinmonth_cosweekday_cosweekday_sin
99112012-11-16 03:00:000112.3658.998160.7071070.707107-0.8660250.5-0.900969-0.433884
99122012-11-16 04:00:000112.3657.001550.8660250.500000-0.8660250.5-0.900969-0.433884
99132012-11-16 05:00:000112.3656.0032360.9659260.258819-0.8660250.5-0.900969-0.433884
\n", "
" ], "text/plain": [ " datetime holiday workingday temp humidity windspeed \\\n", "9911 2012-11-16 03:00:00 0 1 12.3 65 8.9981 \n", "9912 2012-11-16 04:00:00 0 1 12.3 65 7.0015 \n", "9913 2012-11-16 05:00:00 0 1 12.3 65 6.0032 \n", "\n", " count hour_sin hour_cos month_sin month_cos weekday_cos \\\n", "9911 6 0.707107 0.707107 -0.866025 0.5 -0.900969 \n", "9912 5 0.866025 0.500000 -0.866025 0.5 -0.900969 \n", "9913 36 0.965926 0.258819 -0.866025 0.5 -0.900969 \n", "\n", " weekday_sin \n", "9911 -0.433884 \n", "9912 -0.433884 \n", "9913 -0.433884 " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data= bike_data.reset_index().iloc[:last_five_percent_ind, :]\n", "data.tail(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Features and Target Variable\n", "Having splitted the hold-out dataset, now time to create the features (X) and the target (y) datasets and fit a Linear Regression model" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Target data\n", "y=data[\"count\"]\n", "\n", "# Features data\n", "X=data.drop([\"datetime\", \"count\"], axis=1)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
holidayworkingdaytemphumiditywindspeedhour_sinhour_cosmonth_sinmonth_cosweekday_cosweekday_sin
0009.84810.00.0000001.0000000.01.0-0.222521-0.974928
1009.02800.00.2588190.9659260.01.0-0.222521-0.974928
2009.02800.00.5000000.8660250.01.0-0.222521-0.974928
\n", "
" ], "text/plain": [ " holiday workingday temp humidity windspeed hour_sin hour_cos \\\n", "0 0 0 9.84 81 0.0 0.000000 1.000000 \n", "1 0 0 9.02 80 0.0 0.258819 0.965926 \n", "2 0 0 9.02 80 0.0 0.500000 0.866025 \n", "\n", " month_sin month_cos weekday_cos weekday_sin \n", "0 0.0 1.0 -0.222521 -0.974928 \n", "1 0.0 1.0 -0.222521 -0.974928 \n", "2 0.0 1.0 -0.222521 -0.974928 " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X.head(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Notes about Cross Validation of Timeseries Data\n", "\n", "Here is the the notes from [sklearn official page](https://scikit-learn.org/stable/modules/cross_validation.html#cross-validation) :\n", "\n", "Time series data is characterised by the **correlation between observations that are near in time (autocorrelation)**.\n", "\n", "However, classical cross-validation techniques such as `KFold` and `ShuffleSplit` \n", "- assume the samples are **independent** and **identically distributed**, and \n", "- would result in unreasonable correlation between training and testing instances (yielding poor estimates of generalisation error) on time series data. \n", "\n", "Therefore, it is very important to evaluate our model for time series data on the “future” observations least like those that are used to train the model. \n", "\n", "### Time Series Split\n", "\n", "`TimeSeriesSplit` is a variation of k-fold which \n", "- returns first `k` folds as train set and \n", "- the `k+1 th` fold as test set. \n", "\n", "\n", "- We should **not shuffle** our data when making predictions with timeseries.\n", "- Unlike standard cross-validation methods, successive training sets are **supersets** of those that come before them. \n", "- Also, it adds all **surplus data** to the **first training partition**, which is always used to train the model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sklearn Pipeline \n", "Whenever possible, using Sklearn Pipeline object is always a smart practice because they \n", "- are powerfull tools to standardise our operations,\n", "- create an easy-to-understand workflow with clear order of steps,\n", "- are reproducable\n", "\n", "We will create a pipeline with `StandartScaler` and `LinearRegression` objects" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Instantiate the pipeline with the StandardScaler and LinearRegression\n", "pipeline=Pipeline(steps= [(\"scaler\", StandardScaler()),\n", " (\"linreg\", LinearRegression())])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cross validation scores\n", "\n", "Let's use the `cross_val_score` from `sklearn.metrics` to get the scores of each split " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "R^2 scores of each split: [0.27492746 0.32817026 0.19586317 0.37070607 0.41382806]\n" ] } ], "source": [ "# Split the timeseries data\n", "split = TimeSeriesSplit(n_splits=5)\n", "\n", "# Fit and score the model with cross-validation\n", "scores = cross_val_score(pipeline, X, y, cv=split)\n", "\n", "print(\"R^2 scores of each split:\", scores)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualizing CV Splits for Timeseries: \"TimeSeriesSplit\"\n", "\n", "Visualize the splits of cross validation with `TimeSeriesSplit` object" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAFNCAYAAACqr6PiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYHXWd7/F3d0JYQhIgaS/BhYjIl0VkR7yCcBXksokLwgxeWWQRRhQdQEEJMLiCouIIBEkQ3J0ZcEVk9TJwUdkEleUrKiDIFkJICARC0n3/qEpoQndyutPnnF+636/nycOpOrV8T/2eIp/8qupXHT09PUiSJKn9OttdgCRJkioGM0mSpEIYzCRJkgphMJMkSSqEwUySJKkQBjNJkqRCjG53AZLaIyK+Dry1ntwUuA+YX0+/GfgNsEtmPjVE+/sEcCDQAYwCfgV8KjMXDHA7vwSOz8y7BlnHIcBHqP7/N5rqdx6XmXOWs979wH715ImZuV9EbAcclplH9bd8Zt4ygNpOAyZl5jGNrrOMbW27uM4V3Zak1jGYSSNUZn508ec6RLx/qRCx5VDtKyLeB7wbeHNmzo+I1YD/Ak4DPjWQbWXmnitQx3bAKcC2mflkRIwCzgHOowqNjez/Fl4MaJsBrxpsPc20VJ2SVhIGM0l9iogeoAvYG3gv1a0P6wMPARcAxwAbAV/JzLPqdQ4D/qVedhZwTGbeA0ym6iVbHZifmc9FxDHAK+r1xgBnADvXy/0e+Ghmzq1D4++AN1KFuK9S90RFxD7AycAY4FmqnrTfRMTGwAxgNaoeuumZeW5dRyewBvBkZi6KiFOoAtbiHqsNgVfXy94OHJ6Zc3sdl12AbwB7AKcDEyLiW5l5aB+H8cMRsQWwKnBWZl5Yb6PPuut1No6IX9f7fwz4p8x8JCL2rn//mPq4XZyZUyPi+8CtvdrgaGAXqrD5jcx8Q0RMoAqgWwI9wOVUvZULI+J54KfAFrw8nEtqMe8xk9SInYCjqMLRq4F/At4O7Al8NiI6I2Jn4GBgp8zcCjgT+HG9/sXAU8CjEfGbiDgLeE1m3lR/fyKwENgmM7cAHga+2Gv/f8rMTTJz8faIiNcDnwf2rPd3JHBpRIwFTgB+npnb1DW+NSI6qQLJ/wPuj4jbIuIbwHbA/+21r52B/YGN65pO6euAZOaD9XfX9xPKoAqhWwO7AV+IiM2WUzfABsD+mbkxMBs4PCI6gOOAgzNzW2AH4KSImEQVkg/ptc9D6nm9fZ0qKG8ObEsVwo6vvxtTH6swlEntZ4+ZpEbcXAcRIuI+4MrM7I6Iv1L1Sq0B7EXV23RjRCxeb+2IWCcznwTeEREbAP+Lqkfnsog4NzM/SdUrtxawW73uGODxXvu/vo+adqPqVbqm1/666xp+DHw7IrYHrqbqfeuuv39/RJxQ17EzVWi8Bjig3sZ/ZuZj9W+dAXyNF0PMQJ0PkJkPR8SVVGF24TLqBrgqM2fWn+8AXpGZPXUv294RcSCwCVVP4FiqULlafU/Zs1S9nNfUv22xPYC3ZGYP8HxETAM+xovht6/jK6kNDGaSGvH8UtMv9LHMKOA7ddCi7qFaD5hd3/h/Q2beCPwNmBERO1I9APDJet1jM/Pyet01qQLfYvP62d81mbk4UBERrwYezsw76p6p3ajC0KkRsQ3wDuCJzPwZ8D3gexHxWaoetA/Xm1nYax+dwKJlHZjl6L1uJ9VxG91f3VT34fU+tj1AR92b9nuqwHk9cCHwLqCjDm0zgIOo2mlGPa93HZ31tnpPr9Jruq/jK6kNvJQpaahcAfxzREyup4+i6rmBqkftixGxTq/lNwdu67XuMRExpg50FwBfWM7+rqHqhdsYICL2BP4ArF7fd3VAZv6Q6p63ucDrqHqmzoiI3jfsbwY8QHXZEGDfiJhQ13EE8PNl1LCQlwacpR1S1/YaYNe65n7rXsZ2Xg+MB07OzJ9T9TiuShVOAS4C3gm8D/hWH+svPr4dEbEq1eXTq5axP0ltYjCTNCQy80qqG/iviog/UD3l+J768tlnqC4p3hgRd0fEn4Edqe7lov7+fqpeobuoLtMdt5z93UUVMH4YEXfU23hnZs6rP7+/nv87qp6m/87Mi4B/B34ZERkR9wBHA7tn5uLerceAXwJ3A3Oo7gfrz2+BDSLi0n6+Xy0ibqu395HM/PNy6u7PH4BfAPdExN3APlTHacP6WDxKFXL/kJkP97H+R6keGPhj/SeBzy1jf5LapKOnp2f5S0nSCDCU44hJ0mDYYyZJklQIe8wkSZIKYY+ZJElSIVamYDYamIJDfEiSpPINKresTCFnfeAvVCOQP9TmWiRJkpblVVTjDm4I/LXRlVamYLZ4bCRHqJYkSSuLyQzTYPYIwOzZz9Dd3ZoHFiZOXJNZsxwQuyS2SZlsl/LYJuWxTcrTzDbp7Oxg7bXHQp1fGrUyBbNFAN3dPS0LZov3p7LYJmWyXcpjm5THNilPC9pkQK91W5lu/pckSRrWDGaSJEmFWJkuZUqSpEFYtGghs2fPZOHCBe0upSiPP95Jd3f3Cm1j9OgxrL12F6NGDU2kMphJkjTMzZ49k9VWW4OxY9elo6Oj3eUUY/ToThYuHHww6+np4Zln5jJ79kwmTZq8/BUa4KVMSZKGuYULFzB27HhD2RDr6Ohg7NjxQ9oTaTCTJGkEMJQ1x1AfV4OZJElSIbzHTJIktdRZZ53BH/94BwsXvsBDDz3IlCkbAPC+9/0Te+31zuWuP336NDbeeBN23HHnfpc55JADueii7w9Zza1iMJMkSS9z7Flz+px/9nETVnjbxx33SQAeeeRhPvKRDw04QB1++FHLXWZlDGVgMJMkSYWYMeN87rzzTzz++KO8970HMGXKa/nmN8/l+eef4+mn5/HRj36cnXbahc997jS22mobttpqGz71qePZYIPX8ec/J+usM5HPfOaLjB8/gR133JYbbriFGTPO54knZvLgg3/nscceZe+99+Xggw9j4cKFnHHGF7j99t/T1fUKOjo6OPjgw9h6623begwMZpIkqRgLFjzPd7/7nwCcfPInOPHEqay//hRuvfVmzj77y+y00y4vWf4vf7mXk046hY022phPf/oErrzycvbb759etsy5505n3ryn2X//d/Ge9+zPFVdcxvz58/n+9y/hscce5aCDXrpOuxjMJElSMTbd9A1LPk+d+hluvPF6fv3rq7nzzj8yf/78ly2/9trrsNFGGwOwwQYbMnfu3Jcts/XW27LKKquw9trrMH78eJ55Zh433/w73vWu99DR0cG6605mm222a96PGgCfypQkScVYddVVl3z+8IeP4O677yRiYw466IP09Lz8heNjxox5yfTyluno6KCnp4fOzlF9LttuBjNJklScuXPn8OCDD3DYYUexww5v4frrr1vh1yf1tu2223PVVVfQ09PDE0/M5Pe/v7WIsd68lClJkoozfvwE9t57Xz7wgf0ZPXo0W2+9Hc8991yflzMHY99938Pf/nYvBx10ABMnTmLddSe/pLeuXTpK7MbrxxTgvlmz5tHd3Zqau7rGMXPm0y3Zlxpjm5TJdimPbVKedrbJo48+wLrrrj+gdZo5XEYJbrzxBjo7YYcddmTevHkceuj7mTHj24wfP/Df19fx7ezsYOLENQFeC9zf6LbsMZMkSS8zXAJYf6ZMeS2f/eypTJt2LgCHH/6hQYWyoWYwW8pL/4XQ978WSjT13kPbXcKQmTTt0naXIEka5tZb75V885sXsnDh0N23NhTaEswi4kDgZGAV4GuZeU476pAkSSpJy5/KjIhXAp8DdgS2BI6MiE1bXYckSVJp2jFcxq7AtZn5ZGY+A/wXsF8b6pAkSSpKOy5lrgc80mv6EWD7Rleun3BoopXnvrLhqqtr3Ap9r/awXcpjm5SnXW3y+OOdjB7t0KV9GYrj0tnZOWRt245g1gn0Hu+iA2j4zrtWDpeh9ljW4+QOAVAm26U8tkl52tkm3d3dxd3kXoLRozuH5Lh0d3e/rG17DZcxIO2Izw8Bk3tNrws83IY6JElSP5446j19/llRRx99GFdffcVL5s2fP58993w7Tz31VJ/rHHPMkdx22y3cc89dfPGLn3nZ94888jD77bfPMvd7111/4txzvw7ADTdcx/Tp0wb5C5qrHT1mVwOnRUQX8AzwXuDINtQhSZJabK+93smVV/6KXXfdfcm86667lq233pa11lprmetuvPGmnHji4J4XvP/++5g9+0kAdtxxZ3bccedBbafZWh7MMvMfEfFp4NfAGGB6Zt7U6jr603tAvZXrUoBjf0mSyve2t+3GOeeczdy5c5YM6HrFFb9k//0P5Nprr+aHP/wuzz//PC+8sICTTjqFzTffYsm6t912Cxde+E2+8Y1v8uc/37Ok92zDDTdasszf/vYXvvrVLzF//nxmz36SD3zgEN7+9t2ZPn0a8+fP5+KLZ9DV9Qp+//tbOfXU0/nTn/7I2Wd/mQULFrDWWmtxwgmf4lWvejXHHHMkm266GXfccTtPPTWbj33sBN785rc0/fi05U7AzPx+Zr4hMzfKzDPbUYMkSWq9NdZYg5122plrr70agCeemMnf//4A22+/Az/96SWceebXuPjiH3DggQfxne9c1O92PvvZUzn66I9w4YXfY731Xrlk/s9//lMOPvgwpk//Nl//+jTOOefrjBs3jsMPP4odd3wrBx982JJlX3jhBU477VP8679+gosv/gH77vteTjvt072+X8j553+Lj3zkX7nggvOG/mD0wUc0JElSS+255z5L7jO78srL2X33PRk1ahSf//yXuOmm3zB9+jQuv/wXzJ//bJ/rP/XUUzzxxBNst90OAOyxx95LvjvmmI+xYMECvvOdb3HBBef1uw2Av//9AcaNG8cmm2wGwNvetisPPfQg8+bNA+BNb3ozABts8Dqefnruiv/wBhjMJElSS2255dbMmvUEjz32KFdccTl77fVOnn32WY444mAefvgfbLHFVuy33wH09PQ9CkNHBy/5btSoF+/MOuWUE/nv//41U6a8liOP/Jdl1tHd3dcTmT10dy8CYMyYMfX+OvqtZagZzCRJUsv97/+9F9/+9oWMHz+eV77yVTz44N/p6OjgoIM+yNZbb8t11/26n+AEEyasxbrrrsuNN94AwFVX/WrJdzfffBOHH34UO+20C7/97Y0ALFq0iFGjRrFo0aKXbGf99acwZ84c7r77TgCuueYq/sf/mNzWl5n7EnNJktRye+65D/vttw8nnXQKABtu+Ho23HAjDjxwPzo7O9h++zfzhz/c3u/6U6d+hi984d+44IJz2WyzNy6Z/8EPHsHRRx/OqquO4XWvez2TJ6/HI488zCabbMaFF36T8877d9ZffwpQ9YidfvoX+MpXzuS55+YzfvwETj/9C0393cvT0aquuSEwBbivlQPMrlxPZY4MtkmZbJfy2CblaWebPProA6y77vpt2XfJhmqA2b6Ob68BZl8L3N/otryUKUmSVAiDmSRJUiEMZpIkjQAr0a1LK5WhPq4GM0mShrnRo8fwzDNzDWdDrKenh2eemcvo0WOGbJs+lSlJ0jC39tpdzJ49k3nz+n5J+EjV2dnZ75AcjRo9egxrr901RBUZzCRJGvZGjRrNpEmT211GcUp8etlLmZIkSYUwmEmSJBXCYCZJklQIg5kkSVIhvPl/KceeNafX1Jx+l9OKm3rvoQNe54km1DFUJk27tN0lSJJWcvaYSZIkFcJgJkmSVAiDmSRJUiEMZpIkSYUwmEmSJBXCYCZJklQIg5kkSVIhHMdsKWcfN2HJ5xJfbjq8DHzcL9tEkjSc2WMmSZJUCIOZJElSIQxmkiRJhTCYSZIkFcJgJkmSVAiDmSRJUiEMZpIkSYUwmEmSJBXCYCZJklQIg5kkSVIhDGaSJEmFMJhJkiQVwmAmSZJUCIOZJElSIQxmkiRJhTCYSZIkFWJ0O3YaEacC+9eTl2XmJ9pRR1+OPWtOr6k5/S6ndml/m0y999B2lzCkJk27tN0lSJJqLe8xi4hdgXcAWwFbAttExLtbXYckSVJp2tFj9ghwXGYuAIiIu4HXtKEOSZKkorQ8mGXmnYs/R8TrqS5pvqXVdUiSJJWmLfeYAUTEZsBlwAmZeW+j602cuGbzigJKuIdJaqWurnFFbUdDxzYpj21SntLapF03/78FuAT4WGb+cCDrzpo1j+7unuYUJo1AM2c+vcLb6OoaNyTb0dCxTcpjm5SnmW3S2dkxqM6klgeziHg18BPggMy8ttX7lyRJKlU7esyOB1YDvhIRi+dNy8xpbahFkiSpGO24+f9Y4NhW77dRZx83Yclnu53LU0abOO6XJKk5HPlfkiSpEAYzSZKkQhjMJEmSCmEwkyRJKoTBTJIkqRAGM0mSpEIYzCRJkgphMJMkSSqEwUySJKkQBjNJkqRCGMwkSZIKYTCTJEkqhMFMkiSpEAYzSZKkQhjMJEmSCjG63QWU5tiz5vSamtPvcmoX22Swpt57aNO2/UTTtty3SdMubfEeJak17DGTJEkqhMFMkiSpEAYzSZKkQhjMJEmSCmEwkyRJKoTBTJIkqRAGM0mSpEI4jtlSzj5uwpLPXV3jmDnz6TZWo6XZJiuieWN/2S6SNDTsMZMkSSqEwUySJKkQBjNJkqRCGMwkSZIKMeBgFhGrNKMQSZKkkW65T2VGxI7ALsCZwPXA5hFxaGb+qMm1SZIkjSiN9Jh9Cfgt8C5gFrApcFwzi5IkSRqJGglmozLzamA34CeZeT8wqqlVSZIkjUANBbOI2B7YC7gqIt4AeJ+ZJEnSEGskmH0O+D4wIzPvA34OnNzUqiRJkkag5d78n5mX8tJ3uWyYmYuaV5IkSdLI1MhTmesCM4DXAzsB346IQzLzkWYXJ0mSNJI0cinzXOAnwHzgSeB2YHozi5IkSRqJGglmUzLzAqA7M1/IzE8Cr2lyXZIkSSNOI8GsOyKWLBcR4xpcT5IkSQPQSMC6FPgeMCEiPgRcC/xnU6uSJEkagRp5KvPzEfEBqhC3G/DN+tLmsHTsWXN6Tc3pdzm1i21Spua0y9R7D23Kdttl0rRLl7+QpBGtkacyj87M84Dv9Jr3ycw8Y0V2HBFfBiZl5iErsh1JkqThot9gFhFHAWsAH4+I1Xt9tQpwFDDoYBYRbwcOBi4b7DYkSZKGm2X1mL0AbE4VzjbvNX8hK/AS84hYh+ptAp8HthjsdiRJkoabfoNZZs4AZkTEuzLzJ0O4z/OBTwOvHszKEyeuOYSl9MV7mCQ1R1fXuGG9Py2fbVKe0tpkufeYATdExMeBNYEOYBTVa5neP9CdRcThwIOZeU1EHDLQ9QFmzZpHd3fPYFaVpLaaOfPplu2rq2tcS/en5bNNytPMNuns7BhUZ1Ijwew/qEb93wy4iurJzOsHvKfKAcDkiLgdWAdYMyK+mpkfH+T2JEmSho1Ggtn6mfm6iDiX6jLkaVSvaBqwzNxt8ee6x2wXQ5kkSVKlkWD2aP3fe4E3ZOb3ImKVJtbUVmcfN2HJZ7udy2OblKl57eK4X5JGlkaC2eMRcQLwG+DfImIu1ZOaKyQzLwIuWtHtSJIkDReNvJLpQ8DzmXkDcAtwOvDJplYlSZI0AjXSY/blzDwIIDM/iaFMkiSpKRrpMdsyIjqaXokkSdII10iP2cPAnRHxW2De4pmZ+dGmVSVJkjQCNRLMflP/kSRJUhMtN5hl5r/VLzHfELgTWC0zn216ZZIkSSPMcu8xi4g3AX8FLgPWAx6MiP/Z7MIkSZJGmkZu/v8ysCswKzMfAj4AnN3UqiRJkkagRoLZGpl51+KJzPwljd2bJkmSpAFoJJi9EBFrAz0AERHNLUmSJGlkaqTn63PAdcC6EfED4B3AkU2tSpIkaQRq5KnMn0fE3cBuwCjg9My8u+mVSZIkjTDLDWYRMSMzDwP+0mvef2Xmfk2tTJIkaYTpN5hFxHnAK4GdIqKr11erABs0uzBJkqSRZlk9ZjOANwBbAJf0mr8Q+G0zi5IkSRqJ+g1mmXkLcEtEXJWZ/2hhTW117Flzek3N6Xc5tYttUibbpbep9x7a7hJ4Yoi2M2napUO0JUmNWNalzP/IzP2BX0VEz9LfZ+Ybm1qZJEnSCLOsS5ln1P89phWFSJIkjXTLupR5a/3f61pXjiRJ0sjVyMj/kiRJagGDmSRJUiH6DWYRsW9EdLSyGEmSpJFsWT1mnwHui4gTI2JSqwqSJEkaqTp6el42EsYSEfFm4AhgX+AXwDmZeVOLalvaFOC+WbPm0d3df81DqatrHDNnPt2SfakxtkmZbJfy2CblsU3K08w26ezsYOLENQFeC9zf8HrL+jIzf5OZH6w3ehNwfkTcHBEHrUCtkiRJ6kNDN/9n5tzMPAd4K3A9cGFTq5IkSRqBljXA7BIR8VbgMGBP4OfADs0sSpIkaSRa1iuZJgOHAB+sZ50PfDwzn2xBXZIkSSPOsnrM7geuBD4K/CozW3PHvSRJ0gi1rGAWmXn/S2ZErJqZzze3JEmSpJFpWTf/PxwRF0fEu3vNuyQivhURDd2bJkmSpMYtK5idDowH/l+veR8C1gZOa2JNkiRJI9KygtnewIGZ+fjiGZn5D+Ag4N39riVJkqRBWVYwW5CZ85eemZlzAe8zkyRJGmLLCmaLImLc0jPreas0ryRJkqSRaVnB7AfA9IgYu3hG/Xk6cEmzC5MkSRpplvV05deAacCjEXEnVYjbBPge1YMBkiRJGkL9BrPM7AaOjIjPAdsA3cDvMvORVhUnSZI0kix3PLLMfAB4oAW1SJIkjWgOFLuUY8+a02tqTr/LqV1skzLZLuUZWJtMvffQJtXRHpOmXdruEqRBaUswi4h9gFOBscCVmXlsO+qQJEkqybKeymyKiNiA6qGCdwFvBLaOiD1aXYckSVJp2tFj9m7gR5n5EEBEHAA814Y6JEmSitKOYLYhsCAifga8BvgFMLXRlSdOXLNZddW8V0aSVnZdXS8bH70IpdY1kpXWJu0IZqOBtwK7APOAnwEHAxc1svKsWfPo7u5pVm2SpGFg5syn213Cy3R1jSuyrpGsmW3S2dkxqM6klt9jBjwKXJ2ZM+t3cf4Y2L4NdUiSJBWlHT1mvwAujoi1gKeBPYCftKEOSZKkorQ8mGXm7yLiTOAGqpehXwV8q9V19Ofs4yYs+Wy3c3lskzLZLuUZeJs47pdUgraMY5aZFwIXtmPfkiRJpWrHPWaSJEnqg8FMkiSpEAYzSZKkQhjMJEmSCmEwkyRJKoTBTJIkqRAGM0mSpEIYzCRJkgphMJMkSSqEwUySJKkQBjNJkqRCGMwkSZIKYTCTJEkqhMFMkiSpEAYzSZKkQoxudwGlOfasOb2m5vS7nNrFNimT7VKekdEmU+89tN0lNOyJBpaZNO3SptehstljJkmSVAiDmSRJUiEMZpIkSYUwmEmSJBXCYCZJklQIg5kkSVIhDGaSJEmFcByzpZx93IQln7u6xjFz5tNtrEZLs03KZLuUZ+S0ycoz7tfIaROtCHvMJEmSCmEwkyRJKoTBTJIkqRAGM0mSpEIYzCRJkgphMJMkSSqEwUySJKkQBjNJkqRCGMwkSZIKYTCTJEkqhMFMkiSpEAYzSZKkQhjMJEmSCmEwkyRJKoTBTJIkqRAGM0mSpEKMbncBpTn2rDm9pub0u5zaxTYpk+1SHtukPC+2ydR7D21jHUNr0rRL213CsNKWYBYR/wc4qZ68PDOPb0cdkiRJJWn5pcyIWAP4OrAzsAWwU0Ts2uo6JEmSStOOe8xG1fsdC6xS/5nfhjokSZKK0vJLmZn5dERMBe4BngWuA25sdP2JE9dsVmk178uQJKlRXV3j2l3CCimt/pYHs4h4I/BBYH2qFPRd4HjgS42sP2vWPLq7e5pXoCRJatjMmU+3u4RB6+oa17T6Ozs7BtWZ1I5LmbsD12Tm45n5PHARsEsb6pAkSSpKO57KvAM4MyLGUl3K3Ae4uQ11SJIkFaUd95hdGRFbAbcCLwA3AV9sdR39Ofu4CUs+N7OLU4Njm5TJdimPbVKel7aJY3+pb20ZxywzzwDOaMe+JUmSSuUrmSRJkgphMJMkSSqEwUySJKkQBjNJkqRCGMwkSZIKYTCTJEkqhMFMkiSpEAYzSZKkQhjMJEmSCmEwkyRJKoTBTJIkqRAGM0mSpEIYzCRJkgphMJMkSSqEwUySJKkQo9tdQGmOPWtOr6k5/S6ndrFNymS7lMc2Kc/K3SZT7z203SUMqUnTLm13CX2yx0ySJKkQBjNJkqRCGMwkSZIKYTCTJEkqhMFMkiSpEAYzSZKkQhjMJEmSCuE4Zks5+7gJSz53dY1j5syn21iNlmablMl2KY9tUp6Vv03KHPdruLHHTJIkqRAGM0mSpEIYzCRJkgphMJMkSSqEwUySJKkQBjNJkqRCGMwkSZIKYTCTJEkqhMFMkiSpEAYzSZKkQhjMJEmSCmEwkyRJKoTBTJIkqRAGM0mSpEIYzCRJkgphMJMkSSrE6GZuPCLGAzcCe2fm/RGxK/AVYHXgR5l5cjP3L0mStDJpWo9ZRLwJuAHYqJ5eHbgQ2BfYBNguIvZo1v4lSZJWNs3sMTsC+DDwnXp6e+DezLwPICK+C7wPuLzB7Y0C6OzsGOIyl63V+9Py2SZlsl3KY5uUxzYpT7PapNd2Rw1kvaYFs8w8HCAiFs9aD3ik1yKPAK8awCYnA6y99tihKK9hEyeu2dL9aflskzLZLuWxTcpjm5SnBW0yGfhrows39R6zpXQCPb2mO4DuAax/M7ATVaBbNIR1SZIkDbVRVKHs5oGs1Mpg9hB1r1dtXeDhAaz/PNU9a5IkSSuDhnvKFmtlMPsdEBGxIXAfcCDVwwCSJEmiheOYZeZzwCHAJcBdwD3Af7Vq/5IkSaXr6OnpWf5SkiRJajpH/pckSSqEwUySJKkQBjNJkqRCGMwkSZIKYTCTJEkqRCvHMVtpRMSBwMnAKsDXMvOcNpc0rEXEqcD+9eRlmfmJiNgV+AqwOvCjzDy5XnZLYDowHvhv4KjMXBgRrwG+C7wCSOD9mTmvxT9l2ImILwOTMvOQgR77iFgL+B6wATAT2D8zH23LDxn011mlAAAGe0lEQVQmImIf4FRgLHBlZh7rudJeEfF/gJPqycsz83jPlfaIiPHAjcDemXn/UJ0brW4fe8yWEhGvBD4H7AhsCRwZEZu2t6rhqz5x3gFsRXW8t4mIf6YafHhfYBNgu4jYo17lu8AxmbkR1Wu9jqjnnwucm5kbA7cAU1v3K4aniHg7cHCvWQM99p8Frs/MTYALgLNbUvgwFREbANOAdwFvBLauzwvPlTaJiDWArwM7A1sAO9X/T/NcabGIeBPV24E2qqdXZ+jOjZa2j8Hs5XYFrs3MJzPzGapBcPdrc03D2SPAcZm5IDNfAO6mOrHuzcz7MnMh1Un0vohYH1g9M39br3tRPX8V4K28OGDxRcD7Wvgbhp2IWIfqHyifr6cHc+z3ovpXJsAPgD3q5TU476b6V/9D9blyAPAsnivtNIrq79GxVFdYVgFewHOlHY4APsyLr3rcnqE7N1raPgazl1uPKiws9gjwqjbVMuxl5p2LT5CIeD3VJc1u+m6D/tpmEjC3Pvl6z9fgnQ98GphdTw/m2C9Zp/5+LtDV3LKHtQ2BURHxs4i4HfgX+m8Xz5UWyMynqXpV7qF6H/T9wAI8V1ouMw/PzOt7zRrKc6Ol7WMwe7lOoPfrEDqogoKaKCI2A64CTgD+Rt9t0F/bLD0fbLNBi4jDgQcz85peswdz7DuWmu+5tGJGU/XoHwa8GXgT1T0vnittEhFvBD4IrE/1l/ciqlszPFfar9FzoLj2MZi93EPA5F7T6/Ji16iaICLeAlwDnJiZF9N/G/Q3/3FgQkSMqudPxjZbEQcA76h7ZU4H3gkczsCP/T/q5YiI0cA4YFbTqx++HgWuzsyZmTkf+DFVUPNcaZ/dgWsy8/HMfJ7q8tcueK6UYCj/Hmlp+xjMXu5q4O0R0VXf2Ple4FdtrmnYiohXAz8BDszMH9azf1d9FRvWJ8mBVE87PQA8Vwc5gA/U818ArqcKFAAHAZe37EcMM5m5W2a+ITO3BE4BfpaZhzLwY//Lepr6++vr5TU4vwB2j4i16vNiD6r7YTxX2ucOYNeIGBsRHcA+wHV4rpRgKP8eaWn7OFzGUjLzHxHxaeDXwBhgembe1OayhrPjgdWAr0TE4nnTgEOAS+rvfsmLN2S+H7igfiz6NqonoqC63+biiDgZ+Dvwz60ofoQZ6LGfClwUEXcCT9Xra5Ay83cRcSbVk2erUF36P4/q/ibPlTbIzCsjYivgVqqb/m8CvkjVm+m50kaZ+VxEHMLQnBstbZ+Onp6lL6lKkiSpHbyUKUmSVAiDmSRJUiEMZpIkSYUwmEmSJBXCYCZJklQIh8uQVJSImAL8FfhjPasTmAd8LTP/o4H1TwHuyMyfDmCfpwN/ycxvD2Cd44E3ZOYhja4jSctjMJNUovn1ALfAkpeoXxMRizLzkuWs+zbgroHsLDNPGUSNkjTkDGaSipeZD9Q9YScAl0TERsA5VK9GmQzcTjUi92HAtsCXImIRcGdfy2Xmc723HxEXAX/KzC9HxHNUg4S+o17nzMw8LyJWoRqIcjeq17c8Bsyp158AnA1sTjX46zV1ra8HbgT+V2beHhHfBl7IzMOG/ihJGg68x0zSyuIOquADcARwcWbuAGwIvBbYKzPPAW4BTsjMH/e33HL2syrwRGb+T2A/4KsRsRrVqOAbAZtShbPX9Frnq8CtmbkNsBUwCfjXzLwb+ATVaOKHAVsAx6zAMZA0zNljJmll0QM8W3/+JLBbRHyCKiytB6zZxzqNLre0xfen3UYV1MZSvTD8+5m5AFgQEd8D3lgvtzewfR2+AFZfvKHMvCAidgf+HdiifgG5JPXJYCZpZbEdLz4Q8AOq/3/9B3AZVe9VRx/rNLrc0uYDZGZP/Q7Xxev0Xndhr8+jgPfVPWRExFpUQZKIWJWqt24OsCVwbwP7lzRCeSlTUvHqe8qmAmfVs3YHTs/MH9XTb6IKR1AFplUaWG6gLgcOiojV6kubB/T67grg4xHRUQexn/HiJcsvAX+iumft3+sHGSSpT/aYSSrR6hFxe/25G3gOOCkzL6vnfQr4cUQ8Q9UTdR1VrxRUoegLETFmOcsN1Pn1un8CZvHSnq+PUt38/0eqUHg1cGZE7AW8G9g8M5+KiK8CP4iIt2bmQiRpKR09PT3trkGSJEl4KVOSJKkYBjNJkqRCGMwkSZIKYTCTJEkqhMFMkiSpEAYzSZKkQhjMJEmSCvH/AYoMKJIj/012AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Initialize the cross-validation iterator with 10 splits\n", "cv = TimeSeriesSplit(n_splits=10)\n", "\n", "fig, ax = plt.subplots(figsize=(10, 5)) \n", "\n", "# Loop over the cross validation splits\n", "# cv.split() method creates train and test arrays for each split\n", "for idx, (train, test) in enumerate(cv.split(X, y)):\n", "\n", " # Plot training and test indices\n", " indeces1 = ax.scatter(train, [idx] * len(train), c=[plt.cm.coolwarm(.1)], marker='_', lw=8)\n", " indeces2 = ax.scatter(test, [idx] * len(test), c=[plt.cm.coolwarm(.9)], marker='_', lw=8)\n", " \n", " ax.set(ylim=[10, -1], title='TimeSeriesSplit behavior', xlabel='Data index', ylabel='CV iterates')\n", " ax.legend([indeces1, indeces2], ['Training', 'Validation'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scoring Metric: \"RMSLE\"\n", "\n", "Lets try to understand the metric which we will use for the Kaggle evaluation. Here is the screen shot from the Kaggle evaluation page.\n", "\n", "![scoring](/images/scoring.png)\n", "\n", "So we need to take into account the metric `Root Mean Squared Log Error(RMSLE)` for this project. For the interpretation of RMSLE take a look to [this page](https://hrngok.github.io/posts/metrics/)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## RSMLE Calculator Function" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# Define the RMSLE function for error calculation: rmsle_calculator\n", "# Using the vectorized numpy functions instead of loops always better for computation\n", "def rmsle_calculator(predicted, actual):\n", " assert len(predicted) == len(actual)\n", " return np.sqrt(\n", " np.mean(\n", " np.power(np.log1p(predicted)-np.log1p(actual), 2)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Custom Scoring Function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Having defined our function for the RMSLE calculation, now we should define a scoring function in order to use as a scoring parameter for `model_selection.cross_val_score`\n", "\n", "- We need this parameter for model-evaluation tools which rely on a scoring strategy when using cross-validation internally (such as `model_selection.cross_val_score` and `model_selection.GridSearchCV`) \n", " \n", " \n", " We will use `make_scorer` function of Sklearn to generate a callable object (from our `rsmle_calculator` function) for scoring.
\n", "When defining a custom scorer via `sklearn.metrics.make_scorer`, the convention is that \n", "- custom functions ending in `_score` return a value to maximize and\n", "- for scorers ending in `_loss` or `_error`, a value is returned to be minimized. \n", "- We can use this functionality by setting the `greater_is_better` parameter inside `make_scorer`. \n", "- This parameter would be \n", " - `True` for scorers where higher values are better, and\n", " - `False` for scorers where lower values are better. (this will be our choise since the lower RSMLE values are better)\n", "\n", "**NOTE:**\n", "- If a loss, the output of the python function, is negated by the scorer object, conforming to the cross validation convention; that scorers return higher values for better models i.e\n", "- when `greater_is_better` is `False`, the scorer object will sign-flip the outcome of the `score_func`." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# Make a custom scorer \n", "# rmsle_error will negate the return value of rmsle_calculator,\n", "rmsle_error = make_scorer(rmsle_calculator, greater_is_better=False)\n", "\n", "# Fit and score the model with cross-validation\n", "scores = cross_val_score(pipeline, X, y, cv=split, scoring=rmsle_error)\n", "scores" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As explained the scores are negative due to the Sklearn customs. We can multiply the results by `-1` for our interpretation" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The cross validation scores: [1.06115722 0.97511162 1.02056734 0.93509628 0.92447915]\n" ] } ], "source": [ "print(\"The cross validation scores:\", scores*-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hold-out Data Prediction Scores" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "R^2 score of hold-out: 0.3889173616885714\n", "Root mean squared error of hold-out: 120.48857034926871\n", "RMSLE value for linear regression of hold-out: 0.715555147546642\n" ] } ], "source": [ "# Hold-out target data \n", "y_hold=hold_out_df[\"count\"]\n", "\n", "# Hold-out features data\n", "X_hold=hold_out_df.drop([\"datetime\", \"count\"], axis=1)\n", "\n", "# Fit the pipeline to train data\n", "pipeline.fit(X, y)\n", "\n", "# Generate predictions for hold-out data\n", "predictions_hold = pipeline.predict(X_hold)\n", "\n", "# R^2 score\n", "score = r2_score(y_hold, predictions_hold)\n", "\n", "print(\"R^2 score of hold-out:\", score)\n", "print(\"Root mean squared error of hold-out:\", np.sqrt(mean_squared_error(y_hold, predictions_hold)))\n", "print (\"RMSLE value for linear regression of hold-out: \", rmsle_calculator(y_hold , abs(predictions_hold)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Kaggle submission\n", "Lets predict the given Kaggle test set and submit the predictions to Kaggle to get our score.
\n", "First it would be better if we combine our `X` and `X_hold` and `y` and `y_hold` datasets then train our model with more data in order to let our model learn better before predicting on a new test data i.e Kaggle test dataset which we have not yet uploaded." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Pipeline(memory=None,\n", " steps=[('scaler', StandardScaler(copy=True, with_mean=True, with_std=True)), ('linreg', LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False))])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Combine X and X hold: combined_train\n", "combined_X=pd.concat([X, X_hold])\n", "\n", "# Combine the y and y hold: combined_test\n", "combined_y =pd.concat([y, y_hold])\n", "\n", "# Fit the model to the combined datasets\n", "pipeline.fit(combined_X, combined_y)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
seasonholidayworkingdayweathertempatemphumiditywindspeed
datetime
2011-01-20 00:00:00101110.6611.3655626.0027
2011-01-20 01:00:00101110.6613.635560.0000
2011-01-20 02:00:00101110.6613.635560.0000
\n", "
" ], "text/plain": [ " season holiday workingday weather temp atemp \\\n", "datetime \n", "2011-01-20 00:00:00 1 0 1 1 10.66 11.365 \n", "2011-01-20 01:00:00 1 0 1 1 10.66 13.635 \n", "2011-01-20 02:00:00 1 0 1 1 10.66 13.635 \n", "\n", " humidity windspeed \n", "datetime \n", "2011-01-20 00:00:00 56 26.0027 \n", "2011-01-20 01:00:00 56 0.0000 \n", "2011-01-20 02:00:00 56 0.0000 " ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read the Kaggle test data\n", "kaggle_test=pd.read_csv(\"bike_kaggle_test.csv\", parse_dates=[\"datetime\"], index_col=\"datetime\")\n", "kaggle_test.head(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets also apply all the steps that we followed for train data to the test data in order to conform the train and test datasets" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
holidayworkingdaytemphumiditywindspeedhour_sinhour_cosmonth_sinmonth_cosweekday_cosweekday_sin
datetime
2011-01-20 00:00:000110.665626.00270.0000001.0000000.01.0-0.9009690.433884
2011-01-20 01:00:000110.66560.00000.2588190.9659260.01.0-0.9009690.433884
2011-01-20 02:00:000110.66560.00000.5000000.8660250.01.0-0.9009690.433884
\n", "
" ], "text/plain": [ " holiday workingday temp humidity windspeed \\\n", "datetime \n", "2011-01-20 00:00:00 0 1 10.66 56 26.0027 \n", "2011-01-20 01:00:00 0 1 10.66 56 0.0000 \n", "2011-01-20 02:00:00 0 1 10.66 56 0.0000 \n", "\n", " hour_sin hour_cos month_sin month_cos weekday_cos \\\n", "datetime \n", "2011-01-20 00:00:00 0.000000 1.000000 0.0 1.0 -0.900969 \n", "2011-01-20 01:00:00 0.258819 0.965926 0.0 1.0 -0.900969 \n", "2011-01-20 02:00:00 0.500000 0.866025 0.0 1.0 -0.900969 \n", "\n", " weekday_sin \n", "datetime \n", "2011-01-20 00:00:00 0.433884 \n", "2011-01-20 01:00:00 0.433884 \n", "2011-01-20 02:00:00 0.433884 " ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kaggle_test[\"month\"]=kaggle_test.index.month\n", "kaggle_test[\"weekday\"]=kaggle_test.index.dayofweek\n", "kaggle_test[\"hour\"]=kaggle_test.index.hour\n", "\n", "# Create the sin and cos components of hour, month and weekday columns\n", "kaggle_test['hour_sin'] = np.sin(kaggle_test.hour*(2.*np.pi/24))\n", "kaggle_test['hour_cos'] = np.cos(kaggle_test.hour*(2.*np.pi/24))\n", "kaggle_test['month_sin'] = np.sin((kaggle_test.month-1)*(2.*np.pi/12))\n", "kaggle_test['month_cos'] = np.cos((kaggle_test.month-1)*(2.*np.pi/12))\n", "kaggle_test[\"weekday_cos\"]= np.cos((kaggle_test.weekday)*(2.*np.pi/7))\n", "kaggle_test['weekday_sin'] = np.sin((kaggle_test.weekday)*(2.*np.pi/7))\n", "\n", "# Drop the wrong represantation of time: month, weekday, hour\n", "# Drop the correlated features: atemp, season\n", "# Drop the weather column (not relaible data records)\n", "X_kaggle_test= kaggle_test.drop([\"season\", \"month\", \"weekday\", \"hour\", \"atemp\", \"weather\"], axis=1)\n", "X_kaggle_test.head(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our Kaggle test dataset is in the same structure with our `combined_X` dataset. Time to predict and submit!" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "final_predictions=pipeline.predict(X_kaggle_test)" ] }, { "attachments": { "resim.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAApgAAADxCAIAAADgN1AlAAActUlEQVR4nO3d7Wsb2aHH8f2P+ncMmEFggV44aUhCWEGUEFoE6V0S0/qCKMZsMJu7t7i7vYQrGu7V1ku99d2ycbVJllgsTkyrOm6eupGLBKaG2S2zL5fcF3o6Z87DnJE0ts/u9/PK0mjO0zz8NA8evfUGAAB4662TbgAAAJgcQQ4AgMcIcgAAPEaQAwDgMYIcAACPEeQAAHiMIAcAwGMEOQAAHiPIAQDwGEEOAIDHCHIAADxGkAMA4LGMQf7XX/7oLdnZjyep9uOz08z9faaO8NSDDQD4HssQ5MPwVf3ol3/NWG3eQT4sP3vLjqc8M4IcAJCFc5CPYlzIsnG0Zw04gtxkFOR+hvZgoPxsPAD4yDHIR/GSSLJJY4dT6yYEOQAgi6xBPqM9NEFu4neQD1vvZeMBwEtZgzz15PLwk9LnlNwW3xCvCsulC5+SLh2PPiVetpeyQ/tFIXH5OdkT22TDFw97icJQyB+0jaFzkLvVLQxw/xOTj2lyYuIDulso8r8SAQA/eK7XyJP3YBl30RmD/Ky6+x+Hw+jK9I+UG8DOfqy5K0ydU33H0AvHyUJwGe5KE/o9Ggq1+eaYdgpy97rP/jLxFWzyMU2plCAHgBMxi7vWtYe9bkGufU9zcDh6S0wTdU7tEb+2VfIHUiZrOqDeNGB5x9rUBMtd68OCM9Ytzz7xmCr3++m+cmiXPgAgR1kfCJN6LJgxyLUnxJMHeWJZ2iNW66l74bUpX1ImO3wx0DXN2lRTXelBnq1u0+WK7GNqbKruJARBDgDHZKonu+kup2a/Rj6WDBTdpyYrX4lHQ8X6yUp5hvPfiZC0Jm5qkKcesjvVndYTc5v0iyjlBDpBDgDHbdpHtCondacJ8uTbswty05GueqpcP9nbIFeqmX5MFQQ5AJwgtyC3nA9OTjqdR+SaOvQ9MkwmyMcxrnz5IcgB4AS5Bfk43pLxosSKbleuZo/9GrnxFjND+VmCPFmO4ROJyaaWWTM6nyCfom5dT4wfTrkdUN8QghwAjpvrqXX901g1h2lqRltv3s5077lYWqYgVyNHKiZlsq6CLHeOzzjIJ697mjFVxki7Rvj9OBsA8JH7NXLzb6YY7osWpp89mzj8dSssx5vd3pKmp0ye9v/IZxzkE9c91ZhqFtn4f9FN9xqQ5wCQt4w3u6l7c92uWnlkmHIeW0gJqUj9dd/ZnFpPeaaNdbLbk90MU2cf5BPWnTYyaV1WF5WuNOuj4QAAMzbtXesAAOAEEeQAAHiMIAcAwGMEOQAAHiPIAQDwGEEOAIDHCHIAADxGkAMA4DGCHAAAjxHkAAB4jCAHAMBjBDkAAB4jyAEA8BhBDgCAxwhyAAA8RpADAOAxghwAAI8R5AAAeIwgBwDAYwQ5AAAeI8gBAPAYQQ4AgMcIcgAAPEaQAwDgMYIcAACPOQR5a2Vubu5KoyO92WlcUd88Vp3GlbmV1pQFZO/BZHNptVbm5jRD++ZNa+Vkh9Zu6oGfnXyacrIdzLDwZ9hQc1GG9pzudTR/5mGZapFMN6yDPYq1BeKSns1CnLbPok7jimGnCCunI/JO44q8qGaYZhM7RXkyEXkMpVenYidpbMQpGniC/PsU5DMuMt/FOKvGJlo5VbFOiTr7IJ+hEwmWU7RDm5zjqXU5ymf5FWxivg+/ZZ09FdsXQX4yCPIZ+UEGefrMpz7Ij33TO0U7tMk5XyMf544c48NzIdLpEHlolFWn1biiPf8zPDM0J39pENY14VW/1HH140+1VuauNBorUrNGRZsPgpWqtX1L5K92Pl0DjD0dj6rhiLxlKyhbd/XLRO2E2DxlOSUHXpwuzDioOXk2R1h70jomjr9+S7M1RZjZ3sLE+/3VM9Fi/WBqLzqZmt1pXJlbaegbJfbAtHPVrGmO3ddvP280W1OG9vTfHjdKWj6pS9ZhjUsOsWX9VUdcfEOzgWm6Li1LQ/vdhmXYLtsCV5sttVKz+aaNqm6TNe7U9EGuflzeeBOXU1srws5GKtGw6xNWloYmOzWjYVuLhRBx3AMad0+J91or2qV1umW42a2/vOQgSyxDaR02rTqmUdLPb93+5qRX0maY2CaFVVszi7bqRJHCktdVam2AbpNXvxHogtxURWLcHLurXSaG5WY9ItfndKKk8eAmvmR1nDqWHFHNSuPUlOR3B7WFyjjJWalpprzyOTVbKldsq9Ju00LWrqGZuu8U5E7tSWzIyY7al6zTGpccYvs+RfxbGBDjsa5xR2Jf5E7DIo6jPhIFaitT91bGk3lSF207Ks2IGz6e2C3N6fYnyfHXblbJv11Gw7IWS/M77gENG79yRP59D/LkklRWqPFLt41OlhzOTkeNAdOaJ0+S61BeqQXoq04k2aA9+viXX9oaIPU3NcjVfZt1t5jWXd0yMYy7Pci168BwzmR1QkmjPx06JhWnb42pKZq9Y/9ATd9C44I0NtO0K7U029bW1GWsX0jZi3QIcqf2KGMmDWbK7G5rXHJO2/pr/EblGuSGDVGYmH1YjNtJSiu1G7PLjkBJXW25+iC37cvGH1lptVbUYkx5a94lGmNAbJx1LVaSN30PaNo9Keujl7L9+5l1rUscbrmczFEKV04EOR1DJF4bZ1ELkA6wNFWrp2fMO3NLAyYNcqFd2pNV9g7axit5BKZ03h7kckHy7mpMWejZOiYuAf0HTE2R5pPbommheUEam2nf9rXNtrU1fdy1Cyl7kQ5rh1t7lLeHr12WrNMaZ96+7e01zJBlw9C1331YdCcbsgS5btfhtr1YDnP0X9yl7THl4/3PJl4m+2xrv/lbmX40XNdi5z2gafdEkBu/GU0U5GIp4xX1WIJcX7VQ8Hi524Ncd7g5ZZCnDdfUQW7ofPYgb62Yto5BWcK34/SOdRpX5jRD4tQU085d30L7XszlaM+l2dMFudQD/dcJa5Fpm8ZMgzx1ldV15hQFeUpCa4rTfWzWQZ4+qilBrq401iAXNtWV1ptO48po+a60jAf/+Qa5di123QMad08Euebkh/qtz7agHSpSvlGZVlWHw4431nVd38fkfOYzPJYGTBrk5r2Urmnp3TWdf9J1PnOQyx2SK2itzF1pNAwbv7Fb4084LSvzdyxtKdKOwfD91NRM87ZvbLZjWx2/4hi+LxuL1C1Ufesc25MYGYdgTelMtiAf1+L6Zcuh67qOvdHMYak1pyB3G1VTukqT9IdVxo8PZmgMcvxNP9HFTdklyC2bmHE0HNdiXaW6V+bdE0HeP/zQfIMS3+8f0KYHeXJohTVM+MqevClJnKT/Zu8Q5NqqpQWv32+IlSben1GQS1XIL97o5kzprm6ZGMbd9UhHjg8xAJU9w5UriVGxdkx4T17UDk2RV8zRC2MLlZVMWqk1zbQHubbZTm019VS/kNy6n9jPmrcm48ZrSizt9u62ZNPXONu+Rtqn6Je1bgNL7bqt/e7DkjnIbZuvbhU17AjkmU07J32QGz8+HC15+SZvIk8NcssmlhwN6WuwaS2eLMj1uyfNIWGm76KnwnRB/ma4TivXu8S3G/pVR1u65grQYD2am5MP6oZfFNXaJzgi11c9rtm41MXGpW+Nyd6mB7ncCu3a5dhd8zIxjPvw7eTyMu+gxgVd6V9NS2yUiZIcOjaerN132faVQq90byZaKE4w/vuZ3B3Tpm5otq2t4iyGLUS3kBy7LzXUtjW11A/Z//3M9C+WaUvWZY1Th9i0/srlaZactPtO6bqt/W7Dki3IE6207TrSRlUXcPpBNseubmkmky0ZdE5BLi6k/v8+pQb5G9NaPFGQ23ZP8oryAwhyAACm8r04m32qEOQAgDzJJ9SMl+0wKYIcAJAz8UKl8fIqJkSQAwDgMYIcAACPEeQAAHiMIAcAwGMEOQAAHiPIAQDwGEEOAIDHCHIAADxGkAMA4DGCHAAAjxHkAAB4jCAHAMBjBDkAAB4jyAEA8BhBDgCAxwhyAAA8RpADAOAxghwAAI8R5AAAeIwgBwDAY65B/t1333377bfffPPN1wAA4NRwCvLvvvvum2+++de//hUDAIDTxCnIv/32W1IcAIBTyCnISXEAAE4npyD/+uuvT7qdAABAgyAHAMBjBDkAAB4jyAEA8BhBDgCAxwhyAAA8RpADAOAxghwAAI8R5AAAeIwgBwDAYwQ5AAAeI8gBAPAYQQ4AgMcIcgAAPEaQAwDgMYIcAACPEeQAAHiMIAcAwGPHFOT79XJQa86kxVkc3vtFGP7i3uGxVwwAwPE4bUG+Xy8H5fp+lrKbteAkviQcowkGBQDwA0GQe4AgBwCY5Bjk0d76zXOFIAjnL797Z+XSOGy7O/Wb5wpBEASFczfrO93B281akDCO5+jF1ruX58MgCOcv3Fzfi/rv7tfLiTmGcSdOEUJ+v14OytffWQiDwrl3P7j9diEI53/6v/vWSqxd7Dy8fbVUCIKgcK76wcPOaBapqEb7SOjiuDnCt5v9ejmofdZuDAfs9qNu+qAAAJBfkB8+WC6FlbXtl73eweM7lXCUQK/Wq2FpcaPdOzw8eHy3WgwX/yhcw9YefHabtVJpcaN90Ou93F6rhKXlB8IstiPyxLT9ejko/2b34OXHN8L+Hx/9LDiztpteiU60u3Y+OH/r82e9Xq+9sVgaNbzbrA37/nJ7rRKG1fVXmuYkgjyc/+mHo+GSRoUjcgCASV5Bfrh5IyjeejQ8RN27MzoiP+q9ft0bHqLGzVoQLj8Yz6fLrL98eD64sTnKtd21M8HS1vh4OWOQ15riH6MPpFSiETVrodDFw63lH//b757HcdxZrwbnP/zL8HOvGteGL61BfunO3rDkraWguNqyDgoAAHGcX5Dv18tiLorXyKPOl43a1YX5UHeyWJdZ6ull6SOzCfKUSgxd1H4k2aDxtxhrkAtl6U4kEOQAAI3jD/L9ejkoXa9vP3v9+vXr15/83CnIL/xn67Xo4HDWR+QplRi6SJADAE5WbqfW/7gYak+tP79bEUIp2lqSQ6+zXlUya3ftjPReFEkB21otziDIUyrRSJxaj3Z/t/Jf9zux7dT63p1Lwgl74XqDPch1gwIAQBzneLNb1Fod3+x2t1ocJVNrtRiU399+2esdtDeWFsLE4XRrtRiUf/O4d3h48Pplp39I/2q9GoYXlz9r93q9Z9v166W3//vpeI7OejUoLv3fy8PD3utnX3VjkXuQp1SiO+zv3+xW+6zd6/Xany1fLFTuam52u1Mtjm52ix7dKobV/+lXsSbcAWgPcu2gAAAQn8i/nw3/+yycv3Cz/tF/XP3xv3/6D2G27qPbl+fDIAgKpZ999Gzw5tGzz94dvnv13a0Xka6iIJy/8N4Xh3Gs+8e0oNa0B7m9Ev35e/Hfz4R/pDP++9m4d+H85duN9x1PrRsGBQAAnrXuhpPbAIDTiSB38Y9Pb55L/bdyAACOH0EOAIDHCHIAADxGkAMA4DGCHAAAjxHkAAB4jCAHAMBjBDkAAB4jyAEA8Fh+Qd7dqd+8MHgaae0Pe8LzTqXnl67vKT9OcvTVk43mnzVFGidopFQSHT69v9nquJY0fA5r9YNH3cQkW0+yNGv45NqgULp6+1HXXER6WdmbZa59hs0CAOQgryB/tV4NS4sb7YNe7+X2WiUsjR6Mlvw1lfGU+OirJxu3f7IQBsnnmhsnmJgriQ6f3q/X3i6k/964XNKzXq/X3lgshcNfRrH2JHOzhuPV6/WefX7rvFBL5rKyN8tc+wybBQDIRU5B/vxuJbjWeDV4FW0tBeHyg/6LB8vij3/urp0Z/W75fr0yf7lWv//bxG+UmycYmStp1grnqrc3Pn3vkluQyyXt3bkUVO4+T6skc7Pk8TrcvBGcWdudsIvZm2WsfZbNAgDkIqcg7+z+6d7jv49eNmujIJd/56u78355FA7DnwBXfvvLOMEkvZLEz40Z/fNv2/e+2Bsl0X69PApySyWZmyX1LHqxcaOYyL/O76+HQWm1FaWWld6sL2+XgtEvq1prT21Whi4CAHJxHDe7dbeWiqOjtPGePnqxubQwXygo4WzM6wmC3FSJa5BLor075fERqkNPnJs17ll358Nr84VCmGhd55N3CuHCbV2QJ8pKb9aXv1oIC9d/rwvyRO2pzcrQRQBALvIP8m6zVhKO//p7+r3Ow9WLhYurDztbajjPKshtlUwQ5NFevRIKB8UOPXFuVr9nR+3G9VLpeqPdTm+dsazszTLXPsNmAQDykXOQ98Ov1hzf7bxfLwfFhYVCaXHzRRRrwzlDkDdrgWiYMumVqEG+Xy9LZSVa0G3WSmGlLtyE7dCTBPMczVoQlBYWCpW1na62de5lZW+WufYZNgsAkI9cg1wNv/6ePhwEQxzHcbOmXETNEORHvdeig8Px+Xt7JWomRYcHUlm9I2Ga8n3ErScJ5jmatSAYZt/ggw4Xo7VlZW+WufYZNgsAkI/8glwbfnEct1aLxdXW8JUuG6Y/tZ5eSZZT67rvI06VuDers14NquvDf2uPmrUwNf6MZWVvlrH2WTYLAJCLvILcFH5xHO2unR/+o3F7Y7FUrDUTO/oZBHlqJc5Bbvo+4tQT9zm6mzdG/7C9vVYJk21T7lo3l5XaLOWudXPt9mbpBjH7oAAAppFXkCcuXosXsOVHfzXaR7p5pwzy1Eqcgzx56Vy6fJ7WkyzNkh6h9rCT+Aak3LVuKyulWcpd67babc3SD2L2QQEATI5nrWNirdViWGvyDFYAOEkEOSb15NcXrqU+shUAkC+CHAAAjxHkAAB4jCAHAMBjBDkAAB4jyAEA8BhBDgCAxwhyAAA8RpADAOCx/IK8u1O/eaH/pM7LtT+ID12XnuG5rj6O/eirJxvNP2uKNE7QSKkkOnx6f7PV0c6qKelqqRAEQeFc9YNH3cQkW08yN2s2XczeLPlJrGIfjROyNwsAkIO8gvzVejUsLW60D3q9l9trlbC0/GDw2xlRa7U0/FWNx3erxfGU+OirJxu3f7IQKr8GbpxgYq4kOnx6v157uyA+/N2lpGe9Xq+9sVgKK6OHmVl6krlZM+xi9maNllav13v2+a3z4z4aJ0zSRQBADnIK8ud3K8G1xvBnOaKtpSBcftB/8WA5LN56NDxO2107M/qdy/16Zf5yrX7/tz9PhJlxgpG5kmatcK56e+PT9y65Bblc0t6dS0Hl7vO0SrI3a4ZdzN4seWkdbt4Y/oa4ccIkXQQA5CGnIO/s/une47+PXjZroyCXfzGru/N+eRQOURQNPy2HmXGCSXolrr9+9s+/bd/7Ym+URPv18ijILZVM3CxjF5WfMTWWld4s5WdMpUqjFxs3ioP0NU6YpIsAgDwcx81u3a2l4ugobbynj15sLi3MFwpKck3/M6bplTj/jKko2rtTHh+hOvQkc7PMQZ78GVNjWemVKD9jOq60u/PhtflCIRyUYJwwTRcBALOUf5B3m7WScPzX39PvdR6uXixcXH3Y2VKTa1ZBbqtkgiCP9uqVUDgoduhJ5mbNoIvZmzWo9KjduF4qXW+026OxMU6YposAgFnKOcj74Vdrju923q+Xg+LCQqG0uPkiirXJlSHIm7VANEyZ9ErUIN+vl6WyEi3oNmulsFIXbsJ26EnCNH13Lit7s+JmLQhKCwuFytpONxbHxjhhmi4CAGYp1yBXw6+/pw8HwRDHcdysKRdRMwT5Ue+16OBwfP7eXomaSdHhgVRW70iYpnwfcetJwjR9dy4re7P6eT1I3kEJ42vk2gnTdBEAMEv5Bbk2/OI4bq0Wi6ut4StdNkx/aj29kiyn1nXfR5wqyd6sWXQxe7M669Wguj78p/qoWQsH4WucMEGzAAC5yCvITeEXx9Hu2vnhPxq3NxZLxVozsaOfQZCnVuIc5KbvI049yT6H813r5rJSK1HuWo+7mzdG/y6+vVYJRyNjnBDHsX4Qsw8KAGAaeQV54uK1eAFbfvRXo32km3fKIE+txDnIk5fOpcvnaT3J3Cz3u9ZtZaVUoty1Hice4PawE6VPiE2DmH1QAACT41nrmFhrtRjWmjyDFQBOEkGOST359YVrqY9sBQDkiyAHAMBjBDkAAB4jyAEA8BhBDgCAxwhyAAA8RpADAOAxghwAAI8R5AAAeCy/IO/u1G9e6D+p83LtD+JD16VneK6rj2M/+urJRvPPmiKNEzRSKokOn97fbHW0s2pKuloqBEFQOFf94FE3McnWE/dmWSrJ3sXszZKfxCpWb5yQvVkAgBzkFeSv1qthaXGjfdDrvdxeq4Sl5QeD386IWqul4a9qPL5bLY6nxEdfPdm4/ZOFUPk1cOMEE3Ml0eHT+/Xa2wXx4e8uJT3r9XrtjcVSWBk9zMzSk4zNslSSvYvZmzVaWr1e79nnt86PqzdOyN4sAEAucgry53crwbXG8Gc5oq2lIFx+0H/xYDks3no0PE7bXTsz+p3L/Xpl/nKtfv+3P0/ktXGCkbmSZq1wrnp749P3LrkFuVzS3p1LQeXu87RKsjbLUkn2LmZvlry0DjdvDH9D3DhhgmYBAHKRU5B3dv907/HfRy+btVGQy7+Y1d15vzwKhyiKhp+W89o4wSS9EtdfP/vn37bvfbE3SqL9enmUsZZKMjbLUsmQ8jOmxtrTm6X8jKk0rtGLjRvFQfoaJ2TuIgAgJ8dxs1t3a6k4Okob7+mjF5tLC/OFghLO0/+MaXolzj9jKor27pTHR6gOPcncLKWSIeVnTI1lpVei/IzpeFy7Ox9emy8UwkEJxgnTdREAMDv5B3m3WSsJx3/9Pf1e5+HqxcLF1YedLTWcZxXktkomCPJor14JhYNih55kbpZSSfaysjdrMK5H7cb1Uul6o90ejY1xwjRdBADMUs5B3s+lWnN8t/N+vRwUFxYKpcXNF1GsDecMQd6sBaJhyqRXogb5fr0slZVoQbdZK4WVunATtkNPElLnUCvJXlb2ZsXNWhCUFhYKlbWdrjQ2xgmTdxEAMFu5Brkul/br5SAcBEMcx3GzplxEzRDkR73XooPD8fl7eyVqJkWHB1JZvSNhmvJ9xK0nCfY5tJVkLyt7s/p5PUjeQQnja+TaCRN2EQAwc/kFuSmXWqvF4mpr+EqXDdOfWk+vJMupddNxcnpPMjQrw8F4SlnZm9VZrwbV9eE/1UfNWjgIX+OECZoFAMhFXkFuzqVod+388B+N2xuLpWKtmdjRzyDIUytxDnLLcXJ6T1znSD8YV+5aN9ee2izlrvW4u3lj9O/i22uVcDQyxglxHOsHMfugAACmkVeQJy5eixew5Ud/NdpHunmnDPLUSpyDPHnpXLp8ntYT12ZZK4njWHPXuq32lGYpd63HiQe4PexE6RNi0yBmHxQAwOR41jom1lothrUmz2AFgJNEkGNST3594Vr6o2QBALkiyAEA8BhBDgCAxwhyAAA8RpADAOAxghwAAI8R5AAAeIwgBwDAYwQ5AAAeI8gBAPAYQQ4AgMcIcgAAPEaQAwDgMYIcAACPEeQAAHiMIAcAwGMEOQAAHiPIAQDwWM5BHrVWS0HhnU86zhPiV+vVMFz41ZcT1ggAwA8JQQ4AgMc4tQ4AgMcIcgAAPEaQAwDgMYIcAACPEeQAAHiMIAcAwGMEOQAAHiPIAQDwGEEOAIDHCHIAADxGkAMA4DGCHAAAjxHkAAB4jCAHAMBjBDkAAB4jyAEA8BhBDgCAxwhyAAA8RpADAOCxnIM8aq2WgsI7n3ScJ8Sv1qthuPCrLyesEQCAHxKCHAAAj3FqHQAAjxHkAAB4jCAHAMBjBDkAAB4jyAEA8BhBDgCAxwhyAAA8RpADAOAxghwAAI8R5AAAeIwgBwDAYwQ5AAAeI8gBAPAYQQ4AgMcIcgAAPPb/AIAYWyVo6mEAAAAASUVORK5CYII=" } }, "cell_type": "markdown", "metadata": {}, "source": [ "Having predicted the test targets now we need to create a dataframe complying with Kaggle submission format like shown in the screenshot\n", "![resim.png](attachment:resim.png)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
count
datetime
2011-01-20 00:00:0029.463588
2011-01-20 01:00:0032.955071
2011-01-20 02:00:0022.045803
\n", "
" ], "text/plain": [ " count\n", "datetime \n", "2011-01-20 00:00:00 29.463588\n", "2011-01-20 01:00:00 32.955071\n", "2011-01-20 02:00:00 22.045803" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kaggle_sub=pd.DataFrame({\"datetime\":kaggle_test.index, \"count\":final_predictions}).set_index(\"datetime\")\n", "kaggle_sub[\"count\"]=kaggle_sub[\"count\"].abs()\n", "kaggle_sub.head(3)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# Save the submission dataframe\n", "kaggle_sub.to_csv(\"kaggle_sub.csv\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Kaggle Submission Score\n", "\n", "When we fit our model to the Kaggle test set and submited we got the score below. This competition is closed so there is no ranking anymore but to get a bit of idea about our model's performance the mean ranking and the two scores in the ranking that our score fall in between were below. \n", "
\n", "\n", "For the start of our project this is a pretty good score. \n", "- We have just tried Linear Regression so far and also \n", "- we did not used new features except from the time features. \n", "- Later we can work more on this and i think our ranking can get better" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our score:\n", "![kaggle1](/images/kaggle1.jpg)\n", "\n", "Mean value benchmark:\n", "\n", "![kaggle2](/images/kaggle2.jpg)\n", "\n", "The ranking that our score falls in\n", "\n", "![kaggle3](/images/kaggle3.jpg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wrap Up\n", "\n", "This was the second part of our analysis on bikeshare dataset
\n", "\n", "In this notebook we\n", "\n", "- continued to prepare our dataset for machine learning\n", "- used the trigonometric function for time data transformation in order to better represent the cyclic nature of time features\n", "- splitted our data set and keep an hould-out test set for final evaluation\n", "- used `TimeSeriesSplit` iterators for the cross validation of timeseries data. This was important in order to only using past data to evaluate our model on the “future” observations \n", "- used Sklearn Pipeline for better workflow\n", "- understood the mechanism and the functionality of the RSMLE metric\n", "- created a custom scoring function which uses RSMLE for the sklearn model evaluation tools \n", "- prepared and predicted the Kaggle testset and made a submission\n", "\n", "In the next post we will continue our analysis of bikeshare data with tree-based models and focus on model performance especially with visualization tools" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sources:
\n", "https://scikit-learn.org/stable/modules/model_evaluation.html#defining-your-scoring-strategy-from-metric-functions
\n", "http://blog.davidkaleko.com/feature-engineering-cyclical-features.html
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.6" }, "nikola": { "category": "", "date": "2018-11-16 22:08:04 UTC+02:00", "description": "", "link": "", "slug": "bikeshare part2", "tags": "timeseries data ,Kaggle, TimeSeriesSplit, cyclic time features", "title": "Cyclic Nature of Time- Kaggle Submission", "type": "text" } }, "nbformat": 4, "nbformat_minor": 2 }