Skip to content

Instantly share code, notes, and snippets.

@ledmaster

ledmaster/MultipleTimeSeriesForecasting.ipynb

Last active September 24, 2024 15:14
Show Gist options
  • Save ledmaster/031e531ea35494d2357926405532e254 to your computer and use it in GitHub Desktop.
Save ledmaster/031e531ea35494d2357926405532e254 to your computer and use it in GitHub Desktop.
Download ZIP
How To Predict Multiple Time Series With Scikit-Learn (With a Sales Forecasting Example)
Display the source blob
Display the rendered blob
Raw
MultipleTimeSeriesForecasting.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You got a lot of time series and want to predict the next step (or steps). What should you do now? Train a model for each series? Is there a way to fit a model for all the series together? Which is better?\n",
"\n",
"I have seen many data scientists think about approaching this problem by creating a single model for each product. Although this is one of the possible solutions, it's not likely to be the best. \n",
"\n",
"Here I will demonstrate how to train a single model to predict multiple time series at the same time. This technique usually creates powerful models that help teams win machine learning competitions and can be used in your project.\n",
"\n",
"<h1>Individual Models vs Big Model for Everything</h1>\n",
"\n",
"In machine learning, more data usually means better predictions. If you try to create one model for each series, you will have some trouble with series that have little to no data. When you concatenate all your series into a single dataset, to train a single model, you are using a lot more data. This may help the model perform better!\n",
"\n",
"Considering more than one series at a time, the model will be able to learn more subtle patterns that repeat across series. The product becomes a variable in our model.\n",
"\n",
"First, let's import our basic tools:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"%matplotlib inline\n",
"\n",
"from sklearn.metrics import mean_squared_log_error\n",
"from sklearn.ensemble import RandomForestRegressor\n",
"from lightgbm import LGBMRegressor"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<h1>Preparing the Data</h1>\n",
"\n",
"The data was obtained from the [UCI Repository](https://archive.ics.uci.edu/ml/datasets/Sales_Transactions_Dataset_Weekly). \n",
"\n",
"Here we will try to forecast sales for 811 products. Our dataset has records of sales for 52 weeks for each of the products.\n",
"\n",
"This means we have, originally, 811 time series with 52 data points each. Here I take only the Product Code and non-normalized weekly sales for each product.\n",
"\n",
"This is what the data looks like:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style>\n",
" .dataframe thead tr:only-child th {\n",
" text-align: right;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Product_Code</th>\n",
" <th>W0</th>\n",
" <th>W1</th>\n",
" <th>W2</th>\n",
" <th>W3</th>\n",
" <th>W4</th>\n",
" <th>W5</th>\n",
" <th>W6</th>\n",
" <th>W7</th>\n",
" <th>W8</th>\n",
" <th>...</th>\n",
" <th>W42</th>\n",
" <th>W43</th>\n",
" <th>W44</th>\n",
" <th>W45</th>\n",
" <th>W46</th>\n",
" <th>W47</th>\n",
" <th>W48</th>\n",
" <th>W49</th>\n",
" <th>W50</th>\n",
" <th>W51</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>P1</td>\n",
" <td>11</td>\n",
" <td>12</td>\n",
" <td>10</td>\n",
" <td>8</td>\n",
" <td>13</td>\n",
" <td>12</td>\n",
" <td>14</td>\n",
" <td>21</td>\n",
" <td>6</td>\n",
" <td>...</td>\n",
" <td>4</td>\n",
" <td>7</td>\n",
" <td>8</td>\n",
" <td>10</td>\n",
" <td>12</td>\n",
" <td>3</td>\n",
" <td>7</td>\n",
" <td>6</td>\n",
" <td>5</td>\n",
" <td>10</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>P2</td>\n",
" <td>7</td>\n",
" <td>6</td>\n",
" <td>3</td>\n",
" <td>2</td>\n",
" <td>7</td>\n",
" <td>1</td>\n",
" <td>6</td>\n",
" <td>3</td>\n",
" <td>3</td>\n",
" <td>...</td>\n",
" <td>2</td>\n",
" <td>4</td>\n",
" <td>5</td>\n",
" <td>1</td>\n",
" <td>1</td>\n",
" <td>4</td>\n",
" <td>5</td>\n",
" <td>1</td>\n",
" <td>6</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>P3</td>\n",
" <td>7</td>\n",
" <td>11</td>\n",
" <td>8</td>\n",
" <td>9</td>\n",
" <td>10</td>\n",
" <td>8</td>\n",
" <td>7</td>\n",
" <td>13</td>\n",
" <td>12</td>\n",
" <td>...</td>\n",
" <td>6</td>\n",
" <td>14</td>\n",
" <td>5</td>\n",
" <td>5</td>\n",
" <td>7</td>\n",
" <td>8</td>\n",
" <td>14</td>\n",
" <td>8</td>\n",
" <td>8</td>\n",
" <td>7</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>P4</td>\n",
" <td>12</td>\n",
" <td>8</td>\n",
" <td>13</td>\n",
" <td>5</td>\n",
" <td>9</td>\n",
" <td>6</td>\n",
" <td>9</td>\n",
" <td>13</td>\n",
" <td>13</td>\n",
" <td>...</td>\n",
" <td>9</td>\n",
" <td>10</td>\n",
" <td>3</td>\n",
" <td>4</td>\n",
" <td>6</td>\n",
" <td>8</td>\n",
" <td>14</td>\n",
" <td>8</td>\n",
" <td>7</td>\n",
" <td>8</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>P5</td>\n",
" <td>8</td>\n",
" <td>5</td>\n",
" <td>13</td>\n",
" <td>11</td>\n",
" <td>6</td>\n",
" <td>7</td>\n",
" <td>9</td>\n",
" <td>14</td>\n",
" <td>9</td>\n",
" <td>...</td>\n",
" <td>7</td>\n",
" <td>11</td>\n",
" <td>7</td>\n",
" <td>12</td>\n",
" <td>6</td>\n",
" <td>6</td>\n",
" <td>5</td>\n",
" <td>11</td>\n",
" <td>8</td>\n",
" <td>9</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>5 rows × 53 columns</p>\n",
"</div>"
],
"text/plain": [
" Product_Code W0 W1 W2 W3 W4 W5 W6 W7 W8 ... W42 W43 W44 W45 \\\n",
"0 P1 11 12 10 8 13 12 14 21 6 ... 4 7 8 10 \n",
"1 P2 7 6 3 2 7 1 6 3 3 ... 2 4 5 1 \n",
"2 P3 7 11 8 9 10 8 7 13 12 ... 6 14 5 5 \n",
"3 P4 12 8 13 5 9 6 9 13 13 ... 9 10 3 4 \n",
"4 P5 8 5 13 11 6 7 9 14 9 ... 7 11 7 12 \n",
"\n",
" W46 W47 W48 W49 W50 W51 \n",
"0 12 3 7 6 5 10 \n",
"1 1 4 5 1 6 0 \n",
"2 7 8 14 8 8 7 \n",
"3 6 8 14 8 7 8 \n",
"4 6 6 5 11 8 9 \n",
"\n",
"[5 rows x 53 columns]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = pd.read_csv('Sales_Transactions_Dataset_Weekly.csv')\n",
"data = data.filter(regex=r'Product|W')\n",
"data.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point, many data scientists would try to fit a model for each product. And while this can work well, we may have problems due to having only 52 data points for each model, which is really low!\n",
"\n",
"As machine learning models tend to improve with more data, why not try to concatenate all the series and train a single model with much more data?\n",
"\n",
"To do this we have to \"melt\" the data. This simply means transforming the data in a way that each line becomes: **Product Code, Week, Sales**\n",
"\n",
"To make our lives easier, I decided to remove the letters from the product and week code."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style>\n",
" .dataframe thead tr:only-child th {\n",
" text-align: right;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Product_Code</th>\n",
" <th>Week</th>\n",
" <th>Sales</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1</td>\n",
" <td>0</td>\n",
" <td>11</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>2</td>\n",
" <td>0</td>\n",
" <td>7</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>3</td>\n",
" <td>0</td>\n",
" <td>7</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>4</td>\n",
" <td>0</td>\n",
" <td>12</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>5</td>\n",
" <td>0</td>\n",
" <td>8</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Product_Code Week Sales\n",
"0 1 0 11\n",
"1 2 0 7\n",
"2 3 0 7\n",
"3 4 0 12\n",
"4 5 0 8"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"melt = data.melt(id_vars='Product_Code', var_name='Week', value_name='Sales')\n",
"\n",
"melt['Product_Code'] = melt['Product_Code'].str.extract('(\\d+)', expand=False).astype(int)\n",
"melt['Week'] = melt['Week'].str.extract('(\\d+)', expand=False).astype(int)\n",
"\n",
"melt = melt.sort_values(['Week', 'Product_Code'])\n",
"melt.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we have **42172 rows** to train our model."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Basic Feature Engineering\n",
"\n",
"We need to have variables to send to our model and get the predictions. For now, besides the product code and the week, I will create two features that usually help a lot with time series: lags and differences.\n",
"\n",
"**Last Week Sales**: this is simply the amount of sales that a product had in the previous week \n",
"**Last Week Diff**: the difference between the amount of sales in the previous week and the week before it (t-1 - t-2)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style>\n",
" .dataframe thead tr:only-child th {\n",
" text-align: right;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Product_Code</th>\n",
" <th>Week</th>\n",
" <th>Sales</th>\n",
" <th>Last_Week_Sales</th>\n",
" <th>Last_Week_Diff</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>1622</th>\n",
" <td>1</td>\n",
" <td>2</td>\n",
" <td>10</td>\n",
" <td>12.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1623</th>\n",
" <td>2</td>\n",
" <td>2</td>\n",
" <td>3</td>\n",
" <td>6.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1624</th>\n",
" <td>3</td>\n",
" <td>2</td>\n",
" <td>8</td>\n",
" <td>11.0</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1625</th>\n",
" <td>4</td>\n",
" <td>2</td>\n",
" <td>13</td>\n",
" <td>8.0</td>\n",
" <td>-4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1626</th>\n",
" <td>5</td>\n",
" <td>2</td>\n",
" <td>13</td>\n",
" <td>5.0</td>\n",
" <td>-3.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Product_Code Week Sales Last_Week_Sales Last_Week_Diff\n",
"1622 1 2 10 12.0 1.0\n",
"1623 2 2 3 6.0 -1.0\n",
"1624 3 2 8 11.0 4.0\n",
"1625 4 2 13 8.0 -4.0\n",
"1626 5 2 13 5.0 -3.0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"melt2 = melt.copy()\n",
"melt2['Last_Week_Sales'] = melt2.groupby(['Product_Code'])['Sales'].shift()\n",
"melt2['Last_Week_Diff'] = melt2.groupby(['Product_Code'])['Last_Week_Sales'].diff()\n",
"melt2 = melt2.dropna()\n",
"melt2.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Evaluating the Result\n",
"\n",
"To know if our model is good we need to have an evaluation metric. A metric I really like for sales forecasting is the Root Mean Squared Log Error.\n",
"\n",
"This is our well-known RMSE applied to the log of the target and the prediction output.\n",
"\n",
"It works as an approximation to the percentage error between our predictions and the target, which is a nice way to understand the errors our model is making."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def rmsle(ytrue, ypred):\n",
" return np.sqrt(mean_squared_log_error(ytrue, ypred))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Setting a Baseline and a Validation Split\n",
"\n",
"Let's say that in real life your model will run every Sunday night, so that early on Monday the managers can make decisions based on the predictions for the next week.\n",
"\n",
"To see if our model will work in this scenario, weeks it has not seen before, I will use a sliding window validation. This means that we are going to simulate training the model in all the weeks up to the one we want to predict, and evaluate our score in the new week.\n",
"\n",
"And, to avoid having a very good model in a small amount of weeks simply due to luck, I will use every week from 40 to 52, repeating the process for one at a time, and compute the score.\n",
"\n",
"As this is only a demonstration of the method, to keep it simple I am not separating a test set. In your projects, always keep some periods out of your validation to evaluate your model when you finish developing it. **This is very important**, and helps you be more secure about your model performing well when deployed.\n",
"\n",
"To make sure the model is worth using I like to set a baseline score that it has to beat. In this case, a reasonably strong baseline is using the last week amount of sales as a prediction for the sales this week."
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Week 40 - Error 0.51952\n",
"Week 41 - Error 0.51691\n",
"Week 42 - Error 0.51026\n",
"Week 43 - Error 0.50792\n",
"Week 44 - Error 0.53409\n",
"Week 45 - Error 0.52347\n",
"Week 46 - Error 0.50018\n",
"Week 47 - Error 0.49138\n",
"Week 48 - Error 0.50585\n",
"Week 49 - Error 0.50547\n",
"Week 50 - Error 0.52220\n",
"Week 51 - Error 0.55242\n",
"Mean Error = 0.51581\n"
]
}
],
"source": [
"mean_error = []\n",
"for week in range(40,52):\n",
" train = melt2[melt2['Week'] < week]\n",
" val = melt2[melt2['Week'] == week]\n",
" \n",
" p = val['Last_Week_Sales'].values\n",
" \n",
" error = rmsle(val['Sales'].values, p)\n",
" print('Week %d - Error %.5f' % (week, error))\n",
" mean_error.append(error)\n",
"print('Mean Error = %.5f' % np.mean(mean_error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we know that our model has to have a better error than about 0.51. This would be approximately 50% error in each prediction, which seems huge! But is it really?\n",
"\n",
"If we look at the distribution of sales in the dataset, we see that a lot of items sell very little amounts. So we should expect the error to look \"high\"."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x11ac50160>"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAEyCAYAAABQ2xz2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHV9JREFUeJzt3X+w3XV95/Hna0EtNVVA9E6a0A3ORLcoLcodZMetcyMK\nER2hO3UXhpVo6aQ60NEZdlbsdgdXywzdLbqFunTTkhWmlMiKNlnE0pT1ru2MKGCpAZESMNVINlkN\nAlGGTtz3/nG+tz3mfpMbcu/N/dx7no+ZM+d73ufz/X4/5z25zIvvj3NSVUiSJGnh/ZOFnoAkSZIG\nDGaSJEmNMJhJkiQ1wmAmSZLUCIOZJElSIwxmkiRJjTCYSZIkNcJgJkmS1AiDmSRJUiOOXegJHKmT\nTjqpVq1aNa/7+OEPf8iLX/zied3HYmRfprMn/ezLdPakn32Zzp70W6x9uf/++79XVS+fadyiDWar\nVq3ivvvum9d9TE5OMjExMa/7WIzsy3T2pJ99mc6e9LMv09mTfou1L0n+7nDGeSpTkiSpEQYzSZKk\nRhjMJEmSGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhphMJMkSWrEjMEsyclJvpjk4SQPJflAVz8x\nydYkj3bPJ3T1JLkuyfYkX0/y+qFtrevGP5pk3VD9jCTbunWuS5L5+LCSJEktO5wjZvuBK6rq54Gz\ngMuSnApcCdxdVauBu7vXAG8DVneP9cANMAhywFXAG4Azgaumwlw3Zv3Qemtn/9EkSZIWlxmDWVXt\nqqqvdcvPAA8DK4DzgZu6YTcBF3TL5wM318A9wPFJlgPnAluram9VPQlsBdZ2772kqr5cVQXcPLQt\nSZKkkfG8fiszySrgdcBXgLGq2gWD8JbkFd2wFcB3hlbb2dUOVd/ZU19w2777FO+58vPzvp8d17x9\n3vchSZLad9jBLMky4Hbgg1X19CEuA+t7o46g3jeH9QxOeTI2Nsbk5OQMs56dsePgitP2z+s+gHn/\nHHNt3759i27O882e9LMv09mTfvZlOnvSb6n35bCCWZIXMAhlt1TVZ7vy7iTLu6Nly4E9XX0ncPLQ\n6iuBJ7r6xAH1ya6+smf8NFW1AdgAMD4+XvP96/LX37KZa7c9r4OKR2THxRPzvo+5NDk5yXz3frGx\nJ/3sy3T2pJ99mc6e9FvqfTmcuzID3Ag8XFUfH3prCzB1Z+U6YPNQ/ZLu7syzgKe6U553AeckOaG7\n6P8c4K7uvWeSnNXt65KhbUmSJI2Mwzkc9Ebg3cC2JA90td8ErgFuS3Ip8G3gXd17dwLnAduBHwHv\nBaiqvUk+BtzbjftoVe3tlt8PfAo4DvhC95AkSRopMwazqvor+q8DAzi7Z3wBlx1kWxuBjT31+4DX\nzjQXSZKkpcxv/pckSWqEwUySJKkRBjNJkqRGGMwkSZIaYTCTJElqhMFMkiSpEQYzSZKkRhjMJEmS\nGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhphMJMkSWqEwUySJKkRBjNJkqRGGMwkSZIaYTCTJElq\nhMFMkiSpEQYzSZKkRhjMJEmSGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhoxYzBLsjHJniQPDtU+\nneSB7rEjyQNdfVWSZ4fe+4Ohdc5Isi3J9iTXJUlXPzHJ1iSPds8nzMcHlSRJat3hHDH7FLB2uFBV\n/7qqTq+q04Hbgc8Ovf3Y1HtV9b6h+g3AemB195ja5pXA3VW1Gri7ey1JkjRyZgxmVfUlYG/fe91R\nr38F3HqobSRZDrykqr5cVQXcDFzQvX0+cFO3fNNQXZIkaaRkkJNmGJSsAu6oqtceUH8T8PGqGh8a\n9xDwt8DTwG9V1V8mGQeuqaq3dON+CfhQVb0jyQ+q6vihbT5ZVb2nM5OsZ3DUjbGxsTM2bdr0/D7t\n87Rn71PsfnZedwHAaSteOv87mUP79u1j2bJlCz2NptiTfvZlOnvSz75MZ0/6Lda+rFmz5v6pvHQo\nx85yPxfxk0fLdgE/V1XfT3IG8KdJXgOkZ92ZE+GBK1RtADYAjI+P18TExPOf8fNw/S2buXbbbFs0\nsx0XT8z7PubS5OQk8937xcae9LMv09mTfvZlOnvSb6n35YhTR5JjgX8JnDFVq6rngOe65fuTPAa8\nCtgJrBxafSXwRLe8O8nyqtrVnfLcc6RzkiRJWsxm83UZbwG+WVU7pwpJXp7kmG75lQwu8n+8qnYB\nzyQ5q7su7RJgc7faFmBdt7xuqC5JkjRSDufrMm4Fvgy8OsnOJJd2b13I9Iv+3wR8PcnfAJ8B3ldV\nUzcOvB/4I2A78Bjwha5+DfDWJI8Cb+1eS5IkjZwZT2VW1UUHqb+np3Y7g6/P6Bt/H/Danvr3gbNn\nmockSdJS5zf/S5IkNcJgJkmS1AiDmSRJUiMMZpIkSY0wmEmSJDXCYCZJktQIg5kkSVIjDGaSJEmN\nMJhJkiQ1wmAmSZLUCIOZJElSIwxmkiRJjTCYSZIkNcJgJkmS1AiDmSRJUiMMZpIkSY0wmEmSJDXC\nYCZJktQIg5kkSVIjDGaSJEmNMJhJkiQ1wmAmSZLUCIOZJElSI2YMZkk2JtmT5MGh2keSfDfJA93j\nvKH3Ppxke5JHkpw7VF/b1bYnuXKofkqSryR5NMmnk7xwLj+gJEnSYnE4R8w+BaztqX+iqk7vHncC\nJDkVuBB4TbfOf01yTJJjgE8CbwNOBS7qxgL8Tret1cCTwKWz+UCSJEmL1YzBrKq+BOw9zO2dD2yq\nqueq6lvAduDM7rG9qh6vqr8HNgHnJwnwZuAz3fo3ARc8z88gSZK0JBw7i3UvT3IJcB9wRVU9CawA\n7hkas7OrAXzngPobgJcBP6iq/T3jp0myHlgPMDY2xuTk5CymP7Ox4+CK0/bPPHCW5vtzzLV9+/Yt\nujnPN3vSz75MZ0/62Zfp7Em/pd6XIw1mNwAfA6p7vhb4VSA9Y4v+I3N1iPG9qmoDsAFgfHy8JiYm\nntekn6/rb9nMtdtmk10Pz46LJ+Z9H3NpcnKS+e79YmNP+tmX6exJP/synT3pt9T7ckSpo6p2Ty0n\n+UPgju7lTuDkoaErgSe65b7694DjkxzbHTUbHi9JkjRSjujrMpIsH3r5y8DUHZtbgAuTvCjJKcBq\n4KvAvcDq7g7MFzK4QWBLVRXwReBXuvXXAZuPZE6SJEmL3YxHzJLcCkwAJyXZCVwFTCQ5ncFpxx3A\nrwNU1UNJbgO+AewHLquqH3fbuRy4CzgG2FhVD3W7+BCwKclvA38N3Dhnn06SJGkRmTGYVdVFPeWD\nhqequhq4uqd+J3BnT/1xBndtSpIkjTS/+V+SJKkRBjNJkqRGGMwkSZIaYTCTJElqhMFMkiSpEQYz\nSZKkRhjMJEmSGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhphMJMkSWqEwUySJKkRBjNJkqRGGMwk\nSZIaYTCTJElqhMFMkiSpEQYzSZKkRhjMJEmSGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhoxYzBL\nsjHJniQPDtX+c5JvJvl6ks8lOb6rr0rybJIHuscfDK1zRpJtSbYnuS5JuvqJSbYmebR7PmE+Pqgk\nSVLrDueI2aeAtQfUtgKvrapfAP4W+PDQe49V1end431D9RuA9cDq7jG1zSuBu6tqNXB391qSJGnk\nzBjMqupLwN4Dan9eVfu7l/cAKw+1jSTLgZdU1ZerqoCbgQu6t88HbuqWbxqqS5IkjZQMctIMg5JV\nwB1V9dqe9/4n8Omq+uNu3EMMjqI9DfxWVf1lknHgmqp6S7fOLwEfqqp3JPlBVR0/tL0nq6r3dGaS\n9QyOujE2NnbGpk2bns9nfd727H2K3c/O6y4AOG3FS+d/J3No3759LFu2bKGn0RR70s++TGdP+tmX\n6exJv8XalzVr1txfVeMzjTt2NjtJ8u+B/cAtXWkX8HNV9f0kZwB/muQ1QHpWnzkRHrhC1QZgA8D4\n+HhNTEwc0bwP1/W3bObabbNq0WHZcfHEvO9jLk1OTjLfvV9s7Ek/+zKdPelnX6azJ/2Wel+OOHUk\nWQe8Azi7Oz1JVT0HPNct35/kMeBVwE5+8nTnSuCJbnl3kuVVtas75bnnSOckSZK0mB3R12UkWQt8\nCHhnVf1oqP7yJMd0y69kcJH/41W1C3gmyVnd3ZiXAJu71bYA67rldUN1SZKkkTLjEbMktwITwElJ\ndgJXMbgL80XA1u5bL+7p7sB8E/DRJPuBHwPvq6qpGwfez+AOz+OAL3QPgGuA25JcCnwbeNecfDJJ\nkqRFZsZgVlUX9ZRvPMjY24HbD/LefcC0mweq6vvA2TPNQ5Ikaanzm/8lSZIaYTCTJElqhMFMkiSp\nEQYzSZKkRhjMJEmSGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhphMJMkSWqEwUySJKkRBjNJkqRG\nGMwkSZIaYTCTJElqhMFMkiSpEQYzSZKkRhjMJEmSGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhph\nMJMkSWqEwUySJKkRhxXMkmxMsifJg0O1E5NsTfJo93xCV0+S65JsT/L1JK8fWmddN/7RJOuG6mck\n2datc12SzOWHlCRJWgwO94jZp4C1B9SuBO6uqtXA3d1rgLcBq7vHeuAGGAQ54CrgDcCZwFVTYa4b\ns35ovQP3JUmStOQdVjCrqi8Bew8onw/c1C3fBFwwVL+5Bu4Bjk+yHDgX2FpVe6vqSWArsLZ77yVV\n9eWqKuDmoW1JkiSNjNlcYzZWVbsAuudXdPUVwHeGxu3saoeq7+ypS5IkjZRj52GbfdeH1RHUp284\nWc/glCdjY2NMTk4e4RQPz9hxcMVp++d1H8C8f465tm/fvkU35/lmT/rZl+nsST/7Mp096bfU+zKb\nYLY7yfKq2tWdjtzT1XcCJw+NWwk80dUnDqhPdvWVPeOnqaoNwAaA8fHxmpiY6Bs2Z66/ZTPXbpuP\n7PqTdlw8Me/7mEuTk5PMd+8XG3vSz75MZ0/62Zfp7Em/pd6X2ZzK3AJM3Vm5Dtg8VL+kuzvzLOCp\n7lTnXcA5SU7oLvo/B7ire++ZJGd1d2NeMrQtSZKkkXFYh4OS3MrgaNdJSXYyuLvyGuC2JJcC3wbe\n1Q2/EzgP2A78CHgvQFXtTfIx4N5u3EerauqGgvczuPPzOOAL3UOSJGmkHFYwq6qLDvLW2T1jC7js\nINvZCGzsqd8HvPZw5iJJkrRU+c3/kiRJjTCYSZIkNcJgJkmS1AiDmSRJUiMMZpIkSY0wmEmSJDXC\nYCZJktQIg5kkSVIjDGaSJEmNMJhJkiQ1wmAmSZLUCIOZJElSIwxmkiRJjTCYSZIkNcJgJkmS1AiD\nmSRJUiMMZpIkSY0wmEmSJDXCYCZJktQIg5kkSVIjDGaSJEmNMJhJkiQ1wmAmSZLUCIOZJElSI444\nmCV5dZIHhh5PJ/lgko8k+e5Q/byhdT6cZHuSR5KcO1Rf29W2J7lyth9KkiRpMTr2SFesqkeA0wGS\nHAN8F/gc8F7gE1X1u8Pjk5wKXAi8BvhZ4C+SvKp7+5PAW4GdwL1JtlTVN450bpIkSYvREQezA5wN\nPFZVf5fkYGPOBzZV1XPAt5JsB87s3tteVY8DJNnUjTWYSZKkkZKqmv1Gko3A16rq95N8BHgP8DRw\nH3BFVT2Z5PeBe6rqj7t1bgS+0G1ibVX9Wld/N/CGqrq8Zz/rgfUAY2NjZ2zatGnWcz+UPXufYvez\n87oLAE5b8dL538kc2rdvH8uWLVvoaTTFnvSzL9PZk372ZTp70m+x9mXNmjX3V9X4TONmfcQsyQuB\ndwIf7ko3AB8Dqnu+FvhVoO9QWtF/nVtvWqyqDcAGgPHx8ZqYmJjN1Gd0/S2buXbbXB1UPLgdF0/M\n+z7m0uTkJPPd+8XGnvSzL9PZk372ZTp70m+p92UuUsfbGBwt2w0w9QyQ5A+BO7qXO4GTh9ZbCTzR\nLR+sLkmSNDLm4usyLgJunXqRZPnQe78MPNgtbwEuTPKiJKcAq4GvAvcCq5Oc0h19u7AbK0mSNFJm\ndcQsyU8zuJvy14fK/ynJ6QxOR+6Yeq+qHkpyG4OL+vcDl1XVj7vtXA7cBRwDbKyqh2YzL0mSpMVo\nVsGsqn4EvOyA2rsPMf5q4Oqe+p3AnbOZiyRJ0mLnN/9LkiQ1wmAmSZLUCIOZJElSIwxmkiRJjTCY\nSZIkNcJgJkmS1AiDmSRJUiMMZpIkSY0wmEmSJDXCYCZJktQIg5kkSVIjDGaSJEmNmNWPmGturLry\n8/O+jx3XvH3e9yFJkmbHI2aSJEmNMJhJkiQ1wmAmSZLUCIOZJElSIwxmkiRJjTCYSZIkNcJgJkmS\n1AiDmSRJUiMMZpIkSY0wmEmSJDVi1sEsyY4k25I8kOS+rnZikq1JHu2eT+jqSXJdku1Jvp7k9UPb\nWdeNfzTJutnOS5IkabGZqyNma6rq9Koa715fCdxdVauBu7vXAG8DVneP9cANMAhywFXAG4Azgaum\nwpwkSdKomK9TmecDN3XLNwEXDNVvroF7gOOTLAfOBbZW1d6qehLYCqydp7lJkiQ1KVU1uw0k3wKe\nBAr4b1W1IckPqur4oTFPVtUJSe4Arqmqv+rqdwMfAiaAn6qq3+7q/wF4tqp+94B9rWdwpI2xsbEz\nNm3aNKu5z2TP3qfY/ey87uKoOW3FS+dsW/v27WPZsmVztr2lwJ70sy/T2ZN+9mU6e9JvsfZlzZo1\n9w+dWTyoY+dgX2+sqieSvALYmuSbhxibnlodov6ThaoNwAaA8fHxmpiYOILpHr7rb9nMtdvmokUL\nb8fFE3O2rcnJSea794uNPelnX6azJ/3sy3T2pN9S78usT2VW1RPd8x7gcwyuEdvdnaKke97TDd8J\nnDy0+krgiUPUJUmSRsasglmSFyf5mall4BzgQWALMHVn5Tpgc7e8BbikuzvzLOCpqtoF3AWck+SE\n7qL/c7qaJEnSyJjtebox4HNJprb1J1X1Z0nuBW5LcinwbeBd3fg7gfOA7cCPgPcCVNXeJB8D7u3G\nfbSq9s5ybpIkSYvKrIJZVT0O/GJP/fvA2T31Ai47yLY2AhtnMx9JkqTFzG/+lyRJaoTBTJIkqREG\nM0mSpEYYzCRJkhphMJMkSWqEwUySJKkRBjNJkqRGGMwkSZIaYTCTJElqhMFMkiSpEQYzSZKkRhjM\nJEmSGmEwkyRJaoTBTJIkqREGM0mSpEYYzCRJkhphMJMkSWqEwUySJKkRxy70BHR0rLry83O2rStO\n2897era345q3z9k+JEkaRR4xkyRJaoTBTJIkqREGM0mSpEYYzCRJkhpxxMEsyclJvpjk4SQPJflA\nV/9Iku8meaB7nDe0zoeTbE/ySJJzh+pru9r2JFfO7iNJkiQtTrO5K3M/cEVVfS3JzwD3J9navfeJ\nqvrd4cFJTgUuBF4D/CzwF0le1b39SeCtwE7g3iRbquobs5ibFsBc3vl5KN79KUlaqo44mFXVLmBX\nt/xMkoeBFYdY5XxgU1U9B3wryXbgzO697VX1OECSTd1Yg5kkSRopc3KNWZJVwOuAr3Sly5N8PcnG\nJCd0tRXAd4ZW29nVDlaXJEkaKamq2W0gWQb8b+DqqvpskjHge0ABHwOWV9WvJvkk8OWq+uNuvRuB\nOxmEw3Or6te6+ruBM6vqN3r2tR5YDzA2NnbGpk2bZjX3mezZ+xS7n53XXSxKY8exoH05bcVLF27n\nB7Fv3z6WLVu20NNojn2Zzp70sy/T2ZN+i7Uva9asub+qxmcaN6tv/k/yAuB24Jaq+ixAVe0eev8P\ngTu6lzuBk4dWXwk80S0frP4TqmoDsAFgfHy8JiYmZjP9GV1/y2au3eaPIxzoitP2L2hfdlw8sWD7\nPpjJyUnm+9/jYmRfprMn/ezLdPak31Lvy2zuygxwI/BwVX18qL58aNgvAw92y1uAC5O8KMkpwGrg\nq8C9wOokpyR5IYMbBLYc6bwkSZIWq9kc9ngj8G5gW5IHutpvAhclOZ3BqcwdwK8DVNVDSW5jcFH/\nfuCyqvoxQJLLgbuAY4CNVfXQLOYlSZK0KM3mrsy/AtLz1p2HWOdq4Oqe+p2HWk+SJGkU+M3/kiRJ\njTCYSZIkNcJgJkmS1AiDmSRJUiMMZpIkSY0wmEmSJDXCYCZJktQIg5kkSVIj/CFILTqrrvz8vO9j\nxzVvn/d9SJJ0II+YSZIkNcJgJkmS1AhPZUpLnKd+JWnxMJhJWhSORsAEQ6akheWpTEmSpEZ4xEyS\nhsz3kbkrTtvPxLzuQdJiZjCTFsjROjV3NBzss1xx2n7es4Q+pyTNN09lSpIkNcIjZlKP53s0yyND\nkqS54BEzSZKkRhjMJEmSGuGpTEk6yvxONkkH4xEzSZKkRnjETJKWKH+OS1p8mjlilmRtkkeSbE9y\n5ULPR5Ik6WhrIpglOQb4JPA24FTgoiSnLuysJEmSjq5WTmWeCWyvqscBkmwCzge+saCzkiQd0lyd\nLj3UdwF6ulSjpJVgtgL4ztDrncAbFmgukqSGLKWfLzNkaiapqoWeA0neBZxbVb/WvX43cGZV/cYB\n49YD67uXrwYemeepnQR8b573sRjZl+nsST/7Mp096WdfprMn/RZrX/5pVb18pkGtHDHbCZw89Hol\n8MSBg6pqA7DhaE0qyX1VNX609rdY2Jfp7Ek/+zKdPelnX6azJ/2Wel+auPgfuBdYneSUJC8ELgS2\nLPCcJEmSjqomjphV1f4klwN3AccAG6vqoQWeliRJ0lHVRDADqKo7gTsXeh4HOGqnTRcZ+zKdPeln\nX6azJ/3sy3T2pN+S7ksTF/9LkiSpnWvMJEmSRp7BTJIkqREGs4PwtzshycYke5I8OFQ7McnWJI92\nzycs5BwXQpKTk3wxycNJHkryga4+sr1J8lNJvprkb7qe/MeufkqSr3Q9+XR31/VISXJMkr9Ockf3\n2p4kO5JsS/JAkvu62sj+/UxJcnySzyT5Zvffl38+yn1J8uru38jU4+kkH1zqPTGY9fC3O//Bp4C1\nB9SuBO6uqtXA3d3rUbMfuKKqfh44C7is+/cxyr15DnhzVf0icDqwNslZwO8An+h68iRw6QLOcaF8\nAHh46LU9GVhTVacPfR/VKP/9TPk94M+q6p8Bv8jg383I9qWqHun+jZwOnAH8CPgcS7wnBrN+//Db\nnVX198DUb3eOlKr6ErD3gPL5wE3d8k3ABUd1Ug2oql1V9bVu+RkG//FcwQj3pgb2dS9f0D0KeDPw\nma4+Uj0BSLISeDvwR93rMOI9OYSR/fsBSPIS4E3AjQBV9fdV9QNGvC9DzgYeq6q/Y4n3xGDWr++3\nO1cs0FxaM1ZVu2AQUIBXLPB8FlSSVcDrgK8w4r3pTtk9AOwBtgKPAT+oqv3dkFH8O/ovwL8D/l/3\n+mXYExiE9j9Pcn/3U3sw4n8/wCuB/wv89+7U9x8leTH2ZcqFwK3d8pLuicGsX3pqfq+IfkKSZcDt\nwAer6umFns9Cq6ofd6ccVjI46vzzfcOO7qwWTpJ3AHuq6v7hcs/QkenJkDdW1esZXC5yWZI3LfSE\nGnAs8Hrghqp6HfBDltgpuiPVXYf5TuB/LPRcjgaDWb/D+u3OEbU7yXKA7nnPAs9nQSR5AYNQdktV\nfbYr2xugO/0yyeD6u+OTTH2R9aj9Hb0ReGeSHQwuh3gzgyNoo9wTAKrqie55D4Nrhs7Ev5+dwM6q\n+kr3+jMMgtqo9wUGAf5rVbW7e72ke2Iw6+dvdx7cFmBdt7wO2LyAc1kQ3XVCNwIPV9XHh94a2d4k\neXmS47vl44C3MLj27ovAr3TDRqonVfXhqlpZVasY/Dfkf1XVxYxwTwCSvDjJz0wtA+cADzLCfz8A\nVfV/gO8keXVXOhv4BiPel85F/ONpTFjiPfGb/w8iyXkM/u926rc7r17gKR11SW4FJoCTgN3AVcCf\nArcBPwd8G3hXVR14g8CSluRfAH8JbOMfrx36TQbXmY1kb5L8AoOLcI9h8D98t1XVR5O8ksHRohOB\nvwb+TVU9t3AzXRhJJoB/W1XvGPWedJ//c93LY4E/qaqrk7yMEf37mZLkdAY3irwQeBx4L93fEyPa\nlyQ/zeCa71dW1VNdbUn/WzGYSZIkNcJTmZIkSY0wmEmSJDXCYCZJktQIg5kkSVIjDGaSJEmNMJhJ\nkiQ1wmAmSZLUiP8PDSGyopdbkPgAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x11ac6a470>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"melt2['Sales'].hist(bins=20, figsize=(10,5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Creating the Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have a baseline, let's try to beat it!\n",
"\n",
"As a first model, let's train a Random Forest. Besides being a strong model with structured data (like the one we have), we usually can already get a very good result by just setting a high number of trees.\n",
"\n",
"I usually tell data scientists that a Random Forest is a very good model to use in a lazy day. You only set the number of trees to the maximum your computer can run and get a good model (hopefully).\n",
"\n",
"I will keep the Week as a feature although our model will have never seen the new week value. It may be that after week 35, for example, we have more sales for a product. This type of effect can be captured by the model.\n",
"\n",
"As decision trees can handle categorical features well even in ordinal encoding, I left them in this format.\n",
"\n",
"Anyway, be careful with these types of features."
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Week 40 - Error 0.46242\n",
"Week 41 - Error 0.45964\n",
"Week 42 - Error 0.46166\n",
"Week 43 - Error 0.46655\n",
"Week 44 - Error 0.47122\n",
"Week 45 - Error 0.43100\n",
"Week 46 - Error 0.47183\n",
"Week 47 - Error 0.44945\n",
"Week 48 - Error 0.46080\n",
"Week 49 - Error 0.44696\n",
"Week 50 - Error 0.47465\n",
"Week 51 - Error 0.48651\n",
"Mean Error = 0.46189\n"
]
}
],
"source": [
"mean_error = []\n",
"for week in range(40,52):\n",
" train = melt2[melt2['Week'] < week]\n",
" val = melt2[melt2['Week'] == week]\n",
" \n",
" xtr, xts = train.drop(['Sales'], axis=1), val.drop(['Sales'], axis=1)\n",
" ytr, yts = train['Sales'].values, val['Sales'].values\n",
" \n",
" mdl = RandomForestRegressor(n_estimators=1000, n_jobs=-1, random_state=0)\n",
" mdl.fit(xtr, ytr)\n",
" \n",
" p = mdl.predict(xts)\n",
" \n",
" error = rmsle(yts, p)\n",
" print('Week %d - Error %.5f' % (week, error))\n",
" mean_error.append(error)\n",
"print('Mean Error = %.5f' % np.mean(mean_error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nice!\n",
"\n",
"Applying the Random Forest to our initial variables we got a good reduction of the error (about 10%). This is a very good way to start.\n",
"\n",
"Now, let's try to add some features from more weeks in the past. First, let's add the lag and the difference of sales for the week before the previous."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style>\n",
" .dataframe thead tr:only-child th {\n",
" text-align: right;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Product_Code</th>\n",
" <th>Week</th>\n",
" <th>Sales</th>\n",
" <th>Last_Week_Sales</th>\n",
" <th>Last_Week_Diff</th>\n",
" <th>Last-1_Week_Sales</th>\n",
" <th>Last-1_Week_Diff</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>2433</th>\n",
" <td>1</td>\n",
" <td>3</td>\n",
" <td>8</td>\n",
" <td>10.0</td>\n",
" <td>-2.0</td>\n",
" <td>12.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2434</th>\n",
" <td>2</td>\n",
" <td>3</td>\n",
" <td>2</td>\n",
" <td>3.0</td>\n",
" <td>-3.0</td>\n",
" <td>6.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2435</th>\n",
" <td>3</td>\n",
" <td>3</td>\n",
" <td>9</td>\n",
" <td>8.0</td>\n",
" <td>-3.0</td>\n",
" <td>11.0</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2436</th>\n",
" <td>4</td>\n",
" <td>3</td>\n",
" <td>5</td>\n",
" <td>13.0</td>\n",
" <td>5.0</td>\n",
" <td>8.0</td>\n",
" <td>-4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2437</th>\n",
" <td>5</td>\n",
" <td>3</td>\n",
" <td>11</td>\n",
" <td>13.0</td>\n",
" <td>8.0</td>\n",
" <td>5.0</td>\n",
" <td>-3.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Product_Code Week Sales Last_Week_Sales Last_Week_Diff \\\n",
"2433 1 3 8 10.0 -2.0 \n",
"2434 2 3 2 3.0 -3.0 \n",
"2435 3 3 9 8.0 -3.0 \n",
"2436 4 3 5 13.0 5.0 \n",
"2437 5 3 11 13.0 8.0 \n",
"\n",
" Last-1_Week_Sales Last-1_Week_Diff \n",
"2433 12.0 1.0 \n",
"2434 6.0 -1.0 \n",
"2435 11.0 4.0 \n",
"2436 8.0 -4.0 \n",
"2437 5.0 -3.0 "
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"melt3 = melt.copy()\n",
"melt3['Last_Week_Sales'] = melt3.groupby(['Product_Code'])['Sales'].shift()\n",
"melt3['Last_Week_Diff'] = melt3.groupby(['Product_Code'])['Last_Week_Sales'].diff()\n",
"melt3['Last-1_Week_Sales'] = melt3.groupby(['Product_Code'])['Sales'].shift(2)\n",
"melt3['Last-1_Week_Diff'] = melt3.groupby(['Product_Code'])['Last-1_Week_Sales'].diff()\n",
"melt3 = melt3.dropna()\n",
"melt3.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, we train a model, now using the old variables plus the two new."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Week 40 - Error 0.44389\n",
"Week 41 - Error 0.45812\n",
"Week 42 - Error 0.43652\n",
"Week 43 - Error 0.44033\n",
"Week 44 - Error 0.46318\n",
"Week 45 - Error 0.42361\n",
"Week 46 - Error 0.45752\n",
"Week 47 - Error 0.44533\n",
"Week 48 - Error 0.44598\n",
"Week 49 - Error 0.43538\n",
"Week 50 - Error 0.47349\n",
"Week 51 - Error 0.46487\n",
"Mean Error = 0.44902\n"
]
}
],
"source": [
"\n",
"mean_error = []\n",
"for week in range(40,52):\n",
" train = melt3[melt3['Week'] < week]\n",
" val = melt3[melt3['Week'] == week]\n",
" \n",
" xtr, xts = train.drop(['Sales'], axis=1), val.drop(['Sales'], axis=1)\n",
" ytr, yts = train['Sales'].values, val['Sales'].values\n",
" \n",
" mdl = RandomForestRegressor(n_estimators=1000, n_jobs=-1, random_state=0)\n",
" mdl.fit(xtr, ytr)\n",
" \n",
" p = mdl.predict(xts)\n",
" \n",
" error = rmsle(yts, p)\n",
" print('Week %d - Error %.5f' % (week, error))\n",
" mean_error.append(error)\n",
"print('Mean Error = %.5f' % np.mean(mean_error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Great! Another reduction!\n",
"\n",
"Can you guess my next idea?\n",
"\n",
"Let's try a model with one more lag and difference, of course!"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style>\n",
" .dataframe thead tr:only-child th {\n",
" text-align: right;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: left;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Product_Code</th>\n",
" <th>Week</th>\n",
" <th>Sales</th>\n",
" <th>Last_Week_Sales</th>\n",
" <th>Last_Week_Diff</th>\n",
" <th>Last-1_Week_Sales</th>\n",
" <th>Last-1_Week_Diff</th>\n",
" <th>Last-2_Week_Sales</th>\n",
" <th>Last-2_Week_Diff</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>3244</th>\n",
" <td>1</td>\n",
" <td>4</td>\n",
" <td>13</td>\n",
" <td>8.0</td>\n",
" <td>-2.0</td>\n",
" <td>10.0</td>\n",
" <td>-2.0</td>\n",
" <td>12.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3245</th>\n",
" <td>2</td>\n",
" <td>4</td>\n",
" <td>7</td>\n",
" <td>2.0</td>\n",
" <td>-1.0</td>\n",
" <td>3.0</td>\n",
" <td>-3.0</td>\n",
" <td>6.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3246</th>\n",
" <td>3</td>\n",
" <td>4</td>\n",
" <td>10</td>\n",
" <td>9.0</td>\n",
" <td>1.0</td>\n",
" <td>8.0</td>\n",
" <td>-3.0</td>\n",
" <td>11.0</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3247</th>\n",
" <td>4</td>\n",
" <td>4</td>\n",
" <td>9</td>\n",
" <td>5.0</td>\n",
" <td>-8.0</td>\n",
" <td>13.0</td>\n",
" <td>5.0</td>\n",
" <td>8.0</td>\n",
" <td>-4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3248</th>\n",
" <td>5</td>\n",
" <td>4</td>\n",
" <td>6</td>\n",
" <td>11.0</td>\n",
" <td>-2.0</td>\n",
" <td>13.0</td>\n",
" <td>8.0</td>\n",
" <td>5.0</td>\n",
" <td>-3.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Product_Code Week Sales Last_Week_Sales Last_Week_Diff \\\n",
"3244 1 4 13 8.0 -2.0 \n",
"3245 2 4 7 2.0 -1.0 \n",
"3246 3 4 10 9.0 1.0 \n",
"3247 4 4 9 5.0 -8.0 \n",
"3248 5 4 6 11.0 -2.0 \n",
"\n",
" Last-1_Week_Sales Last-1_Week_Diff Last-2_Week_Sales Last-2_Week_Diff \n",
"3244 10.0 -2.0 12.0 1.0 \n",
"3245 3.0 -3.0 6.0 -1.0 \n",
"3246 8.0 -3.0 11.0 4.0 \n",
"3247 13.0 5.0 8.0 -4.0 \n",
"3248 13.0 8.0 5.0 -3.0 "
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"melt4 = melt.copy()\n",
"melt4['Last_Week_Sales'] = melt4.groupby(['Product_Code'])['Sales'].shift()\n",
"melt4['Last_Week_Diff'] = melt4.groupby(['Product_Code'])['Last_Week_Sales'].diff()\n",
"melt4['Last-1_Week_Sales'] = melt4.groupby(['Product_Code'])['Sales'].shift(2)\n",
"melt4['Last-1_Week_Diff'] = melt4.groupby(['Product_Code'])['Last-1_Week_Sales'].diff()\n",
"melt4['Last-2_Week_Sales'] = melt4.groupby(['Product_Code'])['Sales'].shift(3)\n",
"melt4['Last-2_Week_Diff'] = melt4.groupby(['Product_Code'])['Last-2_Week_Sales'].diff()\n",
"melt4 = melt4.dropna()\n",
"melt4.head()"
]
},
{
"cell_type": "code",
"execution_count": 89,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Week 40 - Error 0.43096\n",
"Week 41 - Error 0.43894\n",
"Week 42 - Error 0.42754\n",
"Week 43 - Error 0.41965\n",
"Week 44 - Error 0.44807\n",
"Week 45 - Error 0.40729\n",
"Week 46 - Error 0.43950\n",
"Week 47 - Error 0.43334\n",
"Week 48 - Error 0.43856\n",
"Week 49 - Error 0.43049\n",
"Week 50 - Error 0.45920\n",
"Week 51 - Error 0.45356\n",
"Mean Error = 0.43559\n"
]
}
],
"source": [
"mean_error = []\n",
"for week in range(40,52):\n",
" train = melt4[melt4['Week'] < week]\n",
" val = melt4[melt4['Week'] == week]\n",
" \n",
" xtr, xts = train.drop(['Sales'], axis=1), val.drop(['Sales'], axis=1)\n",
" ytr, yts = train['Sales'].values, val['Sales'].values\n",
" \n",
" mdl = RandomForestRegressor(n_estimators=1000, n_jobs=-1, random_state=0)\n",
" mdl.fit(xtr, ytr)\n",
" \n",
" p = mdl.predict(xts)\n",
" \n",
" error = rmsle(yts, p)\n",
" print('Week %d - Error %.5f' % (week, error))\n",
" mean_error.append(error)\n",
"print('Mean Error = %.5f' % np.mean(mean_error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As expected, the error is reduced.\n",
"\n",
"I will stop the lag and difference creation here, as I am getting repetitive with this, but feel free to try more lags, differences, and even other features to see when they stop improving."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Changing the Loss Function (and Target Distribution)\n",
"\n",
"We are evaluating our results using the log of our predictions and target, but we are training using them in their original form, so we are optimizing the model using the MSE instead of the MSLE. \n",
"\n",
"Usually you get an improvement by using the same function to evaluate and optimize your model. In this case we can do this by simply taking the log of our target variable before passing it to the fit model. \n",
"\n",
"Sometimes it's not possible to optimize the evaluation function directly. In this case, look for an approximation. We can say we are already using an approximation here (of the percentage error).\n",
"\n",
"By taking the log of the target we get the added bonus that its distribution will look closer to a normal distribution. This usually helps the model learn better."
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Week 40 - Error 0.41679\n",
"Week 41 - Error 0.42251\n",
"Week 42 - Error 0.41238\n",
"Week 43 - Error 0.40354\n",
"Week 44 - Error 0.43317\n",
"Week 45 - Error 0.39745\n",
"Week 46 - Error 0.42623\n",
"Week 47 - Error 0.41229\n",
"Week 48 - Error 0.43018\n",
"Week 49 - Error 0.43758\n",
"Week 50 - Error 0.44610\n",
"Week 51 - Error 0.47289\n",
"Mean Error = 0.42593\n"
]
}
],
"source": [
"\n",
"mean_error = []\n",
"for week in range(40,52):\n",
" train = melt4[melt4['Week'] < week]\n",
" val = melt4[melt4['Week'] == week]\n",
" \n",
" xtr, xts = train.drop(['Sales'], axis=1), val.drop(['Sales'], axis=1)\n",
" ytr, yts = train['Sales'].values, val['Sales'].values\n",
" \n",
" mdl = RandomForestRegressor(n_estimators=1000, n_jobs=-1, random_state=0)\n",
" mdl.fit(xtr, np.log1p(ytr))\n",
" \n",
" p = np.expm1(mdl.predict(xts))\n",
" \n",
" error = rmsle(yts, p)\n",
" print('Week %d - Error %.5f' % (week, error))\n",
" mean_error.append(error)\n",
"print('Mean Error = %.5f' % np.mean(mean_error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Yes! Lower error!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# A Better Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Anyone that ever stumbled over a machine learning competition solution may know about the love of competitors for Gradient Boosted Trees. I am guilty of this too!\n",
"\n",
"It's very common to have the best results with this type of model, which is an ensemble of trees like a Random Forest, but the method that it uses to create and combine the trees is different. Worth understanding.\n",
"\n",
"This is a very good scenario for us to try this model and see if we can reduce the error even more!\n",
"\n",
"There are many parameters we can tune for this model, but I will simply put 1000 trees with a small learning rate. If you want to improve the prediction, you can try tuning it and leave a comment about the improvement you got."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Week 40 - Error 0.40960\n",
"Week 41 - Error 0.39913\n",
"Week 42 - Error 0.39525\n",
"Week 43 - Error 0.39711\n",
"Week 44 - Error 0.42211\n",
"Week 45 - Error 0.37741\n",
"Week 46 - Error 0.40374\n",
"Week 47 - Error 0.39017\n",
"Week 48 - Error 0.39724\n",
"Week 49 - Error 0.40813\n",
"Week 50 - Error 0.42336\n",
"Week 51 - Error 0.45235\n",
"Mean Error = 0.40630\n"
]
}
],
"source": [
"mean_error = []\n",
"for week in range(40,52):\n",
" train = melt4[melt4['Week'] < week]\n",
" val = melt4[melt4['Week'] == week]\n",
" \n",
" xtr, xts = train.drop(['Sales'], axis=1), val.drop(['Sales'], axis=1)\n",
" ytr, yts = train['Sales'].values, val['Sales'].values\n",
" \n",
" mdl = LGBMRegressor(n_estimators=1000, learning_rate=0.01)\n",
" mdl.fit(xtr, np.log1p(ytr))\n",
" \n",
" p = np.expm1(mdl.predict(xts))\n",
" \n",
" error = rmsle(yts, p)\n",
" print('Week %d - Error %.5f' % (week, error))\n",
" mean_error.append(error)\n",
"print('Mean Error = %.5f' % np.mean(mean_error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Very good! \n",
"\n",
"This model gave us a sizeable improvement, we are almost below the 0.40 barrier. Let's take a look at our predictions for Week 52."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/mnestevao/anaconda3/lib/python3.5/site-packages/pandas/core/indexing.py:517: SettingWithCopyWarning: \n",
"A value is trying to be set on a copy of a slice from a DataFrame.\n",
"Try using .loc[row_indexer,col_indexer] = value instead\n",
"\n",
"See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
" self.obj[item] = s\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x11c796f60>"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3sAAAJcCAYAAABAE73ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X903Pdd5/vX+1tPRyJSa3UclNSyWoO84QYsRHdaIoq9\n3Rayu1Dsniu4C/cGwTmEwr2HvbSw12J7drMld2GpgHRPzt4LlMC2DlDoWqzN7aE03ZSuk21cGAdZ\n3qakHlDxyMRqPJXoTJGm434/9w+NXU3jHxrNfPT99XyckxPPR/q+9f58f8285js/zDknAAAAAEC6\nBFE3AAAAAADoPcIeAAAAAKQQYQ8AAAAAUoiwBwAAAAApRNgDAAAAgBQi7AEAAABAChH2AACxZ2av\nNTNnZrtatz9iZj+yjTqjZlY3s5f1vstomdn7zezfRd0HACA+CHsAgJ4ws8+Z2VorTC2b2X8yswEf\nf8s598+ccx/YYk/ftWm5i865AefcV3z01S0zGzGzOTO7YmZ/Z2bnzexHo+4LAJBMhD0AQC99n3Nu\nQNLrJL1e0r/+2l+wDdz/3NjjkiqSXiOpIGla0nKkHQEAEos7WwBAzznnLkn6iKRvkSQz+4SZ/YKZ\n/XdJfy/pG8zslWb2W2b2gpldMrN/d+3llWb2MjP7ldYVrr+W9L2b67fqPbjp9o+b2WfMrGZmz5nZ\n68zscUmjkv6/1tXGYzd4OeirzeyPzOwLZlY2sx/fVPPdZvYhMzveqvtpMyveaL5m9utm9itfM3bK\nzH6m9e+Z1hxrZva8mb3lJqvu9ZLe75z7knPuqnPuL5xzH9lU8z+b2eXWVb/TZvbNN9sGZvZWM5s3\ns1Uz+6SZjW/62Vb7AQAkGGEPANBzZrZP0vdI+otNwz8s6e2SBiX9jaQPSLoqaUzSt0m6X9K1APfj\nkt7aGi9K+v5b/K0fkPRubVwFe4WkI5KqzrkflnRRrauNzrnZGyz+QUlLkl7d+hu/+DXB54ik35e0\nW9IfSfqPN2nj9yT9czOzVk9Drfn8vpndI+mnJL3eOTco6Z9I+txN6pyR9P+Y2Q+a2egNfv4RSQck\nfb2kZyX97o2KmNnrJP22pJ/QxhXC35D0R2aW77AfAECCEfYAAL100sxWJT0t6b9J+sVNP3u/c+7T\nzrmrkl4l6Z9JekfrKtbnJb1X0g+2fvd/kfQfnHMV59wXJP37W/zNByXNOuf+3G0oO+f+5naNtgLp\nd0qacc6tO+fmJT2mjVB6zdPOuT9uvcfvcUnfepNyT0lykg61bn+/pGecc38r6SuS8pLuNbOcc+5z\nzrm/ukmdH2jV+jeSFltX5l5/7YfOud92ztWccw1tBNxvNbNX3qDOj0v6Defcp5xzX2m9v7Eh6b4O\n+wEAJBhhDwDQS29zzu12zr3GOfd/OOfWNv2ssunfr5GUk/RC62WGq9q4+vT1rZ+/+mt+/1bhbZ+k\n7YSVV0v6gnOu9jV/Z++m25c3/fvvJfVdewnoZs45p40rgD/UGvpf1brq5pwrS3qHNsLZ583s983s\n1TdqyDm34pz7OefcN0saljSvjQBtrZe2/pKZ/ZWZfVFfvRq35walXiPpZ6+t29b63Sfp1Z30AwBI\nNsIeAGCnuE3/rmjjStOeVjjc7Zx7RSvkSNIL2ggn19zoJY2ba33jFv7m1/pbSa8ys8Gv+TuXbrHM\nrXxQ0veb2WskfbukuetNOPd7zrnv1EYIc5Lec7tizrkrkn5FG6H0VdoIkEclfZekV0p6betX7QaL\nVyT9wqZ1u9s593XOuQ9utx8AQPIQ9gAAO84594KkJyT9qpm9wswCM/tGM/tHrV/5kKT/s/VVBEOS\nfu4W5R6T9C/N7B+2roCNtQKXtPFJlt9wkx4qkj4p6d+bWV/rA0x+TDd5H9wW5vQXkl5s9fNR59yq\nJJnZPWb2ZjPLS1qXtKaNl1K+hJm9x8y+xcx2tULo/y6p7JyrauO9jg1JVUlfp/aXyH6t35T0k2b2\n7a11coeZfa+ZDXbSDwAg2Qh7AICoTEt6uaTnJK1IOiHp7tbPflPSRyWd08YHkfzhzYo45/6zpF/Q\nxoek1CSd1MaVMGnjvX7/uvVSxn95g8V/SBtXyP5W0n+R9G+dcx/rYk4f1MaVt9/bNJaX9EuSrmjj\nZaFfL+ldN1n+61p9rEr6a21ceTvS+tlxbbzM9JI21tmZmzXhnCtp4317/1Eb67Ys6Ue30Q8AIMFs\n420GAAAAAIA04coeAAAAAKSQ97DX+vSwvzCzD7du7zezT5nZBTP7AzN7ue8eAAAAACBrduLK3k9L\n+sym2++R9F7n3AFtvI/gx3agBwAAAADIFK9hz8xGJH2vNj6ZTGZmkt6sjTfhS9IHJL3NZw8AAAAA\nkEUv+WLYHvsPko5p4+OiJakgadU5d7V1e0ntX157nZm9XdLbJemOO+74h9/0Td/kuVUAAAAAiKez\nZ89ecc7d2cky3sKemb1V0uedc2fN7E3Xhm/wqzf8OFDn3PskvU+SisWiK5VKXvoEAAAAgLgzs7/p\ndBmfV/beKOmImX2PpD5Jr9DGlb7dZrardXVvRBvfbQQAAAAA6CFv79lzzv0r59yIc+61kn5Q0sed\nc/+bpD+V9P2tX/sRSad89QAAAAAAWRXF9+zNSPoZMytr4z18vxVBDwAAAACQar4/oEWS5Jz7hKRP\ntP7915LesBN/FwAAAACyKoorewAAAAAAzwh7AAAAAJBChD0AAAAASCHCHgAAAACkEGEPAAAAAFKI\nsAcAAACvqvWGzlVWVa03om4FyJQd+eoFAAAAZNOp+UuamVtQLgjUDEPNTo3ryMTeqNsCMoErewAA\nAPCiWm9oZm5B681QtcZVrTdDHZtb4AofsEMIewAAAPBiaWVNuaD94WYuCLS0shZRR0C2EPYAAADg\nxchQv5ph2DbWDEONDPVH1BGQLYQ9AAAAeFEYyGt2alx9uUCD+V3qywWanRpXYSAfdWtAJvABLQAA\nAPDmyMRevXFsj5ZW1jQy1E/QA3YQYQ8AAABeFQbyhDwgAryMEwAAAABSiLAHAAAAAClE2AMAAACA\nFCLsAQAAAEAKEfYAAAAAIIUIewAAYMdU6w2dq6yqWm9E3QoApB5fvQAAAHbEqflLmplbUC4I1AxD\nzU6N68jE3qjbAoDU4soeAADwrlpvaGZuQevNULXGVa03Qx2bW+AKHwB4RNgDAADeLa2sKRe0P+zI\nBYGWVtYi6ggA0o+wBwAAvBsZ6lczDNvGmmGokaH+iDoCgPQj7AEAAO8KA3nNTo2rLxdoML9LfblA\ns1PjKgzko24NAFKLD2gBAAA74sjEXr1xbI+WVtY0MtRP0AMAzwh7AABgxxQG8oQ8ANghvIwTAAAA\nAFKIsAcAAAAAKUTYAwAAAIAUIuwBAAAAQAoR9gAAAAAghQh7AAAAAJBChD0AAJBo1XpD5yqrqtYb\nUbcCALHC9+wBAIDEOjV/STNzC8oFgZphqNmpcR2Z2Bt1WwAQC1zZAwAAiVStNzQzt6D1Zqha46rW\nm6GOzS1whQ8AWgh7AAAgkZZW1pQL2h/K5IJASytrEXUEAPFC2AMAAIk0MtSvZhi2jTXDUCND/RF1\nBADxQtgDAACJVBjIa3ZqXH25QIP5XerLBZqdGldhIB91awAQC3xACwAASKwjE3v1xrE9WlpZ08hQ\nP0EPADYh7AEAgEQrDOQJeQBwA7yMEwAAAABSiLAHAAAAAClE2AMAAACAFCLsAQAAAEAKEfYAAAAA\nIIUIewAAAACQQoQ9AAAAAEghwh4AAAAApBBhDwAAAABSiLAHAAAAAClE2AMAAACAFCLsAQAAAEAK\nEfYAAAAAIIUIewAAAACQQoQ9AAAAAEghwh4AAAAApBBhDwAAAABSiLAHAAAAAClE2AMAAACAFPIW\n9sysz8z+zMzOmdmnzeznW+PvN7NFM5tv/TfhqwcAAIDtqNYbOldZVbXeiLqVHZfluQNps8tj7Yak\nNzvn6maWk/S0mX2k9bP/yzl3wuPfBgAA2JZT85c0M7egXBCoGYaanRrXkYm9Ube1I7I8dyCNvF3Z\ncxvqrZu51n/O198DAADoVrXe0MzcgtaboWqNq1pvhjo2t5CJq1xZnjuQVl7fs2dmLzOzeUmfl/Qx\n59ynWj/6BTNbMLP3mln+Jsu+3cxKZlZ68cUXfbYJAAAgSVpaWVMuaH94lAsCLa2sRdTRzsny3IG0\n8hr2nHNfcc5NSBqR9AYz+xZJ/0rSN0l6vaRXSZq5ybLvc84VnXPFO++802ebAAAAkqSRoX41w7Bt\nrBmGGhnqj6ijnZPluQNptSOfxumcW5X0CUn/1Dn3Quslng1J/0nSG3aiBwAAgNspDOQ1OzWuvlyg\nwfwu9eUCzU6NqzBwwxcipUqW5w6klbcPaDGzOyU1nXOrZtYv6bskvcfM7nbOvWBmJultkv6Hrx4A\nAAA6dWRir944tkdLK2saGerPVNjJ8tyBNPL5aZx3S/qAmb1MG1cQP+Sc+7CZfbwVBE3SvKSf9NgD\nAABAxwoD+cwGnSzPHUgbb2HPObcg6dtuMP5mX38TAAAAALBhR96zBwAAAADYWYQ9AAAAAEghwh4A\nAAAApBBhDwAAAABSiLAHAAAAAClE2AMAADumWm/oXGVV1Xoj6lYAIPV8fs8eAADAdafmL2lmbkG5\nIFAzDDU7Na4jE3ujbgsAUosrewAAwLtqvaGZuQWtN0PVGle13gx1bG6BK3wA4BFhDwAAeLe0sqZc\n0P6wIxcEWlpZi6gjAEg/wh4AAPBuZKhfzTBsG2uGoUaG+iPqCADSj7AHAAC8KwzkNTs1rr5coMH8\nLvXlAs1OjaswkI+6NQBILT6gBQAA7IgjE3v1xrE9WlpZ08hQP0EPADwj7AEAgB1TGMgT8gBgh/Ay\nTgAAAABIIcIeAAAAAKQQYQ8AAAAAUoiwBwAAAAApRNgDAAAAgBQi7AEAACRUtd7QucqqqvVG1K3c\nko8+kzJ3IEp89QIAAEACnZq/pJm5BeWCQM0w1OzUuI5M7I26rZfw0WdS5g5EjSt7AAAACVOtNzQz\nt6D1Zqha46rWm6GOzS3E7iqXjz6TMncgDgh7AAAACbO0sqZc0P4wLhcEWlpZi6ijG/PRZ1LmDsQB\nYQ8AACBhRob61QzDtrFmGGpkqD+ijm7MR59JmTsQB4Q9AACAhCkM5DU7Na6+XKDB/C715QLNTo2r\nMJCPurU2PvpMytyBODDnXNQ93FaxWHSlUinqNgAAAGKlWm9oaWVNI0P9sQ47PvpMytyBXjGzs865\nYifL8GmcAAAACVUYyCci6PjoMylzB6LEyzgBAAAAIIUIewAAAACQQoQ9AAAAAEghwh4AAAAApBBh\nDwAAAABSiLAHAAAAAClE2AMAAEDilJdrOlGqqLxci7oVILb4nj0AAAAkykMnz+v4mYvXb09Pjurh\nowcj7AiIJ67sAQAAIDHKy7W2oCdJx5+5yBU+4AYIewAAAEiM+cpqR+NAlhH2AAAAkBgT+3Z3NA5k\nGWEPAAAAiTE2PKjpydG2senJUY0ND0bUERBffEALAAAAEuXhowc1fd9rNV9Z1cS+3QQ94CYIewAA\nAEicseFBQh5wG7yMEwAAAABSiLAHAAAAAClE2AMAAACAFCLsAQAAAEAKEfYAAAAAIIUIewAAAACQ\nQoQ9AECqVOsNnausqlpvRN3KTZWXazpRqqi8XIu6FSCxfBxHSTh/AJ3ge/YAAKlxav6SZuYWlAsC\nNcNQs1PjOjKxN+q22jx08ryOn7l4/fb05KgePnowwo6A5PFxHCXh/AF0iit7AIBUqNYbmplb0Hoz\nVK1xVevNUMfmFmL1DH15udb2AFWSjj9zkSt8QAd8HEdJOH8A20HYAwCkwtLKmnJB+91aLgi0tLIW\nUUcvNV9Z7WgcwEv5OI6ScP4AtoOwBwBIhZGhfjXDsG2sGYYaGeqPqKOXmti3u6NxAC/l4zhKwvkD\n2A7CHgAgFQoDec1OjasvF2gwv0t9uUCzU+MqDOSjbu26seFBTU+Oto1NT45qbHgwoo6A5PFxHCXh\n/AFshznnou7htorFoiuVSlG3AQBIgGq9oaWVNY0M9cf2gVp5uab5yqom9u0m6AHb5OM4SsL5A9ll\nZmedc8WOliHsAQAAAEC8bSfs8TJOAAAAAEghwh4AAAAApBBhDwAAAABSiLAHAAAAAClE2AMAAACA\nFPIW9sysz8z+zMzOmdmnzeznW+P7zexTZnbBzP7AzF7uqwcAQPZU6w2dq6yqWm9E3cpNJaFHX7I8\ndx/KyzWdKFVUXq7FuiaAaOzyWLsh6c3OubqZ5SQ9bWYfkfQzkt7rnPt9M/t1ST8m6dc89gEAyIhT\n85c0M7egXBCoGYaanRrXkYm9UbfVJgk9+pLlufvw0MnzOn7m4vXb05OjevjowdjVBBAdb1f23IZ6\n62au9Z+T9GZJJ1rjH5D0Nl89AACyo1pvaGZuQevNULXGVa03Qx2bW4jVFaQk9OhLlufuQ3m51hbK\nJOn4Mxe7uhrnoyaAaHl9z56ZvczM5iV9XtLHJP2VpFXn3NXWryxJuuFTemb2djMrmVnpxRdf9Nkm\nACAFllbWlAva79ZyQaCllbWIOnqpJPToS5bn7sN8ZbWj8ahqAoiW17DnnPuKc25C0oikN0j6n270\nazdZ9n3OuaJzrnjnnXf6bBMAkAIjQ/1qhmHbWDMMNTLUH1FHL5WEHn3J8tx9mNi3u6PxqGoCiNaO\nfBqnc25V0ick3Sdpt5lde6/giKS/3YkeAADpVhjIa3ZqXH25QIP5XerLBZqdGldhIB91a9cloUdf\nsjx3H8aGBzU9Odo2Nj05qrHhwVjVBBAtc+6GF9a6L2x2p6Smc27VzPolPSHpPZJ+RNLcpg9oWXDO\n/b+3qlUsFl2pVPLSJwAgXar1hpZW1jQy1B/bIJGEHn3J8tx9KC/XNF9Z1cS+3T0LZT5qAuiemZ11\nzhU7WsZj2BvXxgewvEwbVxA/5Jx72My+QdLvS3qVpL+Q9IBz7pbvzibsAQAAAMiy7YQ9b1+94Jxb\nkPRtNxj/a228fw8AAAAA4MmOvGcPAAAAALCzCHsAAAAAkEKEPQAAAABIIcIeAAAAAKQQYQ8AAAAA\nUoiwBwApVK03dK6yqmr9lt9sgy1gXQLdKy1W9cgTz6u0WI26lVvieEfaePvqBQBANE7NX9LM3IJy\nQaBmGGp2alxHJvZG3VYisS6B7j3w2Bk9Xd4IeY9+vKxDYwU9/uB9EXf1UhzvSCOu7AFAilTrDc3M\nLWi9GarWuKr1Zqhjcws8S70NrEuge6XF6vWgd81T5WrsrvBxvCOtCHsAkCJLK2vKBe2n9lwQaGll\nLaKOkot1CXTv9IUrHY1HheMdaUXYA4AUGRnqVzMM28aaYaiRof6IOkou1iXQvcMH9nQ0HhWOd6QV\nYQ8AUqQwkNfs1Lj6coEG87vUlws0OzWuwkA+6tYSh3UJdK+4v6BDY4W2sUNjBRX3F26yRDQ43pFW\n5pyLuofbKhaLrlQqRd0GACRGtd7Q0sqaRob6ebDSJdYl0L3SYlWnL1zR4QN7Yhf0NuN4R5yZ2Vnn\nXLGjZQh7AAAAABBv2wl7vIwTAAAAAFKIsAcAAAAAKUTYAwAAAIAUIuwBAAAAQAoR9gAAAAAghQh7\nAGKhWm/oXGVV1Xoj6laQcEnYl5LQo5ScPpMiy+vTx9xPPlvRgx/4c518ttKzmlneRkinXVE3AACn\n5i9pZm5BuSBQMww1OzWuIxN7o24LCZSEfSkJPUrJ6TMpsrw+fcz9vl/8mC5/8cuSpP/6mc/rPX/y\nl3rmXd8duz6BqHFlD0CkqvWGZuYWtN4MVWtc1Xoz1LG5BZ5VRceSsC8loUcpOX0mRZbXp4+5n3y2\ncj3oXfPCF7/c1RW+LG8jpBthD0CkllbWlAvaT0W5INDSylpEHSGpkrAvJaFHKTl9JkWW16ePuX/4\n/OWOxrciy9sI6UbYAxCpkaF+NcOwbawZhhoZ6o+oIyRVEvalJPQoJafPpMjy+vQx97cevKuj8a3I\n8jZCuhH2AESqMJDX7NS4+nKBBvO71JcLNDs1rsJAPurWkDBJ2JeS0KOUnD6TIsvr08fc3/a6fbr7\nFS9vG7v7FS/X2163L1Z9AnFgzrmoe7itYrHoSqVS1G0A8Khab2hpZU0jQ/3cuaIrSdiXktCjlJw+\nkyLL69PH3E8+W9GHz1/WWw/e1VXQ2yzL2wjxZ2ZnnXPFjpYh7AEAAABAvG0n7PEyTgAAAABIIcIe\nAAAAAKQQYQ8AAAAAUoiwBwAAAAApRNgDAAAAgBQi7AEAAABAChH2AKRWtd7QucqqqvVG1K3suCzP\nHciSLB/r5eWaTpQqKi/XelYzy+sT6bQr6gYAwIdT85c0M7egXBCoGYaanRrXkYm9Ube1I7I8dyBL\nsnysP3TyvI6fuXj99vTkqB4+erCrmllen0gvruwBSJ1qvaGZuQWtN0PVGle13gx1bG4hE8/UZnnu\nQJZk+VgvL9fagp4kHX/mYldX+LK8PpFuhD0AqbO0sqZc0H56ywWBllbWIupo52R57kCWZPlYn6+s\ndjS+FVlen0g3wh6A1BkZ6lczDNvGmmGokaH+iDraOVmeO5AlWT7WJ/bt7mh8K7K8PpFuhD0AqVMY\nyGt2alx9uUCD+V3qywWanRpXYSAfdWveZXnuQJZk+VgfGx7U9ORo29j05KjGhge3XTPL6xPpZs65\nqHu4rWKx6EqlUtRtAEiYar2hpZU1jQz1Z+4OO8tzB7Iky8d6ebmm+cqqJvbt7irobZbl9Yn4M7Oz\nzrliR8sQ9gAAAAAg3rYT9ngZJwAAAACkEGEPAAAAAFKIsAcAAAAAKUTYAwAAAIAUIuwBAAAAQAoR\n9gAAuIVqvaFzlVVV642oWwESy8dxVF6u6USpovJyrWc1e90n5w9EbVfUDQAAEFen5i9pZm5BuSBQ\nMww1OzWuIxN7o24LSBQfx9FDJ8/r+JmL129PT47q4aMHY9Un5w/EAVf2AAC4gWq9oZm5Ba03Q9Ua\nV7XeDHVsboFn6IEO+DiOysu1tqAnScefudjVFb5e98n5A3FB2AMA4AaWVtaUC9rvJnNBoKWVtYg6\nApLHx3E0X1ntaHwret0n5w/EBWEPAIAbGBnqVzMM28aaYaiRof6IOgKSx8dxNLFvd0fjW9HrPjl/\nIC4IewAA3EBhIK/ZqXH15QIN5nepLxdodmpchYF81K0BieHjOBobHtT05Gjb2PTkqMaGB2PTJ+cP\nxIU556Lu4baKxaIrlUpRtwEAyKBqvaGllTWNDPXzQA3YJh/HUXm5pvnKqib27e4q6G3W6z45f6CX\nzOysc67Y0TKEPQAAAACIt+2EPV7GCQAAAAApRNgDAAAAgBQi7AEAAABAChH2AAAAACCFCHsAAAAA\nkEKEPQAAAABIIW9hz8z2mdmfmtlnzOzTZvbTrfF3m9klM5tv/fc9vnoAgCQoL9d0olRRebnWs5rV\nekPnKquq1hs9q5lVrEvEmY/908c5yYff+eSifuDXP6nf+eRiz2pyvCNtdnmsfVXSzzrnnjWzQUln\nzexjrZ+91zn3Kx7/NgAkwkMnz+v4mYvXb09Pjurhowe7qnlq/pJm5haUCwI1w1CzU+M6MrG321Yz\niXWJOPOxf/o4J/nwre/+E/3d+lckSX/+uRX98hPP69y7/2lXNTnekUberuw5515wzj3b+ndN0mck\nccQAQEt5udb2oEqSjj9zsatn06v1hmbmFrTeDFVrXNV6M9SxuQWepd4G1iXizMf+6eOc5MPvfHLx\netC75u/Wv9LVFT6Od6TVjrxnz8xeK+nbJH2qNfRTZrZgZr9tZkM3WebtZlYys9KLL764E20CwI6a\nr6x2NL4VSytrygXtp/ZcEGhpZW3bNbOKdYk487F/+jgn+XBq4YWOxreC4x1p5T3smdmApDlJ73DO\nfVHSr0n6RkkTkl6Q9Ks3Ws459z7nXNE5V7zzzjt9twkAO25i3+6OxrdiZKhfzTBsG2uGoUaG+rdd\nM6tYl4gzH/unj3OSD0fH7+5ofCs43pFWXsOemeW0EfR+1zn3h5LknFt2zn3FORdK+k1Jb/DZAwDE\n1djwoKYnR9vGpidHNTY8uO2ahYG8ZqfG1ZcLNJjfpb5coNmpcRUG8t22mzmsS8SZj/3TxznJhwe+\nY79e2feytrFX9r1MD3zH/m3X5HhHWplzzk9hM5P0AUlfcM69Y9P43c65F1r/fqekb3fO/eCtahWL\nRVcqlbz0CQBRKy/XNF9Z1cS+3T17UFWtN7S0sqaRoX4erHSJdYk487F/+jgn+fA7n1zUqYUXdHT8\n7q6C3mYc74gzMzvrnCt2tIzHsPedkp6SdF7Stevi75L0Q9p4CaeT9DlJP3Et/N0MYQ8AAABAlm0n\n7Hn76gXn3NOS7AY/+mNffxMAAAAAsGFHPo0TAAAAALCzCHsAAAAAkEKEPQAAAABIIcIeAAAAAKQQ\nYQ8AAAAAUoiwBwAdqNYbOldZVbXeiLqVHZfluWcV27y3yss1nShVVF6u9axmUrbRk89d1syJc3ry\nucs9q1larOqRJ55XabHas5o+thEQJW9fvQAAaXNq/pJm5haUCwI1w1CzU+M6MrE36rZ2RJbnnlVs\n89566OR5HT9z8frt6clRPXz0YFc1k7KN7n/vJ/TZ5S9Jkv6gtKR7hu/QR9/5pq5qPvDYGT1d3gh5\nj368rENjBT3+4H1d1fSxjYCocWUPALagWm9oZm5B681QtcZVrTdDHZtbiP2z6b2Q5blnFdu8t8rL\ntbYQIUnHn7nY1dWjpGyjJ5+7fD3oXfP88pe6usJXWqxeD3rXPFWudnWFz8c2AuKAsAcAW7C0sqZc\n0H7KzAWBllbWIupo52R57lnFNu+t+cpqR+NbkZRt9MRzyx2Nb8XpC1c6Gt8KH9sIiAPCHgBswchQ\nv5ph2DbWDEONDPVH1NHOyfLcs4pt3lsT+3Z3NL4VSdlG99873NH4Vhw+sKej8a3wsY2AOCDsAcAW\nFAbymp0aV18u0GB+l/pygWanxlUYyEfdmndZnntWsc17a2x4UNOTo21j05OjGhse3HbNpGyjt9x7\nl+4ZvqOH9x5RAAAgAElEQVRt7J7hO/SWe+/ads3i/oIOjRXaxg6NFVTcX7jJErfnYxsBcWDOuah7\nuK1isehKpVLUbQCAqvWGllbWNDLUH7sHVb5lee5ZxTbvrfJyTfOVVU3s292zEJGUbfTkc5f1xHPL\nuv/e4a6C3malxapOX7iiwwf2dBX0NvOxjYBeMbOzzrliR8sQ9gAAAAAg3rYT9ngZJwAAAACkEGEP\nAAAAAFKIsAcAAAAAKUTYAwAAAIAUIuwBAAAAQAoR9gCgA9V6Q+cqq6rWG1G3kgpZXZ9ZnbcvWV6f\nSZm7jz5Li1U98sTzKi1We1YTSJtdUTcAAElxav6SZuYWlAsCNcNQs1PjOjKxN+q2Eiur6zOr8/Yl\ny+szKXP30ecDj53R0+WNkPfox8s6NFbQ4w/e14t2gVThyh4AbEG13tDM3ILWm6Fqjatab4Y6NrcQ\n+2fT4yqr6zOr8/Yly+szKXP30WdpsXo96F3zVLnKFT7gBgh7ALAFSytrygXtp8xcEGhpZS2ijpIt\nq+szq/P2JcvrMylz99Hn6QtXOhoHsoywBwBbMDLUr2YYto01w1AjQ/0RdZRsWV2fWZ23L1len0mZ\nu48+Dx/Y09E4kGWEPQDYgsJAXrNT4+rLBRrM71JfLtDs1LgKA/moW0ukrK7PrM7blyyvz6TM3Uef\nxf0FHRortI0dGiuouL9wkyWA7DLnXNQ93FaxWHSlUinqNgBA1XpDSytrGhnqj92DqiTK6vrM6rx9\nyfL6TMrcffRZWqzq9IUrOnxgD0EPmWBmZ51zxY6WIewBAAAAQLxtJ+zxMk4AAAAASCHCHgAAAACk\nEGEPAAAAAFKIsAcAAAAAKUTYAwAAAIAUIuwBAAAAQAoR9gCkVrXe0LnKqqr1Rqxr+pDluScB6xK9\nkpRjvbRY1SNPPK/SYjXWNYG02RV1AwDgw6n5S5qZW1AuCNQMQ81OjevIxN7Y1fQhy3NPAtYleiUp\nx/oDj53R0+WNQPbox8s6NFbQ4w/eF7uaQBpxZQ9A6lTrDc3MLWi9GarWuKr1ZqhjcwtdPUvto6YP\nWZ57ErAu0StJOdZLi9Xroeyap8rVrq7G+agJpBVhD0DqLK2sKRe0n95yQaCllbVY1fQhy3NPAtYl\neiUpx/rpC1c6Go+qJpBWhD0AqTMy1K9mGLaNNcNQI0P9sarpQ5bnngSsS/RKUo71wwf2dDQeVU0g\nrQh7AFKnMJDX7NS4+nKBBvO71JcLNDs1rsJAPlY1fcjy3JOAdYleScqxXtxf0KGxQtvYobGCivsL\nN1kimppAWplzLuoebqtYLLpSqRR1GwASplpvaGllTSND/T17MO2jpg9ZnnsSsC7RK0k51kuLVZ2+\ncEWHD+zpWSjzUROIMzM765wrdrQMYQ8AAAAA4m07YY+XcQIAAABAChH2AAAAACCFCHsAAAAAkEKE\nPQAAAABIIcIeAAAAAKQQYQ+IkWq9oXOVVVXrjZ7VLC/XdKJUUXm51rOaPvr0wUefSamZ5e2O7EnK\nvpnlPpNS08e5E4jSrqgbALDh1PwlzcwtKBcEaoahZqfGdWRib1c1Hzp5XsfPXLx+e3pyVA8fPRi7\nPn3w0WdSamZ5uyN7krJvZrnPpNT0ce4EosaVPSAGqvWGZuYWtN4MVWtc1Xoz1LG5ha6erSwv19ru\ntCTp+DMXu3q20kefPvjoMyk1s7zdkT1J2Tez3GdSavo4dwJxQNgDYmBpZU25oP1wzAWBllbWtl1z\nvrLa0fhW+OjTBx99JqVmlrc7sicp+2aW+0xKTR/nTiAOCHtADIwM9asZhm1jzTDUyFD/tmtO7Nvd\n0fhW+OjTBx99JqVmlrc7sicp+2aW+0xKTR/nTiAOCHtADBQG8pqdGldfLtBgfpf6coFmp8ZVGMhv\nu+bY8KCmJ0fbxqYnRzU2PBirPn3w0WdSamZ5uyN7krJvZrnPpNT0ce4E4sCcc50tYBZIGnDOfdFP\nSy9VLBZdqVTaqT8HRKZab2hpZU0jQ/09exBQXq5pvrKqiX27e3an5aNPH3z0mZSaWd7uyJ6k7JtZ\n7jMpNX2cO4FeMbOzzrliR8tsJeyZ2e9J+klJX5F0VtIrJT3inPvl7TTaKcIeAAAAgCzbTtjb6ss4\n721dyXubpD+WNCrphzvsDwAAAACwQ7Ya9nJmltNG2DvlnGtK6uz1nwAAAACAHbPVsPcbkj4n6Q5J\np83sNZJ27D17AAAAAIDO7NrKLznnHpX06KahvzGzf+ynJQAAAABAt7Z0Zc/Mhs3st8zsI63b90r6\nEa+dAQAAAAC2basv43y/pI9KenXr9mclvcNHQwAAAACA7m017O1xzn1IUihJzrmr2vgahpsys31m\n9qdm9hkz+7SZ/XRr/FVm9jEzu9D6/1BXMwBwS9V6Q+cqq6rWG7GumRTMvbdzLy1W9cgTz6u0WO1Z\nzV5jm2fz/JGUPpOivFzTiVJF5eVaz2qefLaiBz/w5zr5bKVnNX30CURpS+/Zk/QlMyuo9QmcZnaf\npL+7zTJXJf2sc+5ZMxuUdNbMPibpRyU96Zz7JTP7OUk/J2lmW90DuKVT85c0M7egXBCoGYaanRrX\nkYm9sauZFMy9t3N/4LEzerq8EfIe/XhZh8YKevzB+3rRbs+wzbN5/khKn0nx0MnzOn7m4vXb05Oj\nevjowa5q3veLH9PlL35ZkvRfP/N5vedP/lLPvOu7Y9cnELWtXtn7GUl/JOkbzey/Szou6V/cagHn\n3AvOuWdb/65J+oykvZKOSvpA69c+oI2vcwDQY9V6QzNzC1pvhqo1rmq9GerY3EJXz1L7qJkUzL23\ncy8tVq8HvWueKldjdYWPbZ7N80dS+kyK8nKtLUBJ0vFnLnZ15ezks5XrQe+aF7745a6u8PnoE4iD\nLYW9Vmj7R5K+Q9JPSPpm59zCVv+Imb1W0rdJ+pSkYefcC626L0j6+pss83YzK5lZ6cUXX9zqnwLQ\nsrSyplzQfojngkBLK2uxqpkUzL23cz994UpH41Fgm2fz/JGUPpNivrLa0fhWfPj85Y7Gt8JHn0Ac\n3DLsmdn/fO0/SUck3SPpH0j6vtbYbZnZgKQ5Se9wzm35u/mcc+9zzhWdc8U777xzq4sBaBkZ6lcz\nDNvGmmGokaH+WNVMCube27kfPrCno/EosM2zef5ISp9JMbFvd0fjW/HWg3d1NL4VPvoE4uB2V/a+\n7xb/vfV2xc0sp42g97vOuT9sDS+b2d2tn98t6fPbax3ArRQG8pqdGldfLtBgfpf6coFmp8ZVGMjH\nqmZSMPfezr24v6BDY4W2sUNjBRX3F26yxM5jm2fz/JGUPpNibHhQ05OjbWPTk6MaGx7cds23vW6f\n7n7Fy9vG7n7Fy/W21+3bdk0ffQJxYM45P4XNTBvvyfuCc+4dm8Z/WVJ10we0vMo5d+xWtYrFoiuV\nSl76BNKuWm9oaWVNI0P9PXuw4qNmUjD33s69tFjV6QtXdPjAnlgFvc3Y5tk8fySlz6QoL9c0X1nV\nxL7dPQtQJ5+t6MPnL+utB+/qKuht5qNPoFfM7KxzrtjRMlsNe2b2vZK+WVLftTHn3MO3+P3vlPSU\npPNqfWWDpHdp4317H5I0KumipB9wzn3hVn+bsAcAAAAgy7YT9rb01Qtm9uuSvk7SP5b0mKTvl/Rn\nt1rGOfe0JLvJj9/SQY8AAAAAgA5t9asXvsM5Ny1pxTn385ImJfXmejkAAAAAoOe2Gvaufd7w35vZ\nq7Xxhen7/bQEAAAAAOjWll7GKenDZrZb0qyks62xx/y0BAAAAADo1i3Dnpm9XlLFOfd/t24PaOMD\nV/5S0nv9twcAAAAA2I7bvYzzNyR9WZLM7LCkX2qN/Z2k9/ltDQAAAACwXbcLey/b9LUI/1zS+5xz\nc865fyNpzG9rQPZU6w2dq6yqWm9E3UoqsD6zKQnbPQk9Irt87J8+apYWq3rkiedVWqz2rGZ5uaYT\npYrKy7We1QSidLv37L3MzHY5565q4+sS3t7BsgA6cGr+kmbmFpQLAjXDULNT4zoysTfqthKL9ZlN\nSdjuSegR2eVj//RR84HHzujp8kbIe/TjZR0aK+jxB+/rquZDJ8/r+JmL129PT47q4aMHu6oJRO12\nV/Y+KOm/mdkpbXwi51OSZGZj2ngpJ4AeqNYbmplb0HozVK1xVevNUMfmFnjWf5tYn9mUhO2ehB6R\nXT72Tx81S4vV60HvmqfK1a6u8JWXa21BT5KOP3ORK3xIvFuGPefcL0j6WUnvl/Sdzjm3abl/4bc1\nIDuWVtaUC9oPx1wQaGll7SZL4FZYn9mUhO2ehB6RXT72Tx81T1+40tH4VsxXVjsaB5Litt+z55w7\n45z7L865L20a+6xz7lm/rQHZMTLUr2YYto01w1AjQ/0RdZRsrM9sSsJ2T0KPyC4f+6ePmocP7Olo\nfCsm9u3uaBxIiq1+qToAjwoDec1OjasvF2gwv0t9uUCzU+MqDOSjbi2RWJ/ZlITtnoQekV0+9k8f\nNYv7Czo0VmgbOzRWUHF/4SZL3N7Y8KCmJ0fbxqYnRzU2PLjtmkAc2FdfmRlfxWLRlUqlqNsAvKvW\nG1paWdPIUD8P/nqA9ZlNSdjuSegR2eVj//RRs7RY1ekLV3T4wJ6ugt5m5eWa5iurmti3m6CH2DGz\ns865YkfLEPYAAAAAIN62E/Z4GScAAAAApBBhDwAAAABSiLAHAAAAAClE2AMAAACAFCLsAQAAAEAK\nEfaAlKvWGzpXWVW13oi6lR2X5bmjd9iPeivL6zPLcy8v13SiVFF5udazmqXFqh554nmVFqs9qwmk\nza6oGwDgz6n5S5qZW1AuCNQMQ81OjevIxN6o29oRWZ47eof9qLeyvD6zPPeHTp7X8TMXr9+enhzV\nw0cPdlXzgcfO6OnyRsh79ONlHRor6PEH7+uqJpBGXNkDUqpab2hmbkHrzVC1xlWtN0Mdm1vIxDPK\nWZ47eof9qLeyvD6zPPfycq0t6EnS8WcudnWFr7RYvR70rnmqXOUKH3ADhD0gpZZW1pQL2g/xXBBo\naWUtoo52Tpbnjt5hP+qtLK/PLM99vrLa0fhWnL5wpaNxIMsIe0BKjQz1qxmGbWPNMNTIUH9EHe2c\nLM8dvcN+1FtZXp9ZnvvEvt0djW/F4QN7OhoHsoywB6RUYSCv2alx9eUCDeZ3qS8XaHZqXIWBfNSt\neZfluaN32I96K8vrM8tzHxse1PTkaNvY9OSoxoYHt12zuL+gQ2OFtrFDYwUV9xdusgSQXeaci7qH\n2yoWi65UKkXdBpBI1XpDSytrGhnqz8QDi82yPHf0DvtRb2V5fWZ57uXlmuYrq5rYt7uroLdZabGq\n0xeu6PCBPQQ9ZIKZnXXOFTtahrAHAAAAAPG2nbDHyzgBAAAAIIUIewAAAACQQoQ9AAAAAEghwh4A\nAAAApBBhDwAAAABSiLAHAAAAAClE2ANipFpv6FxlVdV6I+pWbok+4y/LcweQDZzngNvbFXUDADac\nmr+kmbkF5YJAzTDU7NS4jkzsjbqtl6DP+Mvy3AFkA+c5YGu4sgfEQLXe0MzcgtaboWqNq1pvhjo2\ntxC7ZyvpM/6yPHcA2cB5Dtg6wh4QA0sra8oF7YdjLgi0tLIWUUc3Rp/xl+W5A8gGznPA1hH2gBgY\nGepXMwzbxpphqJGh/og6ujH6jL8szx1ANnCeA7aOsAfEQGEgr9mpcfXlAg3md6kvF2h2alyFgXzU\nrbWhz/jL8twBZAPnOWDrzDkXdQ+3VSwWXalUiroNwLtqvaGllTWNDPXH+k6LPuMvy3MHkA2c55A1\nZnbWOVfsZBk+jROIkcJAPhF3WPQZf1meO4Bs4DwH3B4v4wQAAACAFCLsAQAAAEAKEfYAAAAAIIUI\newAAAACQQoQ9AAAAAEghwl7GVOsNnausqlpvRN1K4vlYl0nZPj76LC/XdKJUUXm51rOaPiRlG6F3\nsnysA3Hm437jyecua+bEOT353OWe1QSixFcvZMip+UuamVtQLgjUDEPNTo3ryMTeqNtKJB/rMinb\nx0efD508r+NnLl6/PT05qoePHuy21Z5LyjZC72T5WAfizMf9xv3v/YQ+u/wlSdIflJZ0z/Ad+ug7\n39RVTSBqXNnLiGq9oZm5Ba03Q9UaV7XeDHVsboFnlbfBx7pMyvbx0Wd5udZ2hy1Jx5+5GLsrfEnZ\nRuidLB/rQJz5uN948rnL14PeNc8vf4krfEg8wl5GLK2sKRe0b+5cEGhpZS2ijpLLx7pMyvbx0ed8\nZbWj8agkZRuhd7J8rANx5uN+44nnljsaB5KCsJcRI0P9aoZh21gzDDUy1B9RR8nlY10mZfv46HNi\n3+6OxqOSlG2E3snysQ7EmY/7jfvvHe5oHEgKwl5GFAbymp0aV18u0GB+l/pygWanxlUYyEfdWuL4\nWJdJ2T4++hwbHtT05Gjb2PTkqMaGB7ttt6eSso3QO1k+1oE483G/8ZZ779I9w3e0jd0zfIfecu9d\n264JxIE556Lu4baKxaIrlUpRt5EK1XpDSytrGhnq58FFl3ysy6RsHx99lpdrmq+samLf7tgFvc2S\nso3QO1k+1oE483G/8eRzl/XEc8u6/95hgh5ix8zOOueKHS1D2AMAAACAeNtO2ONlnAAAAACQQoQ9\nAAAAAEghwh4AAAAApBBhDwAAAABSiLAHAAAAAClE2AMAAACAFPIW9szst83s82b2PzaNvdvMLpnZ\nfOu/7/H194EkqtYbOldZVbXe6FnN0mJVjzzxvEqL1Z7V9NFnUmomRVLmnpQ+AXSnvFzTiVJF5eVa\nz2r6uH8D0sbb9+yZ2WFJdUnHnXPf0hp7t6S6c+5XOqnF9+whC07NX9LM3IJyQaBmGGp2alxHJvZ2\nVfOBx87o6fJX7wQPjRX0+IP3xa7PpNRMiqTMPSl9AujOQyfP6/iZi9dvT0+O6uGjB7uq6eP+DYi7\nWH3PnnPutKQv+KoPpEm13tDM3ILWm6Fqjatab4Y6NrfQ1dWO0mK17Y5Qkp4qV7t6BtRHn0mpmRRJ\nmXtS+gTQnfJyrS3oSdLxZy52dYXPx/0bkFZRvGfvp8xsofUyz6Gb/ZKZvd3MSmZWevHFF3eyP2DH\nLa2sKRe0H465INDSytq2a56+cKWj8a3w0WdSaiZFUuaelD4BdGe+strR+Fb4uH8D0mqnw96vSfpG\nSROSXpD0qzf7Refc+5xzRedc8c4779yp/oBIjAz1qxmGbWPNMNTIUP+2ax4+sKej8a3w0WdSaiZF\nUuaelD4BdGdi3+6OxrfCx/0bkFY7Gvacc8vOua8450JJvynpDTv594G4KgzkNTs1rr5coMH8LvXl\nAs1OjaswkN92zeL+gg6NFdrGDo0VVNxfuMkS0fSZlJpJkZS5J6VPAN0ZGx7U9ORo29j05KjGhge3\nXdPH/RuQVt4+oEWSzOy1kj686QNa7nbOvdD69zslfbtz7gdvV4cPaEFWVOsNLa2saWSov2cPekuL\nVZ2+cEWHD+zp2R2hjz6TUjMpkjL3pPQJoDvl5ZrmK6ua2Le7q6C3mY/7NyDOtvMBLT4/jfODkt4k\naY+kZUn/tnV7QpKT9DlJP3Et/N0KYQ8AAABAlm0n7O3y1Yxz7oduMPxbvv4eAAAAAOCrovg0TgAA\nAACAZ4Q9AAAAAEghwh4AAAAApBBhDwAAAABSiLAHAAAAAClE2EPXqvWGzlVWVa03Yl0zCXzM+8nn\nLmvmxDk9+dzlntUsLVb1yBPPq7RY7VnN8nJNJ0oVlZdrPavpQ1b3TSnbcwcQPzz+AG7P21cvIBtO\nzV/SzNyCckGgZhhqdmpcRyb2xq5mEviY9/3v/YQ+u/wlSdIflJZ0z/Ad+ug739RVzQceO6Onyxsh\n79GPl3VorKDHH7yvq5oPnTyv42cuXr89PTmqh48e7KqmD1ndN6Vszx1A/PD4A9garuxh26r1hmbm\nFrTeDFVrXNV6M9SxuYWung3zUTMJfMz7yecuXw961zy//KWurvCVFqvXg941T5WrXV3hKy/X2oKe\nJB1/5mLsrvBldd+Usj13APHD4w9g6wh72LallTXlgvZdKBcEWlpZi1XNJPAx7yeeW+5ofCtOX7jS\n0fhWzFdWOxqPSlb3TSnbcwcQPzz+ALaOsIdtGxnqVzMM28aaYaiRof5Y1UwCH/O+/97hjsa34vCB\nPR2Nb8XEvt0djUclq/umlO25A4gfHn8AW0fYw7YVBvKanRpXXy7QYH6X+nKBZqfGVRjIx6pmEviY\n91vuvUv3DN/RNnbP8B16y713bbtmcX9Bh8YKbWOHxgoq7i/cZInbGxse1PTkaNvY9OSoxoYHt13T\nh6zum1K25w4gfnj8AWydOeei7uG2isWiK5VKUbeBm6jWG1paWdPIUH/PToo+aiaBj3k/+dxlPfHc\nsu6/d7iroLdZabGq0xeu6PCBPV0Fvc3KyzXNV1Y1sW937ILeZlndN6Vszx1A/PD4A1ljZmedc8WO\nliHsAQAAAEC8bSfs8TJOAAAAAEghwh4AAAAApBBhDwAAAABSiLAHAAAAAClE2AMAAACAFCLsIZbK\nyzWdKFVUXq71pF613tC5yqqq9UZP6vlSWqzqkSeeV2mx2rOaPuaelJo+JKVPAIgTzp1ANHZF3QDw\ntR46eV7Hz1y8fnt6clQPHz247Xqn5i9pZm5BuSBQMww1OzWuIxN7e9FqTz3w2Bk9Xd4IeY9+vKxD\nYwU9/uB9XdX0Mfek1PQhKX0CQJxw7gSiw5U9xEp5udYW9CTp+DMXt32Fr1pvaGZuQevNULXGVa03\nQx2bW4jdM4ulxer1oHfNU+VqV1f4fMw9KTV9SEqfABAnnDuBaBH2ECvzldWOxm9naWVNuaB9N88F\ngZZW1rZVz5fTF650NL4VPuaelJo+JKVPAIgTzp1AtAh7iJWJfbs7Gr+dkaF+NcOwbawZhhoZ6t9W\nPV8OH9jT0fhW+Jh7Umr6kJQ+ASBOOHcC0SLsIVbGhgc1PTnaNjY9Oaqx4cFt1SsM5DU7Na6+XKDB\n/C715QLNTo2rMJDvRbs9U9xf0KGxQtvYobGCivsLN1ni9nzMPSk1fUhKnwAQJ5w7gWiZcy7qHm6r\nWCy6UqkUdRvYQeXlmuYrq5rYt3vbQW+zar2hpZU1jQz1x/oOprRY1ekLV3T4wJ6ugt5mPuaelJo+\nJKVPAIgTzp1A98zsrHOu2NEyhD0AAAAAiLfthD1exgkAAAAAKUTYAwAAAIAUIuwBAAAAQAoR9gAA\nAAAghQh7AAAAAJBChD0AAAAASCHCXoxV6w2dq6yqWm/0rGZ5uaYTpYrKy7We1fSh13P3Me+Tz1b0\n4Af+XCefrfSspo8+S4tVPfLE8yotVntW0wcf+7sPSekTAACA79mLqVPzlzQzt6BcEKgZhpqdGteR\nib1d1Xzo5HkdP3Px+u3pyVE9fPRgt632XK/n7mPe9/3ix3T5i1++fvvuV7xcz7zru7uq6aPPBx47\no6fLXw15h8YKevzB+7qq6YOP/d2HpPQJAADSh+/ZS4lqvaGZuQWtN0PVGle13gx1bG6hqysJ5eVa\nW5CQpOPPXIzdFb5ez93HvE8+W2kLepL0whe/3NUVPh99lharbUFPkp4qV2N3hc/H/u5DUvoEAAC4\nhrAXQ0sra8oF7ZsmFwRaWlnbds35ympH41Hp9dx9zPvD5y93NL4VPvo8feFKR+NR8bG/+5CUPgEA\nAK4h7MXQyFC/mmHYNtYMQ40M9W+75sS+3R2NR6XXc/cx77cevKuj8a3w0efhA3s6Go+Kj/3dh6T0\nCQAAcA1hL4YKA3nNTo2rLxdoML9LfblAs1PjKgzkt11zbHhQ05OjbWPTk6MaGx7stt2e6vXcfcz7\nba/bp7tf8fK2sbtf8XK97XX7tl3TR5/F/QUdGiu0jR0aK6i4v3CTJaLhY3/3ISl9AgAAXMMHtMRY\ntd7Q0sqaRob6e/aAsrxc03xlVRP7dscu6G3W67n7mPfJZyv68PnLeuvBu7oKepv56LO0WNXpC1d0\n+MCe2AW9zXzs7z4kpU8AAJAu2/mAFsIeAAAAAMQcn8YJAAAAAJBE2AMAAACAVCLsAQAAAEAKEfYA\nAAAAIIUIewAAAACQQoS9GKvWGzpXWVW13oi6lR3X67mXl2s6UaqovFzrST1JevK5y5o5cU5PPne5\nZzV99FlarOqRJ55XabEa65rs7/Gfe1L67LWszhsAkHx89UJMnZq/pJm5BeWCQM0w1OzUuI5M7I26\nrR3R67k/dPK8jp+5eP329OSoHj56sKse73/vJ/TZ5S9dv33P8B366Dvf1FVNH30+8NgZPV3+aiA7\nNFbQ4w/eF7ua7O/xn3tS+uy1rM4bABA/fPXC/9/e/UdHdd53Hv88gwZJAWHkwRGqJRUaKcQkyCrV\npsgpHDuuqZ1Qw1bNad2y2Ocsx+42Pafr3dQQ1+s6Xtt1lV3j7NltS+qmNnZ/pEWpSCAJsJQsdoDU\niiokRwRrGjmMCBqbsQSSK4mBefYPTWAGI5iZe6/mzuj9OicH3S8zXz3PfXTDfPxo7hSJ2NiktrT3\naCKe0OjkBU3EE3qkvWdW/Fdlt+cejo6mBShJ2nHkpKOdswN9Q2lBT5JORN9ztMPnxTg7B2JpoUyS\nXg3HHO3GedGTn3f/z71Qxum22TpvAEDxIOz50ODwuIKB9KUJBgIaHB7P04hmjttz746MZFXPxL6+\naFb1THgxzkP9Z7Kq56snP+/+n3uhjNNts3XeAIDiQdjzoZrKcsUTibRaPJFQTWV5nkY0c9yee1Pt\nwqzqmVi7vCqreia8GOeahkVZ1fPVk593/8+9UMbpttk6bwBA8SDs+VBofqnaWhtVFgyoorREZcGA\n2lobFZpfmu+hec7tuddXVWhTS11abVNLneqrKnIe453LF2tZ1by02rKqebpz+eKce3oxzualIa2u\nD6XVVteH1Lw0NM0z8tOTn3f/z71Qxum22TpvAEDx4AYtPhYbm9Tg8LhqKstn3YsLt+cejo6qOzKi\nplTYS64AACAASURBVNqFjgJUqgN9Q9rXF9Xa5VWOgl4qL8bZORDTof4zWtOwyFEo87onP+/+n3uh\njNNts3XeAAB/yeUGLYQ9AAAAAPA57sYJAAAAAJBE2AMAAACAokTYAwAAAIAiRNgDAAAAgCJE2AMA\nAACAIkTYAwAAAIAi5FnYM8Z8xRjztjHmjZTajcaY/caY/uSflV59/2IQG5vUsciIYmOTrvUMR0e1\nszOicHTU1z3dnrsX57JzIKbn9p1Q50DMtZ4dXRFtful1dXRFXOtZCOsjeTNOL3gxdwAAAC949jl7\nxpg1ksYk7bDWfixZa5P0rrX2WWPMVkmV1tot1+s1Gz9nb1f3KW1p71EwEFA8kVBba6PubbrZUc/H\nO3q14+jJS8ebWur05PoVvuvp9ty9OJcbXziq18KXQ97q+pBe3rzKUc9Vz+zX0Lnzl46rF8zVkUfv\nctSzENbHq3F6wYu5AwAAZMJXn7NnrT0k6d0ryuslvZT8+iVJG7z6/oUsNjapLe09mognNDp5QRPx\nhB5p73G0kxCOjqa9mJakHUdOOtpF8aKn23P34lx2DsTSgp4kvRqOOdrh6+iKpAU9STp97ryjHb5C\nWB+vxukFL+YOAADgpZl+z16Vtfa0JCX//OB0DzTGPGiM6TTGdL7zzjszNkA/GBweVzCQvjTBQECD\nw+M59+yOjGRVz1dPt+fuxbk81H8mq3omdvcOZVXPRCGsz7XG42ScXvBi7gAAAF7y7Q1arLVfttY2\nW2ubb7rppnwPZ0bVVJYrnkik1eKJhGoqy3Pu2VS7MKt6vnq6PXcvzuWahkVZ1TOxbsXirOqZKIT1\nudZ4nIzTC17MHQAAwEszHfaixphqSUr++fYMf/+CEJpfqrbWRpUFA6ooLVFZMKC21kaF5pfm3LO+\nqkKbWurSapta6lRfVeGrnm7P3Ytz2bw0pNX1obTa6vqQmpeGpnnG9W1YWavqBXPTatUL5mrDytqc\nexbC+ng1Ti94MXcAAAAveXaDFkkyxiyRtDvlBi1flBRLuUHLjdbaR67XZzbeoEWaeo/Q4PC4airL\nXXtBGY6Oqjsyoqbaha69mPaip9tz9+Jcdg7EdKj/jNY0LHIU9FJ1dEW0u3dI61YsdhT0UhXC+kje\njNMLXswdAADgenK5QYuXd+P8W0m3S1okKSrpjyR1SPp7SXWSTkr6jLX2ypu4vM9sDXsAAAAAIOUW\n9kq8Goy19r5p/upOr74nAAAAAGCKb2/QAgAAAADIHWEPAAAAAIoQYQ8AAAAAihBhDwAAAACKEGHP\nJbGxSR2LjCg2NulazwN9Q9qy85gO9A251rOjK6LNL72ujq6Iaz07B2J6bt8JdQ7EXOvptu0H+3XP\nlw5p+8F+13qGo6Pa2RlRODrqWk8v1seLcXqx5l5cQwAAALOZp5+z5xa/f/TCru5T2tLeo2AgoHgi\nobbWRt3bdLOjnmu3fUdvRt+7dLysap72Pny7o56rntmvoXPnLx1XL5irI4/e5ajnxheO6rXw5Rf8\nq+tDennzKkc93XbLY9/U+IXLP+flJUbHn/qUo56Pd/Rqx9GTl443tdTpyfUrHPX0Yn28GKcXa+7F\nNQQAAFBMcvnoBXb2HIqNTWpLe48m4gmNTl7QRDyhR9p7HO1OHOgbSgt6knQi+p6jHb6OrkhakJCk\n0+fOO9pB6hyIpb3ol6RXwzFf7fBtP9ifFvQkafyCdbTDF46OpgUoSdpx5KSjnTMv1seLcXqx5l5c\nQwAAACDsOTY4PK5gIP00BgMBDQ6P59xzX180q3omdvdePShOV8/Eof4zWdXzoaPndFb1THRHRrKq\nZ8KL9fFinF6suRfXEAAAAAh7jtVUliueSKTV4omEairLc+65dnlVVvVMrFuxOKt6JtY0LMqqng8b\nGquzqmeiqXZhVvVMeLE+XozTizX34hoCAAAAYc+x0PxStbU2qiwYUEVpicqCAbW1Nio0vzTnnncu\nX6xlVfPSasuq5unO5bm/8N+wslbVC+am1aoXzNWGlbU592xeGtLq+lBabXV9SM1LQ9M8Y+Y9dEeD\nyktMWq28xOihOxpy7llfVaFNLXVptU0tdaqvqsi5pxfr48U4vVhzL64hAAAAcIMW18TGJjU4PK6a\nynLXXqQe6BvSvr6o1i6vchT0UnV0RbS7d0jrVix2FCRSdQ7EdKj/jNY0LPJV0Eu1/WC/OnpOa0Nj\ntaOglyocHVV3ZERNtQsdBahUXqyPF+P0Ys29uIYAAACKRS43aCHsAQAAAIDPcTdOAAAAAIAkwh4A\nAAAAFCXCHgAAAAAUIcIeAAAAABQhwh4AAAAAFCHCHgAAAAAUoVkZ9mJjkzoWGVFsbNLXPcPRUe3s\njCgcHXWt5yuHB/SZPz+sVw4PuNazEObuxbnsHIjpuX0n1DkQ83VPL+buxZoDAADAXbPuc/Z2dZ/S\nlvYeBQMBxRMJtbU26t6mm33X8/GOXu04evLS8aaWOj25foWjnrc+8W2dnbh46fiGsjk69sTdjnoW\nwty9OJcbXziq18KXA9nq+pBe3rzKdz29mLsXaw4AAIBr43P2riM2Nqkt7T2aiCc0OnlBE/GEHmnv\ncbQ74UXPcHQ07QW6JO04ctLRzswrhwfSgp4knZ246GiHrxDm7sW57ByIpYUySXo1HHO0G+dFTy/m\n7sWaAwAAwBuzKuwNDo8rGEifcjAQ0ODwuK96dkdGsqpnYlfP6azqmSiEuXtxLg/1n8mqnq+eXszd\nizUHAACAN2ZV2KupLFc8kUirxRMJ1VSW+6pnU+3CrOqZWN9YnVU9E4Uwdy/O5ZqGRVnV89XTi7l7\nseYAAADwxqwKe6H5pWprbVRZMKCK0hKVBQNqa21UaH6pr3rWV1VoU0tdWm1TS53qqypy7rnxtqW6\noWxOWu2GsjnaeNvSnHsWwty9OJfNS0NaXR9Kq62uD6l5aWiaZ+Snpxdz92LNAQAA4I1Zd4MWaep9\nR4PD46qpLHftRaoXPcPRUXVHRtRUu9DRC/RUrxwe0K6e01rfWO0o6KUqhLl7cS47B2I61H9GaxoW\nOQplXvf0Yu5erDkAAACml8sNWmZl2AMAAACAQsLdOAEAAAAAkgh7AAAAAFCUCHsAAAAAUIQIewAA\nAABQhAh7AAAAAFCECHsAAAAAUIQIez7W0RXR5pdeV0dXxNc9w9FR7eyMKBwdda1nbGxSxyIjio1N\nutLvQN+Qtuw8pgN9Q67080rnQEzP7TuhzoGYaz3dPpezHecTAAAUCj5nz6dWPbNfQ+fOXzquXjBX\nRx69y3c9H+/o1Y6jJy8db2qp05PrVzjquav7lLa09ygYCCieSKittVH3Nt2cc7+1276jN6PvXTpe\nVjVPex++3dEYvbDxhaN6LXw55K2uD+nlzasc9XT7XM52nE8AAJAvfM5ekejoiqSFMkk6fe68o904\nL3qGo6NpQU+Sdhw56WiHLzY2qS3tPZqIJzQ6eUET8YQeae/JeRflQN9QWtCTpBPR93y3w9c5EEsL\nepL0ajjmaIfP7XM523E+AQBAoSHs+dDu3qsHkenq+erZHRnJqp6JweFxBQPpP5bBQECDw+M59dvX\nF82qni+H+s9kVc+E2+dytuN8AgCAQkPY86F1KxZnVc9Xz6bahVnVM1FTWa54IpFWiycSqqksz6nf\n2uVVWdXzZU3DoqzqmXD7XM52nE8AAFBoCHs+tGFlraoXzE2rVS+Yqw0ra33Vs76qQpta6tJqm1rq\nVF9VkXPP0PxStbU2qiwYUEVpicqCAbW1Nio0vzSnfncuX6xlVfPSasuq5unO5bmHXC80Lw1pdX0o\nrba6PqTmpaFpnnF9bp/L2Y7zCQAACg03aPGxjq6IdvcOad2KxY5Cmdc9w9FRdUdG1FS70FHQSxUb\nm9Tg8LhqKstdeTF9oG9I+/qiWru8yndBL1XnQEyH+s9oTcMiR0EvldvncrbjfAIAgHzI5QYthD0A\nAAAA8DnuxgkAAAAAkETYAwAAAICiRNgDAAAAgCJE2AMAAACAIkTYAwAAAIAiNCvDXmxsUsciI4qN\nTbrWs6Mros0vva6OrohrPV85PKDP/PlhvXJ4wLWe2w/2654vHdL2g/2u9TzQN6QtO4/pQN+Qaz3D\n0VHt7IwoHB31ZT9p6mMSntt3Qp0DMV/39GLuXvQEAACAu2bdRy/s6j6lLe09CgYCiicSamtt1L1N\nNzvqueqZ/Ro6d/7ScfWCuTry6F2Oet76xLd1duLipeMbyubo2BN3O+p5y2Pf1PiFy+tdXmJ0/KlP\nOeq5dtt39Gb0vUvHy6rmae/Dtzvq+XhHr3YcPXnpeFNLnZ5cv8I3/SRp4wtH9Vr4ciBbXR/Sy5tX\n+a6nF3P3oicAAACujY9euI7Y2KS2tPdoIp7Q6OQFTcQTeqS9x9EOX0dXJC3oSdLpc+cd7fC9cngg\nLehJ0tmJi452+LYf7E8LepI0fsE62uE70DeUFvQk6UT0PUc7fOHoaFqQkKQdR07mvIPkdj9pavct\nNZRJ0qvhmKPdOC96ejF3L3oCAADAG7Mq7A0OjysYSJ9yMBDQ4PB4zj1391492ExXz8SuntNZ1TPR\nMc1zp6tnYl9fNKt6JrojI1nVZ7qfJB3qP5NVPV89vZi7Fz0BAADgjVkV9moqyxVPJNJq8URCNZXl\nOfdct2JxVvVMrG+szqqeiQ3TPHe6eibWLq/Kqp6JptqFWdVnup8krWlYlFU9Xz29mLsXPQEAAOCN\nWRX2QvNL1dbaqLJgQBWlJSoLBtTW2qjQ/NKce25YWavqBXPTatUL5mrDytqce268baluKJuTVruh\nbI423rY0554P3dGg8hKTVisvMXrojoace965fLGWVc1Lqy2rmqc7l+cedOurKrSppS6ttqmlTvVV\nFb7oJ0nNS0NaXR9Kq62uD6l5aWiaZ+Snpxdz96InAAAAvDHrbtAiTb13b3B4XDWV5Y6CXqqOroh2\n9w5p3YrFjoJeqlcOD2hXz2mtb6x2FPRSbT/Yr46e09rQWO0o6KU60DekfX1RrV1e5SjopQpHR9Ud\nGVFT7UJXgoTb/aSp99kd6j+jNQ2LHIUyr3t6MXcvegIAAGB6udygZVaGPQAAAAAoJNyNEwAAAAAg\nibAHAAAAAEWJsAcAAAAARYiwBwAAAABFiLAHAAAAAEWIsAcAAAAARSgvYc8Y85YxptcY022MmfHP\nVOgciOm5fSfUORBzrWdHV0SbX3pdHV0R13o+9rVj+oX/vk+Pfe2Yr3t6cT696Om22NikjkVGFBub\ndK2nF/P2YpwAAADwv7x8zp4x5i1JzdbaM5k83s3P2dv4wlG9Fr78Qnp1fUgvb17lqOeqZ/Zr6Nz5\nS8fVC+bqyKN3Oeq5ZOue99XeevbTvuvpxfn0oqfbdnWf0pb2HgUDAcUTCbW1Nureppsd9fRi3l6M\nEwAAADOPz9m7js6BWNqLaUl6NRxztIvS0RVJC3qSdPrceUc7fNPtujnZjfOipxfn04uebouNTWpL\ne48m4gmNTl7QRDyhR9p7HO2ceTFvL8YJAACAwpGvsGcl7TPGfN8Y8+DVHmCMedAY02mM6XznnXdc\n+aaH+q++kThdPRO7e4eyqmfiWz+IZlXPV08vzqcXPd02ODyuYCD90gkGAhocHs+5pxfz9mKcAAAA\nKBz5CnufsNaulHSPpM8aY9Zc+QBr7Zettc3W2uabbrrJlW+6pmFRVvVMrFuxOKt6Ju75aFVW9Xz1\n9OJ8etHTbTWV5YonEmm1eCKhmsrynHt6MW8vxgkAAIDCkZewZ639SfLPtyX9o6SPz8T3bV4a0ur6\nUFptdX1IzUtD0zzj+jasrFX1grlpteoFc7VhZW3OPZ/6tVuzquerpxfn04uebgvNL1Vba6PKggFV\nlJaoLBhQW2ujQvNLc+7pxby9GCcAAAAKx4zfoMUYM09SwFo7mvx6v6QnrbXfnu45bt6gRZp6f9Sh\n/jNa07DItRDR0RXR7t4hrVux2FHQS/XY147pWz+I6p6PVjkKZV739OJ8etHTbbGxSQ0Oj6umsty1\nAOXFvL0YJwAAAGZWLjdoyUfY+zlN7eZJUomkv7HWPn2t57gd9gAAAACgkOQS9kq8Gsx0rLU/kuTO\nlhIAAAAA4Kpm1UcvAAAAAMBsQdgDAAAAgCJE2AMAAACAIkTYAwAAAIAiNCvD3oG+IW3ZeUwH+oZc\n67n9YL/u+dIhbT/Y71rPB144og//4R498MIR13o++OL39JHHvqkHX/yeaz07uiLa/NLr6uiKuNYz\nNjapY5ERxcYmXevptnB0VDs7IwpHR13rWQjzBgAAQGGY8Y9eyIWbH72wdtt39Gb0vUvHy6rmae/D\ntzvqectj39T4hcvnsbzE6PhTn3LUc8nWPe+rvfXsp33Xc9Uz+zV07vyl4+oFc3Xk0bsc9dzVfUpb\n2nsUDAQUTyTU1tqoe5tudtTTbY939GrH0ZOXjje11OnJ9Ssc9SyEeQMAACA/cvnohVm1s3egbygt\n6EnSieh7jnb4th/sTwt6kjR+wTra4ZtuJ8/JDt90O3lOdvg6uiJpQU+STp8772iHLzY2qS3tPZqI\nJzQ6eUET8YQeae/x1U5XODqaFvQkaceRk452+Aph3gAAACgssyrs7euLZlXPREfP6azqmTg88G5W\n9UwcCseyqmdid+/VQ/J09UwMDo8rGEj/sQwGAhocHs+5p9u6IyNZ1TNRCPMGAABAYZlVYW/t8qqs\n6pnY0FidVT0Tty29Mat6JtbUh7KqZ2LdisVZ1TNRU1mueCKRVosnEqqpLM+5p9uaahdmVc9EIcwb\nAAAAhWVWhb07ly/Wsqp5abVlVfN05/Lcw8lDdzSovMSk1cpLjB66oyHnni9ubsmqnokvP/CLWdUz\nsWFlraoXzE2rVS+Yqw0ra3PuGZpfqrbWRpUFA6ooLVFZMKC21kaF5pfm3NNt9VUV2tRSl1bb1FKn\n+qqKnHsWwrwBAABQWGbdDVqkqffu7euLau3yKkdBL9X2g/3q6DmtDY3VjoJeqgdeOKLDA+/qtqU3\nOgp6qR588Xs6FI5pTX3IUdBL1dEV0e7eIa1bsdhR0EsVG5vU4PC4airLfRt4wtFRdUdG1FS70FHQ\nS1UI8wYAAMDMy+UGLbMy7AEAAABAIeFunAAAAAAASYQ9AAAAAChKhD0AAAAAKEKEPQAAAAAoQoQ9\nAAAAAChChD0AAAAAKEIFEfYGh8d1oG/ItX7b9h7X7V88qG17j7vW80DfkLbsPObqODu6Itr80uvq\n6Iq41vNzX+3SrU98W5/7apdrPb2Ye2xsUsciI4qNTfqyHwAAAOB3BfE5e6XVDbb6/ue1rGqe9j58\nu6NeDZ/fo3jKlING6v/jTzvquXbbd/Rm9L1Lx26Mc9Uz+zV07vyl4+oFc3Xk0bsc9Vyydc/7am89\n67+57+o+pS3tPQoGAoonEmprbdS9TTf7ph8AAAAw04r+c/ZORN9ztHu0be/xtKAnSXErRzt8B/qG\n0sKO5HycHV2RtKAnSafPnXe0wzfdTp6THT4v5h4bm9SW9h5NxBManbygiXhCj7T35Lwj53Y/AAAA\noFAUVNiTpH190Zyfu6vn6iFkunomphuPk3Hu7r36eKarZ2L/8bezqmfCi7kPDo8rGEj/sQwGAhoc\nHvdFPwAAAKBQFFzYW7u8Kufnrm9cnFU9E9ONx8k41624+nimq2firls+mFU9E17MvaayXPFEIq0W\nTyRUU1nui34AAABAoSiosLesap7uXJ574Hn4V25R0KTXgmaqnqs7ly/Wsqp5aTWn49ywslbVC+am\n1aoXzNWGlbU59/wfv7Eyq3omvJh7aH6p2lobVRYMqKK0RGXBgNpaGxWaX+qLfgAAAEChKIgbtCz+\n0EftX3/jgKMQkWrb3uPa1TOk9Y2LHQW9VAf6hrSvL6q1y6tcG2dHV0S7e4e0bsViR0Ev1ee+2qX9\nx9/WXbd80FHQS+XF3GNjkxocHldNZbkrwcztfgAAAMBMyuUGLQUR9pqbm21nZ2e+hwEAAAAAeVH0\nd+MEAAAAAGSGsAcAAAAARYiwBwAAAABFiLAHAAAAAEWIsAcAAAAARYiwBwAAAABFqCDCXu+ps/rM\nn77qWr/f2v5d1T+6R7+1/buu9Xz6G2+o5Y//r57+xhuu9ewciOm5fSfUORBzref2g/2650uHtP1g\nv2s9w9FR7eyMKBwdda2n22JjkzoWGVFsbDLfQwEAAABmREF8zl5pdYOtvv95SdJbz37aUa8lW/e8\nr+a0589t3aNEynFA0o8c9tz4wlG9Fr4c8lbXh/Ty5lWOet7y2Dc1fuHyepeXGB1/6lOOej7e0asd\nR09eOt7UUqcn169w1NNtu7pPaUt7j4KBgOKJhNpaG3Vv0835HhYAAACQsVnxOXtOdvim28lzssP3\n9DfeSAt6kpRI1nPVORBLC3qS9Go45miHb/vB/rSgJ0njF6yjHb5wdDQt6EnSjiMnfbXDFxub1Jb2\nHk3EExqdvKCJeEKPtPewwwcAAICiV3BhrytyLufn/vOPR7KqZ2L3G0NZ1TNxqP9MVvVMdPSczqqe\nie7I1c/bdPV8GBweVzCQ/mMeDAQ0ODyepxEBAAAAM6Pgwt7K2gU5P/fjP7swq3om1n1scVb1TKxp\nWJRVPRMbGquzqmeiqfbq5226ej7UVJYrnkjfe40nEqqpLM/TiAAAAICZUXBh7x9+d3XOz/2bhz6R\nVT0Tf/irH3vfSQwk67lqXhrS6vpQWm11fUjNS0PTPOP6HrqjQeUlJq1WXmL00B0NOfesr6rQppa6\ntNqmljrVV1Xk3NNtofmlamttVFkwoIrSEpUFA2prbVRofmm+hwYAAAB4qmBu0HLvf/uKo6CX6re2\nf1f//OMRffxnFzoKeqme/sYb2v3GkNZ9bLGjoJeqcyCmQ/1ntKZhkaOgl2r7wX519JzWhsZqR0Ev\nVTg6qu7IiJpqF/oq6KWKjU1qcHhcNZXlBD0AAAAUnFxu0FIQYa+5udl2dnbmexgAAAAAkBez4m6c\nAAAAAIDrI+wBAAAAQBEi7AEAAABAESLsAQAAAEARIuwBAAAAQBEqiLDXe+qsln1+j2v9Vn7hW1qy\ndY9WfuFbrvV85fCAPvPnh/XK4QHXem7be1y3f/Ggtu097lrPcHRUOzsjCkdHXesJAAAAwH8K4qMX\nSqsbbPX9z0uS3nr20456Ldn6/tDotOetT3xbZycuXjq+oWyOjj1xt6OeDZ/fo3jK0gSN1P/Hzsb5\neEevdhw9eel4U0udnly/wlFPAAAAAN6bFR+94GSHb7qdPCc7fK8cHkgLepJ0duKiox2+bXuPpwU9\nSYpbOdrhC0dH04KeJO04cpIdPgAAAKBIFVzYm3SwEfnueCKreiZ29ZzOqp5Zz6Gs6pnojoxkVQcA\nAABQ2Aou7JWa3J97Y/nVpztdPRPrG6uzqmfWc3FW9Uw01S7Mqg4AAACgsBVc2Dvh4H1rXX90T1b1\nTGy8baluKJuTVruhbI423rY0554P/8otCl4RaoNmqp6r+qoKbWqpS6ttaqlTfVVFzj0BAAAA+FfB\n3KBlyQPPOwp6qVZ+4Vt6dzyhG8sDjoJeqlcOD2hXz2mtb6x2FPRSbdt7XLt6hrS+cbGjoJcqHB1V\nd2RETbULCXoAAABAgcjlBi0FEfaam5ttZ2dnvocBAAAAAHkxK+7GCQAAAAC4PsIeAAAAABQhwh4A\nAAAAFCHCHgAAAAAUIcIeAAAAABQhwh4AAAAAFKG8hD1jzN3GmBPGmLAxZuv1Ht976qyWbN3j2vdf\n9dQ+Ldm6R6ue2udaz86BmJ7bd0KdAzHXeoajo9rZGVE4OurrnrGxSR2LjCg2NulaTwAAAADOzPjn\n7Blj5kh6U9JdkgYlvS7pPmtt33TPKa1usNX3Py9JeutZZx+sfrXQ6LTnxheO6rXw5ZC3uj6klzev\nctTz8Y5e7Th68tLxppY6Pbl+he967uo+pS3tPQoGAoonEmprbdS9TTc76gkAAAAgXaF8zt7HJYWt\ntT+y1p6X9HeS1mf6ZCc7fNPt5DnZ4esciKUFPUl6NRxztMMXjo6mhTJJ2nHkpKPdOC96xsYmtaW9\nRxPxhEYnL2gintAj7T3s8AEAAAA+kI+dvV+XdLe1dnPy+D9I+kVr7e9d8bgHJT0oSZpT8gtzb1oy\n9RdWOh8Nfz+X7z23qv4XZK7yFw56zqlY9DNz5i2svrJ+8b2R0xdHz/wkl56BD9wQKlnw0wlfduHc\nO28l/u1sTinSi54mWP6BksrqD5tAYM7FfzurOR+4QTaRuHhh+PSbNj7+b7n0hGcWSTqT70Hgmlgj\nf2N9/I818j/WyN9YH/9bZq2tyOYJJV6N5BqmiVtXFKz9sqQvS5IxpnPydH9WW5aYWcaYzgtn32aN\nfMoY05nttj9mFmvkb6yP/7FG/sca+Rvr43/GmM5sn5OPX+MclFSbclwjKacdMAAAAADA1eUj7L0u\nqcEYs9QYM1fSb0r6eh7GAQAAAABFa8Z/jdNae8EY83uS9kqaI+kr1tofXOdpX/Z+ZHCINfI31sf/\nWCN/Y338jzXyP9bI31gf/8t6jWb8Bi0AAAAAAO/l5UPVAQAAAADeIuwBAAAAQBHyddgzxtxtjDlh\njAkbY7bmezx4P2PMW8aYXmNMdy63g4X7jDFfMca8bYx5I6V2ozFmvzGmP/lnZT7HONtNs0ZPGGNO\nJa+lbmPMp/I5xtnMGFNrjDlojDlujPmBMeb3k3WuIx+4xvpwDfmEMabMGPPPxphjyTX6QrK+1Bjz\nveQ19NXkjfqQB9dYoxeNMQMp11FTvsc6mxlj5hhj/sUYszt5nPU15NuwZ4yZI+n/SLpH0nJJ9xlj\nlud3VJjGHdbaJj6bxTdelHT3FbWtkg5YaxskHUgeI39e1PvXSJK2Ja+lJmvtN2d4TLjsgqT/aq29\nRdIqSZ9N/vvDdeQP062PxDXkF5OSPmmtvVVSk6S7jTGrJP2JptaoQdKwpP+YxzHOdtOtkST9fNvT\nDgAABjBJREFUQcp11J2/IULS70s6nnKc9TXk27An6eOSwtbaH1lrz0v6O0nr8zwmwPestYckvXtF\neb2kl5JfvyRpw4wOCmmmWSP4hLX2tLW2K/n1qKb+ob1ZXEe+cI31gU/YKWPJw2Dyf1bSJyXtTNa5\nhvLoGmsEnzDG1Ej6tKQXksdGOVxDfg57N0uKpBwPiv8z9yMraZ8x5vvGmAfzPRhMq8pae1qaeqEk\n6YN5Hg+u7veMMT3JX/PkVwR9wBizRNLPS/qeuI5854r1kbiGfCP562fdkt6WtF/Sv0oasdZeSD6E\n13V5duUaWWt/eh09nbyOthljSvM4xNnueUmPSEokj0PK4Rryc9gzV6nxXxz85xPW2pWa+nXbzxpj\n1uR7QECB+jNJH9LUr9OclvQ/8zscGGPmS2qX9J+ttefyPR6ku8r6cA35iLX2orW2SVKNpn5b65ar\nPWxmR4VUV66RMeZjkj4v6SOS/p2kGyVtyeMQZy1jzDpJb1trv59avspDr3sN+TnsDUqqTTmukfST\nPI0F07DW/iT559uS/lFT/4cO/4kaY6olKfnn23keD65grY0m/+FNSPoLcS3llTEmqKkg8dfW2q8l\ny1xHPnG19eEa8idr7Yik72jq/ZULjTElyb/idZ1PpKzR3clfk7bW2klJfyWuo3z5hKR7jTFvaeqt\nbJ/U1E5f1teQn8Pe65IaknedmSvpNyV9Pc9jQgpjzDxjTMVPv5a0VtIb134W8uTrku5Pfn2/pF15\nHAuu4qchIunfi2spb5Lvi/hLScettc+l/BXXkQ9Mtz5cQ/5hjLnJGLMw+XW5pF/W1HsrD0r69eTD\nuIbyaJo1+mHKf9Aymno/GNdRHlhrP2+trbHWLtFUBvona+1vK4dryFjr3x305G2Tn5c0R9JXrLVP\n53lISGGM+TlN7eZJUomkv2GN8s8Y87eSbpe0SFJU0h9J6pD095LqJJ2U9BlrLTcIyZNp1uh2Tf36\nmZX0lqSHfvr+MMwsY8wvSXpVUq8uv1fiUU29L4zrKM+usT73iWvIF4wxjZq6ecQcTW0s/L219snk\n64a/09SvB/6LpI3JHSTMsGus0T9JuklTvzLYLel3Um7kgjwwxtwu6XPW2nW5XEO+DnsAAAAAgNz4\n+dc4AQAAAAA5IuwBAAAAQBEi7AEAAABAESLsAQAAAEARIuwBAAAAQBEi7AEACp4x5qIxptsY84Yx\n5h+MMR9w0Ot2Y8zu5Nf3GmO2XuOxC40xv5ty/DPGmJ25fm8AANxE2AMAFINxa22TtfZjks5L+p3U\nvzRTsv43z1r7dWvts9d4yEJJv5vy+J9Ya3/9Go8HAGDGEPYAAMXmVUn1xpglxpjjxpg/ldQlqdYY\ns9YYc8QY05XcAZwvScaYu40xPzTGvCbp137ayBjzgDHmfye/rjLG/KMx5ljyf7dJelbSh5K7il9M\nfs83ko8vM8b8lTGm1xjzL8aYO1J6fs0Y821jTL8xpm1mTw8AYLYg7AEAioYxpkTSPZJ6k6VlknZY\na39e0nuSHpP0y9balZI6Jf0XY0yZpL+Q9KuSVktaPE37/yXp/1lrb5W0UtIPJG2V9K/JXcU/uOLx\nn5Uka+0KSfdJein5vSSpSdJvSFoh6TeMMbXOZg4AwPsR9gAAxaDcGNOtqQB3UtJfJus/ttYeTX69\nStJySd9NPvZ+ST8r6SOSBqy1/dZaK+mVab7HJyX9mSRZay9aa89eZ0y/JOnl5ON/KOnHkj6c/LsD\n1tqz1toJSX3JcQAA4KqSfA8AAAAXjFtrm1ILxhhpajfvUknSfmvtfVc8rkmS9WBM5hp/N5ny9UXx\n7zEAwAPs7AEAZoujkj5hjKmXJGPMB4wxH5b0Q0lLjTEfSj7uvmmef0DSf0o+d44xZoGkUUkV0zz+\nkKTfTj7+w5LqJJ1wYyIAAGSCsAcAmBWste9IekDS3xpjejQV/j6S/FXKByXtSd6g5cfTtPh9SXcY\nY3olfV/SR621MU39WugbxpgvXvH4P5U0J/n4r0p6wFo7KQAAZoiZensCAAAAAKCYsLMHAAAAAEWI\nsAcAAAAARYiwBwAAAABFiLAHAAAAAEWIsAcAAAAARYiwBwAAAABFiLAHAAAAAEXo/wNG9dfYRvLV\n1wAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x11c7ac5f8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"val.loc[:, 'Prediction'] = np.round(p)\n",
"val.plot.scatter(x='Prediction', y='Sales', figsize=(15,10), title='Prediction vs Sales', \n",
" ylim=(0,40), xlim=(0,40))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"Our model is doing reasonably well for items that have sales under about 25 units. \n",
"\n",
"But we can see a \"cluster\" of products that have about 30 sales a day where our model underpredicts. We should definitely investigate this group, understand why this is happening, and see if we can fix it.\n",
"\n",
"Are we missing information? Maybe these are products bought in bulk? Or the units are different (pounds vs ounces)?\n",
"\n",
"Now, to get predictions for each time serie, your just need put the product code, week and compute the new features.\n",
"\n",
"Remember our baseline at RMSLE 0.51581? We improved it to 0.4063, which means a **21% error reduction!**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Improving the Solution\n",
"\n",
"Some ideas you can try to improve this solution:\n",
"\n",
"- Investigate products with sales over 25 units.\n",
"- Try more features (lags, differences) not only at the product level, but at the global.\n",
"- Try using LightGBM native support for categorical features with Product Code.\n",
"- Try tuning the model hyperparameters.\n",
"- Try a neural network and ensemble with GBM.\n",
"\n",
"Remember to let us know how your experiments went by leaving a comment!"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [Root]",
"language": "python",
"name": "Python [Root]"
},
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@sahikatia
Copy link

thank you very much for your clear explanation on this subject. I will test it and see if I can add a more

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment