652 lines
38 KiB
Plaintext
652 lines
38 KiB
Plaintext
{
|
|
"nbformat": 4,
|
|
"nbformat_minor": 0,
|
|
"metadata": {
|
|
"colab": {
|
|
"name": "Lab4_Exercises",
|
|
"provenance": []
|
|
},
|
|
"kernelspec": {
|
|
"name": "python3",
|
|
"display_name": "Python 3"
|
|
},
|
|
"language_info": {
|
|
"name": "python"
|
|
}
|
|
},
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**Last week, we learned how to fit a simple model to a single or multiple predictors to predict a given target variable. However, by fixing the number of predictors and by fixing the model's hyper-parameters beforehand, we restricted ourselves to a single model. By doing so, we omitted to explore a broader range of models, one of which might better explain to relationship between our input and target variables.**\n",
|
|
"\n",
|
|
"**In this lab, we'll experiment with the general methodology of model selection, meaning that we'll define a set of predifined models, and we'll retain the one that minimizes the out-of-sample error.**"
|
|
],
|
|
"metadata": {
|
|
"id": "H3g2zh1W2t_c",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**First, import the necessary libraries**"
|
|
],
|
|
"metadata": {
|
|
"id": "1U5VE1l3YtKA",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {
|
|
"id": "zkt8eCSmYkbz",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import pandas as pd\n",
|
|
"from sklearn.linear_model import LinearRegression\n",
|
|
"from sklearn.impute import SimpleImputer\n",
|
|
"from sklearn.model_selection import train_test_split, cross_validate, LeaveOneOut, GridSearchCV\n",
|
|
"from sklearn.metrics import mean_squared_error, r2_score, make_scorer\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from sklearn.preprocessing import PolynomialFeatures\n",
|
|
"from sklearn.pipeline import make_pipeline\n",
|
|
"from sklearn.feature_selection import SequentialFeatureSelector"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**1) In this lab, we will be working with the 'fish_lab.csv' dataset, which contains several characteristics about fish, such as their weights, lengths, and species, among others. Load the dataset as a pandas dataframe, inspect its properties and check for any missing values. Change the data type of 'Species' to 'category'.**"
|
|
],
|
|
"metadata": {
|
|
"id": "I_O6yz5QY8bH",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"file = '../data/fish_lab.csv'\n",
|
|
"\n",
|
|
"##Read dataframe##\n",
|
|
"\n",
|
|
"df = pd.read_csv(file,index_col=0)\n",
|
|
"df= df.astype({'Species':'category'})\n",
|
|
"print(df.head())\n",
|
|
"print(df.info())"
|
|
],
|
|
"metadata": {
|
|
"id": "TO0Fesd7Zp0C",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 2,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
" Species Weight Length1 Length2 Length3 Height Width\n",
|
|
"0 Bream 242.0 23.2 25.4 30.0 11.5200 4.0200\n",
|
|
"1 Bream 290.0 24.0 26.3 31.2 12.4800 4.3056\n",
|
|
"2 Bream 340.0 23.9 26.5 31.1 12.3778 4.6961\n",
|
|
"3 Bream 363.0 26.3 29.0 33.5 12.7300 4.4555\n",
|
|
"4 Bream 430.0 26.5 29.0 34.0 12.4440 5.1340\n",
|
|
"<class 'pandas.core.frame.DataFrame'>\n",
|
|
"Int64Index: 159 entries, 0 to 158\n",
|
|
"Data columns (total 7 columns):\n",
|
|
" # Column Non-Null Count Dtype \n",
|
|
"--- ------ -------------- ----- \n",
|
|
" 0 Species 154 non-null category\n",
|
|
" 1 Weight 151 non-null float64 \n",
|
|
" 2 Length1 152 non-null float64 \n",
|
|
" 3 Length2 147 non-null float64 \n",
|
|
" 4 Length3 156 non-null float64 \n",
|
|
" 5 Height 159 non-null float64 \n",
|
|
" 6 Width 154 non-null float64 \n",
|
|
"dtypes: category(1), float64(6)\n",
|
|
"memory usage: 9.2 KB\n",
|
|
"None\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**2) Does the dataframe contain any missing values ? If yes, replace the missing values by the sample mean for continuous variables. For missing categorical variables, replace them by the most frequent occurence of the corresponding column. Check the 'SimpleImputer()' class of the scikit.learn library.**"
|
|
],
|
|
"metadata": {
|
|
"id": "ARI4jGJDZ24R",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"imp_mean = SimpleImputer(missing_values=np.nan, strategy='mean')\n",
|
|
"cont = df[['Weight','Length1','Length2','Length3','Width']]\n",
|
|
"imp_mean.fit(cont)\n",
|
|
"df[['Weight','Length1','Length2','Length3','Width']] = imp_mean.transform(cont)\n",
|
|
"\n",
|
|
"imp_cat = SimpleImputer(missing_values=np.nan, strategy='most_frequent')\n",
|
|
"cat = df[['Species']]\n",
|
|
"imp_cat.fit(cat)\n",
|
|
"df[['Species']] = imp_cat.transform(cat)\n",
|
|
"print(df.isna().sum())"
|
|
],
|
|
"metadata": {
|
|
"id": "X6p-MsWdagTK",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 3,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Species 0\n",
|
|
"Weight 0\n",
|
|
"Length1 0\n",
|
|
"Length2 0\n",
|
|
"Length3 0\n",
|
|
"Height 0\n",
|
|
"Width 0\n",
|
|
"dtype: int64\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**3) Unfortunately, categorical variables cannot be fed directly to the model. One approach to deal with the problem is to create 'dummy variables' out of the categorical one, where each dummy variable takes the value '1' if the observation belongs to that class, and '0' otherwise. Check the method 'get_dummies' of the pandas library to generate such dummy variables.** "
|
|
],
|
|
"metadata": {
|
|
"id": "jmndzfo-anRy",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"df = pd.get_dummies(df, columns=['Species'], drop_first=False)\n",
|
|
"\n",
|
|
"print(df)"
|
|
],
|
|
"metadata": {
|
|
"id": "da0wVt_ibvcW",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 4,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
" Weight Length1 Length2 Length3 Height Width Species_Bream \\\n",
|
|
"0 242.0 23.2 25.4 30.0 11.5200 4.0200 1 \n",
|
|
"1 290.0 24.0 26.3 31.2 12.4800 4.3056 1 \n",
|
|
"2 340.0 23.9 26.5 31.1 12.3778 4.6961 1 \n",
|
|
"3 363.0 26.3 29.0 33.5 12.7300 4.4555 1 \n",
|
|
"4 430.0 26.5 29.0 34.0 12.4440 5.1340 1 \n",
|
|
".. ... ... ... ... ... ... ... \n",
|
|
"154 12.2 11.5 12.2 13.4 2.0904 1.3936 0 \n",
|
|
"155 13.4 11.7 12.4 13.5 2.4300 1.2690 0 \n",
|
|
"156 12.2 12.1 13.0 13.8 2.2770 1.2558 0 \n",
|
|
"157 19.7 13.2 14.3 15.2 2.8728 2.0672 0 \n",
|
|
"158 19.9 13.8 15.0 16.2 2.9322 1.8792 0 \n",
|
|
"\n",
|
|
" Species_Parkki Species_Perch Species_Pike Species_Roach \\\n",
|
|
"0 0 0 0 0 \n",
|
|
"1 0 0 0 0 \n",
|
|
"2 0 0 0 0 \n",
|
|
"3 0 0 0 0 \n",
|
|
"4 0 0 0 0 \n",
|
|
".. ... ... ... ... \n",
|
|
"154 0 0 0 0 \n",
|
|
"155 0 0 0 0 \n",
|
|
"156 0 0 0 0 \n",
|
|
"157 0 0 0 0 \n",
|
|
"158 0 0 0 0 \n",
|
|
"\n",
|
|
" Species_Smelt Species_Whitefish \n",
|
|
"0 0 0 \n",
|
|
"1 0 0 \n",
|
|
"2 0 0 \n",
|
|
"3 0 0 \n",
|
|
"4 0 0 \n",
|
|
".. ... ... \n",
|
|
"154 1 0 \n",
|
|
"155 1 0 \n",
|
|
"156 1 0 \n",
|
|
"157 1 0 \n",
|
|
"158 1 0 \n",
|
|
"\n",
|
|
"[159 rows x 13 columns]\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**4) Now let's make some predictions. Select the variable 'Weight' as predictor and 'Height' as the target variable. Generate a scatter plot of the two variables, and fit a simple linear regression model to the data. Split your datasets into a training and test set following a 60/40 partition. Evaluate the model on the Mean Squared Error and on the coefficent of determination. Report such metrics for both the training and test sets.**"
|
|
],
|
|
"metadata": {
|
|
"id": "SNR5NKrlbxbK",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"x = df['Weight'].values[:, np.newaxis]\n",
|
|
"y = df['Height'].values\n",
|
|
"model = LinearRegression(fit_intercept=True)\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"ax.scatter(x.ravel(),y)\n",
|
|
"ax.set_xlabel('Weight')\n",
|
|
"ax.set_ylabel('Height')\n",
|
|
"\n",
|
|
"\n",
|
|
"X_train, X_test, y_train, y_test = train_test_split(x, y, train_size=0.6, test_size=0.4, shuffle=True)\n",
|
|
"\n",
|
|
"model.fit(X_train,y_train)\n",
|
|
"print('Model\\'s coefficients : {}'.format(model.coef_))\n",
|
|
"print('Model\\'s intercept : {}'.format(model.intercept_))\n",
|
|
"\n",
|
|
"#Make predictions for both the training and the test sets##\n",
|
|
"\n",
|
|
"y_pred_train = model.predict(X_train)\n",
|
|
"y_pred_test = model.predict(X_test)\n",
|
|
"\n",
|
|
"##Compute th coefficient of determination and the mean square error on both sets##\n",
|
|
"\n",
|
|
"MSE_train = mean_squared_error(y_train, y_pred_train)\n",
|
|
"MSE_test = mean_squared_error(y_test, y_pred_test)\n",
|
|
"R2_train = model.score(X_train, y_train)\n",
|
|
"R2_test = model.score(X_test, y_test)\n",
|
|
"\n",
|
|
"print('MSE on training set : {}'.format(MSE_train))\n",
|
|
"print('R2 on training set : {}'.format(R2_train))\n",
|
|
"print('MSE on test set : {}'.format(MSE_train))\n",
|
|
"print('R2 on test set : {}'.format(R2_test))\n",
|
|
"\n",
|
|
"xfit = np.linspace(0,2000)\n",
|
|
"xfit = xfit[:, np.newaxis]\n",
|
|
"yfit = model.predict(xfit)\n",
|
|
"\n",
|
|
"ax.plot(xfit, yfit, label='Regression line', color='red')\n",
|
|
"\n",
|
|
"plt.show()"
|
|
],
|
|
"metadata": {
|
|
"id": "ZAQqipascyCd",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 5,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Model's coefficients : [0.00919392]\n",
|
|
"Model's intercept : 5.293929818146502\n",
|
|
"MSE on training set : 8.420729290375146\n",
|
|
"R2 on training set : 0.5162377033178694\n",
|
|
"MSE on test set : 8.420729290375146\n",
|
|
"R2 on test set : 0.39421490361091693\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"text/plain": "<Figure size 432x288 with 1 Axes>",
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAu1UlEQVR4nO3deXiTVfbA8e9pqVhEKCoqVBEFBwUX0LqiI+IIDgiiOO7KOI7oz2XEYSq4giso474z4uCCigtWkFFEUVERsFCQXVRELQiolEUKlPb+/rgJhjRvljbJm7zv+TxPn6Zvs5ym7cnNfc89V4wxKKWU8o8ctwNQSimVXpr4lVLKZzTxK6WUz2jiV0opn9HEr5RSPtPA7QDisccee5jWrVu7HYZSSmWVWbNm/WyMaR5+PCsSf+vWrSktLXU7DKWUyioisjzScZ3qUUopn9HEr5RSPqOJXymlfEYTv1JK+YwmfqWU8hlN/Eop5TOa+JVSymc08SulVCZauRIGDICKiqTftSZ+pZTKJBs2wG23Qdu28MQT8MknSX8ITfxKKZUJtm6Fxx6DNm3gzjuhVy9YtMh+TjJN/Eop5SZj4LXXoH17uPZa6NABZs6EV16xLwIpoIlfKaXc8vHHcOyxcM45kJ8PEyfClClw1FEpfVhN/EoplW4LFtgpnC5dYMUKePZZmDMHevQAkZQ/vCZ+pZRKlx9/hMsug8MOsydthw+Hr76CSy+F3Ny0hZEVbZmVUiqrrVsH994LDz0E1dVw3XVw882w++6uhKOJXymlUmXLFnjySbjrLvjlF7jgArj7bnB5Yymd6lFKqWSrqYGXX4aDD4brr4dOnWDWLBgzxvWkD5r4lVIquaZMgaOPtqP7Jk1g0iSYPBmOOMLtyLbTxK+UUskwdy6cdhqccgqsWQMvvACzZ0O3bm5HVosmfqWUqo/vv4d+/ex0zsyZ8O9/w5IlcNFFkJOZKVZP7iqlVF2sXQvDhsEjj9ivi4th8GBo1szduOKgiV8ppRKxebPtqXPPPbZz5iWXwB13QKtWbkcWt8x8H6KUUpmmpsbO27drZ0f3xxxjV9uOHp1VSR808SulVGzvvWerci65BJo3hw8+gHfesStws5AmfqWUcjJ7Npx6KnTvDuvX29r8mTOha1e3I6sXTfxKKRVu2TK48EI48kgoK7OtFhYtgvPOy9hKnUToyV2llAr65RfbUuHxx23TtJtughtugKZN3Y4sqTTxK6VUZSU8/LDtlrlhg+2WefvtUFjodmQpoYlfKeVf1dXw3HN2j9vyctsjf9gwuwuWh2X/ZJVSSiXKGLvbVceOtj9+YaHdDWv8eM8nfdDEr5Tym5kz4eST4fTT7WKs116D6dPhj390O7K00cSvlPKHr7+Gc8+1C68WLrSrbxcuhLPPTst2h5lE5/iVUt62ejXceSc89RTstBPceqtdebvrrm5H5hpN/Eopb/rtN3jwQbvlYWUl/P3vMGQItGjhdmSu08SvlPKWbdvg2Wdh6FBYuRLOPNNW6rRr53ZkGUMTv1LKG4yxVTmDB8PixXD88fbEbefObkeWcfTkrlIq+33+OZx4IvTpY18A3nwTPv1Uk74DTfxKqez11VfQt68d3X/zjT2BO3++fQHwWaVOIlKW+EVkXxH5UEQWisgCEbkucHw3EZksIksDnzN/uxqlVGb56Sf4v/+D9u1ty+Q77rDlmldcAQ10BjuWVI74twEDjTHtgWOBq0WkPTAY+MAYcyDwQeBrpZSKbcMGe9K2bVt45hm48ko70r/1VthlF7ejyxopS/zGmJXGmNmByxuARUAhcAbwXOBqzwF9UhWDUsojqqrgiSdswr/9dujR4/dFWHvu6XZ0WSct74lEpDXQCZgB7GWMWRn41k/AXg636Q/0B2iVZduaKaWSxBgYNw5uvBGWLoWTToIJE+Doo92OLKul/OSuiDQG3gAGGGPWh37PGGMAE+l2xpiRxpgiY0xR8+bNUx2mUirTfPKJPWl79tl2xe3bb8OHH2rST4KUJn4RycMm/THGmHGBw6tEpEXg+y2A1amMQSmVZRYuhN69bdO0H36AUaNg7lzo2VMrdZIklVU9AowCFhljHgj51nigX+ByP+CtVMWglMoiK1bA5ZfDoYfaFsn33GPLNf/2N7sblkqaVM7xdwYuBuaJyJzAsZuA4cCrInIZsBw4J4UxKKUy3fr1cN998MADtt3CP/4BN98Me+zhdmSelbLEb4z5FHB6X3ZKqh5XKZUltm61C67uvBN+/hnOPx/uugsOOMDtyDxPV+4qpdKrpgbGjoWDD4brroPDD4fSUnjpJU36aaJL3JSqg5KyckZMWsKKikpaFuRT3L0dfTp5c2PupPrwQ7jhBpvoDzsM3n0XunXTk7ZppiN+pRJUUlbOjePmUV5RiQHKKyq5cdw8SsrK3Q4tc82bZxddde0Kq1bZDc5nz4bu3TXpu0ATv0eUlJXTefgU9h88kc7Dp2gSSqERk5ZQWVW9w7HKqmpGTFriUkQZ7Icf4NJL7XTO55/DiBG2UueSS7RSx0U61eMBwRFoMBkFR6CATj8kUXB6p7yiMuL3Vzgc96WKChg+HB5+2K6+HTjQrr7dbTe3I1PoiN8TdASaeqHTO05aFuSnMaIMtWWLLcts08aWaP7lL7BkiR3pa9LPGJr4PcBppKkj0OSJ9OIaKj8vl+LuPt7ar6YGXnzRbm84cCAcdRSUlcHzz8N++7kdnQqjid8DnEaaOgJNnmgvooUF+Qw761D/TqtNngxFRXDxxbD77vD++7Za5/DD3Y5MOdDE7wHF3duRn7fjiTLfj0CTzOlFtLAgn88Gd/Vn0i8rs6WY3brB2rUwZgx88QWcouszM52e3PWAYNLRuvK6c6rLDz2hK+zYSta3L67ffQe33GIT/W672Tn9q66Chg3djkzFSWxn5MxWVFRkSktL3Q5DeVR4VRTYpN73yELemFW+w/Fg8i/044vrL7/YxmmPPQY5OTBgAAwaBAUFbkemHIjILGNMUfhxHfEr33Oqinp5xg9Uhw2Mgkn/s8Fd0xihyyor4dFHbdLfsAH69bN73O6zj9uRqTrSxK98z+nEbXjSj3V9z6muhhdesPvZ/vijXXl7771wyCFuR6bqSU/uKt9zOnHr1EjA89VSxsA770CnTnbVbYsWtsfOxIma9D1CE7/yvUhVUXk5Qk5O7dSflyPePqFbWmqrcnr0gE2bbBfNGTOgSxe3I1NJpIlf+V6fToUMO+tQCgvyEewcfuOdG1BdU3uqp/HODbx5Qvfbb20//KOOsg3VHnnEboF4zjnaRM2DdI5fKWzyD03o+w+eGPF6FZuq0hVSeqxZYzc/efJJyMuzZZrFxdCkiduRqRTSxK9UBC0L8iP25fHM/P6mTfDgg/Zk7aZNcNllMHSonc9XnqdTPUpF4NnV0Nu2wTPPwIEH2tH9KafA/Pnw9NOa9H1ER/zK0+q6U5bnVkMbAxMm2NbICxfCccfZE7cnnOB2ZMoFmviVZ9V3n4Lwef+sNWOGnbf/5BP4wx9g3Djo00dP2vqYTvUoz/L9PgVLl9p++Mcea3e9evJJO61z5pma9H1OR/zKs3y7T8GqVbalwsiRtnHa0KG2R37jxm5HpjKEJn6VMeo6H+/E85U54TZutJ0yR4yAzZuhf3+47TbYay+3I1MZRqd6VEYI3drQYOfjB4ydw8G3vlPnjeM9W5kTrqoKnnoK2raFIUPgtNPsCdzHH9ekryLSEb9yXUlZOQNfnRuxKVplVQ3Fr80FEt843nOVOeGMgZISGDzYzuGfeKL9+thj3Y5MZThN/MpVwZG+UydMgKoaw4hJS+qUsD1TmRPus8/ghhtg2jQ4+GAYPx5OP11P2qq4aOJXrgjd2Soenj8hG6/Fi20tfkkJtGxpF2P16wcN9F9ZxU//WlTaRdrxKhbPnpCN18qVtjpn1Cho1AjuvtvugNWokduRqSykiV+lXaT6+mg83wo5mvXrbZXOAw/Yk7hXX21bLTRv7nZkKotp4ldpl8i0TUF+HkN7d/DmPH00W7fa/jl33mk7aJ53nu2i2aaN25EpD9DEr5Im3jp8p/p68Okm5qGMgddeg5tugm++gZNPhvvug6Ja+2UrVWdiolRTZIqioiJTWlrqdhgZLdmLnxJ97KHjF1BRWbtXfbNGeQzpteOIPdIcf35eLsPOOtS/CR/go49spc4XX9gtDu+7z9bka6WOqiMRmWWMqTVq0BG/B9S3GVkyHzvc2k1VtWLxfH19oubPt7X4EyfCPvvAf/8LF18Mubmxb6tUHWji94BozchSnUzjOVEbKRbP1tcn4scfbUuF556DXXe1m6Jcey3k+7yCSaWcJn4PcKMZmdbh10NFhU3yDz0ENTW2LPOmm2D33V0OTPlFynr1iMizIrJaROaHHBsqIuUiMifw0SNVj+8nTjXuqap9D+2rEy/f1+EDbNlik32bNjB8OPTtC0uWwP33a9JXaZXKJm2jgdMiHH/QGNMx8PG/FD6+b6S7GVmidfixYikpK6fz8CnsP3ginYdPqXNTtoxVUwMvvQQHHQTXXw9HHgmzZ8OLL0Lr1m5Hp3woZVM9xpipItI6Vfevfpfuk6Wxpm1ygKaN8qjYVBUzlpKycopfn0tVta0uK6+opPj1ujVly0jvvw+DBtlE37EjvPcenHqq21Epn3Njjv8aEbkEKAUGGmPWRrqSiPQH+gO0atUqjeFlp3SeLI1Whw9QAzTaqQFlt3WLeV+3T1iwPekHVVUbbp+wILsT/9y5NuFPmgT77QcvvAAXXAA52glduS/df4VPAm2AjsBK4H6nKxpjRhpjiowxRc11eXpGiTS1FC7ek7lrN9Wu/Y92POMtXw6XXAKdOsHMmXb+fvFiuOgiTfoqY6R1xG+MWRW8LCL/Ad5O5+Or5AidWnIa+YeezHVzcVnarF0L99wDjz5qv77hBlubX1DgalhKRZLWxC8iLYwxKwNfngnMj3Z9lT6JJufg1JLTKtzgydxYi8sK8vMirvgtyM9L5o+XOps322R/zz2wbp1tkXzHHbDvvm5HppSjlCV+EXkZ6ALsISI/AkOALiLSETDAd8AVqXp8FV1oom+an8eGLduorkn8BGusE8uxFpcN7d2B4tfmUlXz+zx/Xo4wtHeHpP2sKVFdDWPGwK23wvffQ48etkTz0EPdjkypmFJZ1XN+hMOjUvV4Kn7ho/BII+5ETrBGO7Eca3FZ1rVvMMaesB00CL780jZPGz3aNlNTKkvoyl0fSXS1bTJOsDpVAIWeA8ia9g2zZ9u5+w8+gAMOgFdegb/8RU/aqqwT11+siHwQzzGVOvVd5FSX1bbJkO7FZSmxbJktxTzySFum+fDDsGgRnHuuJn2VlaKO+EVkZ6ARdp6+GRDsD9sEyIIhmjcko/tmoqttkyXrpnJC/fyz3eLw8cftnrY332xH/E2auB2ZUvUSa6rnCmAA0BKYxe+Jfz3wWOrCUqGS0X2zLk3SklVZkzVTOUGVlXZUP2wYbNwIf/ub3e+2MIt+BqWiiJr4jTEPAw+LyLXGmEfTFJMKk4zum9FW2zbKy2FTVU2t46cf3iLu+/eE6mrbIvm226C8HHr1spU67du7HZlSSRXXBKUx5lEROV5ELhCRS4IfqQ5OWcnovnnyQZFXP190bCua7dIw4vc+XLwm7vvPasbA22/D4YfDZZfZzVCmToXx4zXpK0+K9+TuC8C/gROAowIfuglomiTjBKlTEv9w8RpX+vlnjJkzoUsXO7rfuhVefx0+/xxOPNHtyJRKmXjLOYuA9iYbNuj1oGScII2W3OMpufScpUvt5ievvw577mlP4F5+OeRlyYphpeoh3sQ/H9gb21hNuaC+J0ijJffi7u2itl3wlNWrbUuFp5+Ghg3tfP6//mW3PlTKJ2KVc07AtlfYFVgoIjOBLcHvG2N6pzY8lSzRkntWl1zG67ff4IEH4L77bNXO5ZfDkCGw995uR6ZU2sUa8f87LVGolIuV3LOu5DJe27bBqFG2HPOnn+DMM22ZZjsPvptRKk6xyjk/TlcgKvU8m9wjMQbeegtuvNH2w+/cGd54A44/3u3IAJ+0qlYZK645fhHZgJ3yCbWO33fR+jbZgSlnkZIGeHyqJhHTpkFxsf3crh28+SaccQaIxL5tGiRjJbZS9RHvyd2HgB+Bl7Crd8/D7qQ1G3gW235ZpUGkpFH82lwQdti31peJZMkSO8J/8007d//003bVbYPk9yKsz4g9GSuxlaqPeP8jehtjDg/5eqSIzDHGDBKRm1IRmN/Em0giJY3QXvZBwUQSvI2n3wn89JOdw3/mGcjPhzvvhOuvh112ScnDlZSV77CHwPYXX+J7ofX1ugmVEeJtLbhJRM4RkZzAxznA5sD3tLa/nkI7Zxp+H7GHduAMdudMpLtm8H6i3W9W27DBVua0bQujRvHt2RfT87rR7L+xE50fnZGyn3Po+AW1XmyragwDxs6Jq3NqU4ceSE7HlUq2eBP/hcDFwGpgVeDyRSKSD1yToth8I9pbf6h7S+Vckaj3m7WqquCJJ2zCv+MO6NmTya9/SM8Dz2XBtp1T/iIXaeOaoHge1+lUQ4acglA+EG+vnm+NMb2MMXsYY5oHLn9tjKk0xnya6iC9LtZb/1gtlfNyhLzcHbNGfl4u1Q4LrbN2SsEYeO012z/n6qvh4INhxgwYO5ahi7ZmzItcrMetcNjgxum4UskWawHXDcaY+0TkUSJM6Rhj/pGyyHwkVsuEaIm6MEpVj9NuW1nZimHqVNsLf8YM6NABJkyAnj23D5Mzbd482uP6skWGyiixRvyLAp9Lsf34wz9UEsRqwuaUEAoL8vlscFfHE4qe2P1qwQLbQO2kk+DHH+1irLlz4fTTd5gbcXqOmubn1Wvnskji2acgWhL3xO9FZbWoid8YMyHw+TljzHPAa8HLga9VEvTpVMiwsw6lsCAfwSb0YWcduj2hx0oUTieHgaj3m9HKy+Hvf4fDDrOj/WHD4KuvbHlmbm6tq0d6jvJyhN+2bkv6ye2hvTuQl+M8IR8ricf6fSuVahJPw00ROQ4YBTQ2xrQSkcOBK4wxV6U6QICioiJTWlqajofKWNHKPZ2qfYLvCLLKunW2n86DD9qNUa6+2m55uPvuEa8e+rw0zc9DxM6VtyzIZ9PWbRE3jA8+L/WpxY/2uJ4smVVZSURmGWNqtdBPZAFXd2A8gDFmroj8MXnhqViitVvItPntOtm6FZ58Eu66y+51e8EF9vL++zveJHwxW0VlFfl5uTx4bkf6dCpk/8ETI95uRUVlvVfP+qr9hfKcuJc0GmN+kB3rzdK/c3eWS1V/lqw+WVhTA6++anvjL1sGp5xiR/xHHBHzprHKYHNEIlY2tSzI19WzytfireP/QUSOB4yI5InIv/j9xK+KQzyLtILXS/RkZNaeLJwyBY4+Gs4/H5o0gXffhcmT40r64PyOJvjcRkr6wefFE++SlKqjeBP/lcDVQCFQDnQMfK3iFGt0CvG/OITLupOFX34JPXrY0f2aNfD88zB7NnTvntAqJqd3NJEWrgWPB5+XZOxjrFS2imuqxxjzM3b1rqqjeEaY9Zl+yOQ55+AUF98v5+YZr/DnsslI06Z2Sufaa2Hnnet0v5E2lxFwXLhWY8wOlVLFr8/d3tgOIC9XMv9dklJJEGsBV8SFW0G6gCt+TvPwwTrzFYFRfiTZPP1QUlbOsDHTuPSTV7h01gQAnj3mLPYeNpSeXQ6p132Hbi5TXlGJEL1xVK3RfPiVteuU8olYI/7QGsrbgSEpjMXTIo1Og3Xm0Xq/QBZPP2zezA+33MmkKWNosvk33jzkZO4/8SJWNNmTwumr6dml/g8RfKcTq4Fd+DmPEZOWRGy0pid3lR/E2oFr+yItERmgi7bqLtLWh0515qGy4iRtuJoaGDMGbrmFa7//no/2P5J7u/Rj0Z4HbL9KeUXl9nc6yaiDj6etReh96sld79DdzBKXyA4V+ka4nsLn4Z3qzMHOVWfCH3HC/1TvvQeDBsGcOXDEEfzj1GsZv8fBta4msH2EHvqOp66byDhNpTktYsvqEtgo/JYEdTezuom3qkelQLQePMuG94zahycdEqoyKiuDU0+1lTkVFfDSS/DFF3S9+vxapaax5uLr0lUz0ZLWrC2BjaKuVWHZLJ5qOVVb1MQvIhtEZL2IrAcOC14OHk9TjJ6V6cknrn+qZcvgwgtt7X1ZGTz0kN3c/PzzIScnYqlpPG8dE917INGS1qwrgY2DH5OgTtnVTaw5/l3TFYiXxPt2O9K8fya9NY/6T/XLL3D33fD445CTY/e6HTQImjatdf3wKa54dhLLrcOuJImWtGZyCWxd+DEJenXKLtWSvwu1zyU655jJySfSP1XDqi1cOmsCG/Y5l8ZbK5G//hVuvx322Sfu+y3u3o4BY+dEvY5TLb5y5sckGKlaLpPeNWcqneNPMi+93Q6disqpqeYvX07mw/9cweCPRzOjsD1n/P1xSq65I6GkD/bFrlmj6D3tCz2crFIl06cOU8GLU3bpkLIRv4g8C5wOrDbGHBI4thswFmgNfAecY4xZm6oY3JDtb7fDp6n6HtGSzW9N5O8Tn+agn5czp8WBXH/6QGa0OhSgznXvQ3p1qDVSC/J6skqVTJ86TJVMftecqVI51TMaeAx4PuTYYOADY8xwERkc+HpQCmNIu2x+ux0+TbX7orn0emo0xyz/kmXNWnDVGYP5X7vOO/TTqesLWviq29xAJ81INfcqfpoEVTxSlviNMVNFpHXY4TOALoHLzwEf4bHEn81zjsFpqlZrV1I89Xl6Lf6Enxs15YFe1/LW0T1ZvnFbrdvU5wVNk5RS7kj3yd29jDErA5d/AvZyuqKI9Af6A7Rq1SoNocUnVsVONr/d3rziJ4ZMe4ULy95hW24uDx9/Hv85+ix+a9iIB3sekrUvaEqpHblW1WOMMSISrQHcSGAk2K0X0xZYFPFW7GTdSPa33+Chh5g68h4abt3Mq4d148ETLmBN490Ae8Ism1/QVOaLNaDy24rkVEt34l8lIi2MMStFpAWwOs2PXy+ZvGtTnf4xtm2D//4XhgyBlStZ36U7lx3Ul4VNW26/SuioPute0FRWiDWg0rYMyZfucs7xQL/A5X7AW2l+/HqJtuNTIjtmJVvCS/WNgfHj4bDDoH9/aN0aPv2UFh++S//+PbQ0TqVVrBJoL5VIZ4pUlnO+jD2Ru4eI/Iht6TwceFVELgOWA+ek6vFTwaliB9gh4UJ6RyIJvROZPh2Ki+HTT6FdOxg3Dvr02V6po6N6lW6xSqCzvUQ6E6VsxG+MOd8Y08IYk2eM2ccYM8oY84sx5hRjzIHGmD8ZY35N1eOnQqQFMuHcGInE9Y/x1Vdw9tlw3HHw9dfw1FMwfz6ceWZC2x0qlWyxtsHUbTKTT1fuJiB8laCTdI9Eov5jrFoFV10F7dvDpEm2vcLSpXDFFdBAO3Yo98VacezHFcmppv/5CQpfeBRJukcikdYO7G62MmrZB9CmN2zZAldeCbfdBnvumdbYsplWkqRHrIoxrShLPjFZ0AyrqKjIlJaWxr5iGoRXGITLz8ut0wnR0CRTlx2pgrdf/csG+i/9kGs/fYmdf1ljp3fuuQcOPDChePwu0u+5rr9bpdwiIrOMMUXhx3XEn6BIJ1KD6tpuIDzJ1GVHqj4dW9Jn2QzbHvmrr+DEE2HEBDjmmIRiUVYml+4qVV86x58gp/l7gTrtmFVSVs7AV+c6vphAHCeMP/0UOneGvn3tvP348fDxx5r060ErSZSXaeJPUDIrDIIj/Xh6z0dMOIsW2VLME0+k8utvGXbWQNr2HEbnBbtQMmdFwvGo32klifIyneqJIfwE38kHNeeNWeVJ6VkzdPyCqCP9UDsknBUrYOhQGDUKGjdm4dWDuLDxsazF9rjXlY31l83N9pSKRUf8UURaEfvGrHL6HllY79WtJWXlO8zlR7M94axfD7fcAm3bwujRcO218M03XL5Pt+1JP0hXNtaPbvChvExH/FE4neD7cPEaPhvctd73HQ8BzjlsT/p8Og663QE//wznnWf3uz3gAEDno1NFVzErr9LEH0UqE2qszcYBMIaeiz/h8v+8CL+ugK5d4d57oWjH6qxs3vxFKZV+OtUTRSpP8OXGaJNw3PIveev5f/LY+PvYmLsT/O9/8P77tZI+6MpGpVRidMQfRSpO8AVPFjtV8rRb8x2DPhpN129LKd+1OQN7XM/Mzn/mkz+f6niffToVUrr8V16e8QPVxpArQt8jkzNNoatXlfIeTfxRJHupeLRVv3uv/5l/fvoiZ8/7gI0NG3FPl0t57ojTyWnUiGF/bh/zft+YVb79xaTaGN6YVU7RfrvVK0lnWh90fRFSKjm0ZUMadR4+pdZcfJPNG/m/6a9z6azx5Jga3u96Nre078OvOzdGgEY75bJpa3XURBfpfsFWotTnJHSq7rcutIWCUonTlg0pksgoNPSk8E7bqri4bCLXTBtL080bea/jKcz/v4GM+oHtyc0Av22NPdpO1UnoTKoW0hYKSiWPJv4wiSRyp6mQ0uW/8uHiNbXuo2VBPivW/kbvhR/zr09eZN91q5jauhPP9rqS0Y/0587hU6isck6qTokuVVU9mVQtlEkvQkplO038IRKd03YahY6Z/j3BCbTyikquHzuHAWPn0H3FPJ5+7xkOWfUN8/dqw43dr2FG2yMYcfbhQHxJLNJ1op2EjrTyONKLUiSZtHo1k16ElMp2mvhDJDqd4JSow8+aHLzqWwZ/9F/++F0ZPzbZk+tOH8j49idhJAeqDUPHLwCib+0YFCnROZ2EBmq9kL04/fvtt4v1wpZJfdAz6UVIqWyniT9EotMJBY3yWLvJue1C4brVDPzkBfos+Ih1OzfmzpMv44UjTmdrgx3bK1RUVlH82lzOPXrfWn2AQkVLdJFWmXYePiVmL6BY8+SZsno1k16ElMp2mvhDJDKdUFJWzsbN2yLeT9PKDVz9+av0mz0BIzk8fUxfnjz2bNbv3NjxsatqDG/PXcmwsw6t04Yskc5NxDv/vaKiMitKJTPlRUipbOfrcs54O2+GlwwGe+iHL8JquG0r/WZN4OrPX2XXLZt445BTeODEC1nZpHncMX03vGedfo5I0yA75+VEfUcSVJCfx5ZtNbXeHRTk5zG0dwdNtkplKS3nDBPpRG6w82a0k5+Reujn1FRz5oKP+OcnL1K4YQ1TDiji3i5/ZUnz1mn5WZzOTTRskEN+Xm7U6Z78vFxEiHidisoqbe+slAf5NvHXtfPm7RNCeugbw0nLZjP4o/9y8JrvmNPiQAaefj3TWx0W8bYCHN9mNz7/9ldqIrzRatYor/bBODhN6ayrrOLBczvGrOq5fuwcx/vWWnmlvMeziT/WnLVT9Uy0qppbSuZtnzo55KevufGjZ+m8/EuWF+zN1b0HMfGgE8Ch+VrofrwlZeUUvz6Xqurfs39erjCkV4e6/KhRz03EMy8+YtKSqD+31sor5S2eTPzx1OPnikRslObUNfOWknm8OP179q34ieKpz9N70VR+yW/CkD9dwUsdT6Mq13m0Ht7ioK4VKk4vZvUtdYx0+1BaK6+Ut3gy8cdTj+/UHTPS8ZKycv43ZR63TRvLRWX/ozonl0eOO5eRx/RlY8NGMeOJNGJ2Gok7Jfd4XszqWpUTvN7tExbUOhmstfJKeY8nE3889fiFDtMjhSGj25Kycm59aToXzxzPR9NfZ5eqzYw97FQe6nwBq3fdPe544h0xR0vusV7M6lvqGLx9NpR1KqXqx5OJP556/FjTIyWl3zPt1vuZ/OmL7L3xVya3PYZ7T+rH13u0SiiWREbM0ZJ7unrVaK28Ut7nycQfz5y34/RIx5bw9tscftk19Fm9nNkt23HNGYMo3SfxE6/NGuUxpFf8dfDRkrv2qlFKJYsnE3+8c961RrczZkCXC2DqVGjWkiv73Mi7fzjesVLHSa4I959zeMIjZ6cWEAWN8rRXjVIqaTyZ+CHBKYulS+Gmm+D112HPPeGJJ+j3axu+3xB71Wu4+mwO4rSI2pjYL2Y6N6+UipdnE39cVq+GO+6Ap59mS4M8njjhAv5z1JlsWp5PDokn/WCtPtgGaYkm4XWVkR8zeDxaJVAmbZGolMpsOW4H4IqNG23Cb9MGnnqKGX86ixMuG8nDnS9g0052zrwmwbsU2F6rf+O4eZRXVGL4PQmXlJXHvA+n+fpY8/jRTgorpVQ4fyX+bdvg6aehbVsYMgS6d4eFC7mgYz/WNG5Wr7sOJuf6JOHi7u3Iz8vd4Vg88/i6O5VSKhH+SPzGwJtvwiGHwJVXwoEHwrRpdk7/D39wXMwVr7wc2Z6c65OE+3QqZNhZh1JYkI9gp47iOV9Q13cKSil/8v4c/7RpUFxsPx90ELz1FvTqtb1SJ54pmJhCin7qW3ZZlzp6rfhRSiXClRG/iHwnIvNEZI6IJL/RftDAgdC5MyxbBiNHwrx50Lv3DuWZyZgHr6o22++nrtM19VHXdwpKKX9yc8R/sjHm55Q+QteusNtuMGAA7LJLxKskMg8u1N5PN/x+3NoiUFfcKqXi5e2pnp497UcU8WxwDr8vynJqYRw6laNJWCmVydw6uWuA90Rkloj0j3QFEekvIqUiUrpmzZqUBVLcvR2x1uXm5+VuX4nrxlSOUkolk1sj/hOMMeUisicwWUQWG2Omhl7BGDMSGAl2z91kBxC60rVBDlQ5FO6LsMN8uVtTOUoplSyuJH5jTHng82oReRM4Gpga/VbJE77S1Snp5+UKI86u3XNHp3KUUtks7VM9IrKLiOwavAx0A+anM4ZIi6zC5UrkpK+UUtnOjRH/XsCbYksqGwAvGWPeTWcA8VTy1BijSV8p5UlpT/zGmG+Bw9P9uEElZeXkOOy3G0pXvSqlvMofLRsCgnP7sZK+VukopbzMV4nfaW5fsLtl6apXpZQfeHsBV5hoc/tlt3VLYyRKKeUeX434m+bnJXRcKaW8yFeJ32nr3AS31FVKqazmq8QfaSNzgAqH40op5UW+SfwlZeWOPXm0dFMp5Se+SfwjJi2J2FJZQEs3lVK+4pvE71TRY0BLN5VSvuKbxO80nVOo0zxKKZ/xReIvKStn09ZttY7rCl2llB95fgFXSVk5xa/Ppap6xxn+gvw8hvbuoNM8Sinf8fyI//YJC2olfbC1+5r0lVJ+5PnE71S773RcKaW8zvOJXyml1I48n/gb5UX+EQu0P49Syqc8nfhLysrZHGFD3RyBob07uBCRUkq5z9OJ/6ZxXxJpH/WGDXL0xK5Syrc8m/hLysrZFGG0D1DpcFwppfzAs4l/xKQlboeglFIZybOJvzzKbltKKeVnnl+56xclZeWMmLSEFRWVtCzIp7h7Oz2PoZSKSBO/B5SUlXPjuHnbN5Ivr6jkxnHzAF2drJSqzbNTPdFcdGwrt0NIqhGTlmxP+kGVVdV6nkMpFZEvE/9dfQ51O4SkctprwOm4UsrfPJv4nfrse7H/vtNeA7qlpFIqEs8m/uLu7cjPy93hmFf77/vpZ1VK1Z9nT+4GT2r6odLFTz+rUqr+xJhIW5BnlqKiIlNaWup2GEoplVVEZJYxpij8uGenepRSSkWmiV8ppXxGE79SSvmMJn6llPIZTfxKKeUzWVHVIyJrgOV1vPkewM9JDCdZNK7EaFyJ0bgSk6lxQf1i288Y0zz8YFYk/voQkdJI5Uxu07gSo3ElRuNKTKbGBamJTad6lFLKZzTxK6WUz/gh8Y90OwAHGldiNK7EaFyJydS4IAWxeX6OXyml1I78MOJXSikVQhO/Ukr5jKcTv4icJiJLRORrERmcxsfdV0Q+FJGFIrJARK4LHB8qIuUiMifw0SPkNjcG4lwiIt1THN93IjIvEENp4NhuIjJZRJYGPjcLHBcReSQQ25cickSKYmoX8rzMEZH1IjLAjedMRJ4VkdUiMj/kWMLPj4j0C1x/qYj0S1FcI0RkceCx3xSRgsDx1iJSGfK8PRVymyMDv/+vA7FLCuJK+PeW7P9Xh7jGhsT0nYjMCRxP5/PllB/S9zdmjPHkB5ALfAMcAOwEzAXap+mxWwBHBC7vCnwFtAeGAv+KcP32gfgaAvsH4s5NYXzfAXuEHbsPGBy4PBi4N3C5B/AOIMCxwIw0/e5+AvZz4zkD/ggcAcyv6/MD7AZ8G/jcLHC5WQri6gY0CFy+NySu1qHXC7ufmYFYJRD7n1MQV0K/t1T8v0aKK+z79wO3ufB8OeWHtP2NeXnEfzTwtTHmW2PMVuAV4Ix0PLAxZqUxZnbg8gZgERBtV5QzgFeMMVuMMcuAr7Hxp9MZwHOBy88BfUKOP2+s6UCBiLRIcSynAN8YY6Kt1k7Zc2aMmQr8GuHxEnl+ugOTjTG/GmPWApOB05IdlzHmPWPMtsCX04F9ot1HILYmxpjpxmaP50N+lqTFFYXT7y3p/6/R4gqM2s8BXo52Hyl6vpzyQ9r+xryc+AuBH0K+/pHoyTclRKQ10AmYETh0TeDt2rPBt3KkP1YDvCcis0Skf+DYXsaYlYHLPwF7uRQbwHns+A+ZCc9Zos+PG8/b37Ajw6D9RaRMRD4WkRMDxwoDsaQjrkR+b+l+vk4EVhljloYcS/vzFZYf0vY35uXE7zoRaQy8AQwwxqwHngTaAB2Bldi3mm44wRhzBPBn4GoR+WPoNwMjG1fqfEVkJ6A38FrgUKY8Z9u5+fw4EZGbgW3AmMChlUArY0wn4J/ASyLSJI0hZdzvLcz57Di4SPvzFSE/bJfqvzEvJ/5yYN+Qr/cJHEsLEcnD/lLHGGPGARhjVhljqo0xNcB/+H1qIq2xGmPKA59XA28G4lgVnMIJfF7tRmzYF6PZxphVgRgz4jkj8ecnbfGJyF+B04ELAwmDwFTKL4HLs7Dz538IxBA6HZSSuOrwe0vn89UAOAsYGxJvWp+vSPmBNP6NeTnxfwEcKCL7B0aR5wHj0/HAgfnDUcAiY8wDIcdD58bPBILVBuOB80SkoYjsDxyIPaGUith2EZFdg5exJwfnB2IIVgX0A94Kie2SQGXBscC6kLejqbDDSCwTnrOQx0vk+ZkEdBORZoFpjm6BY0klIqcBNwC9jTGbQo43F5HcwOUDsM/Pt4HY1ovIsYG/00tCfpZkxpXo7y2d/69/AhYbY7ZP4aTz+XLKD6Tzb6w+Z6cz/QN7Nvwr7Kv3zWl83BOwb9O+BOYEPnoALwDzAsfHAy1CbnNzIM4l1LNqIEZsB2ArJuYCC4LPC7A78AGwFHgf2C1wXIDHA7HNA4pSGNsuwC9A05BjaX/OsC88K4Eq7LzpZXV5frBz7l8HPi5NUVxfY+d5g39nTwWu2zfw+50DzAZ6hdxPETYRfwM8RmAFf5LjSvj3luz/10hxBY6PBq4Mu246ny+n/JC2vzFt2aCUUj7j5akepZRSEWjiV0opn9HEr5RSPqOJXymlfEYTv1JK+YwmfuVbIvKgiAwI+XqSiDwT8vX9IvJPh9veISJ/inH/Q0XkXxGOF4jIVfUIXal60cSv/Owz4HgAEckB9gA6hHz/eGBapBsaY24zxrxfx8ctADTxK9do4ld+Ng04LnC5A3aRzobASsiGwMGACTTtmhV4RxBcUj9aRM4OXO4htif+LLF9098OeYz2IvKRiHwrIv8IHBsOtBHb931EWn5SpUI0cDsApdxijFkhIttEpBV2dP85trvhccA6bLvcB4EzjDFrRORc4G7sakkARGRn4Gngj8aYZSIS3ub3IOBkbN/1JSLyJLbX+iHGmI4p/QGVcqCJX/ndNGzSPx54AJv4j8cm/nJs/5PJtr0KudgWAKEOwvZ0WRb4+mWgf8j3JxpjtgBbRGQ1v7faVco1mviV3wXn+Q/FTvX8AAwE1gMfAYXGmOMcbx3blpDL1ej/nMoAOsev/G4atqXxr8a2Ef4Ve/L1OOzovbmIHAe2la6IdAi7/RLggMCGGgDnxvGYG7BTP0q5QhO/8rt52Gqe6WHH1hm7X8HZwL0iMhfbRfH40BsbYyqxFTrvisgsbFJfF+0Bje37/pmIzNeTu8oN2p1TqXoSkcbGmI2BPuuPA0uNMQ+6HZdSTnTEr1T9XS4ic7D93Jtiq3yUylg64ldKKZ/REb9SSvmMJn6llPIZTfxKKeUzmviVUspnNPErpZTP/D8EjE4FGUMFmQAAAABJRU5ErkJggg==\n"
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**5) By doing a 60/40 partition for our training and test sets, we are unfortunately neglecting a significant portion of data that could come in handy for training our model. This is where cross-validation becomes highly useful, as it allows us to train on a greater proportion of the total dataset.**\n",
|
|
"\n",
|
|
"**Use cross-validation to train the model. Check the performance of the model for both a 10 folds cross-validation, and a leave-one-out cross-validation. In both cases, evaluate the models on the MSE. Report the mean and the variance of the test MSE across each folds. What do you observe ? Check the classes 'cross_validate()' and 'LeaveOneOut() of the scikit-learn library.**"
|
|
],
|
|
"metadata": {
|
|
"id": "SQivk8PmdDND",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"cv_results = cross_validate(model, x, y, cv=10, scoring= 'neg_mean_squared_error')\n",
|
|
"\n",
|
|
"print(sorted(cv_results.keys()))\n",
|
|
"print(-cv_results['test_score'])\n",
|
|
"print(np.mean(-cv_results['test_score']))\n",
|
|
"print(np.var(-cv_results['test_score']))"
|
|
],
|
|
"metadata": {
|
|
"id": "mr3iWrn3eHfk",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 6,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"['fit_time', 'score_time', 'test_score']\n",
|
|
"[22.50851837 30.97175756 9.30394444 0.7971845 5.09303892 1.17718135\n",
|
|
" 0.39336786 3.97656728 33.84484856 23.54644202]\n",
|
|
"13.161285086064387\n",
|
|
"156.40891319174224\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"cv_results = cross_validate(model, x, y, cv=LeaveOneOut(), scoring='neg_mean_squared_error')\n",
|
|
"\n",
|
|
"print(np.mean(-cv_results['test_score']))\n",
|
|
"print(np.var(-cv_results['test_score']))"
|
|
],
|
|
"metadata": {
|
|
"id": "SgYminmhfaGa",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 7,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"10.020378035919238\n",
|
|
"295.30059627461566\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**6) Maybe a linear model is not the best option to model the relationship between the variables 'Weight' and 'Height'. Instead of fitting a linear model, let's fit a polynomial model of specified degree. First, transform the 'Weight' feature to quadratic, and fit a simple linear regression model to predict the variable 'Height' using a 10 folds cross-validation, and evaluate on the MSE. Report the mean test MSE across each folds. Check the method PolynomialFeatures() to transform the input variable.**\n"
|
|
],
|
|
"metadata": {
|
|
"id": "U_PEwHV0f14c",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"pf = PolynomialFeatures(2)\n",
|
|
"x_quad = pf.fit_transform(x)\n",
|
|
"cv_results = cross_validate(model, x_quad, y, cv=10, scoring= 'neg_mean_squared_error')\n",
|
|
"\n",
|
|
"print(np.mean(-cv_results['test_score']))\n",
|
|
"print(np.var(-cv_results['test_score']))"
|
|
],
|
|
"metadata": {
|
|
"id": "o1a-6zXngJgZ",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 8,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"7.8976789878676\n",
|
|
"33.74661480774945\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**7) Let's now see how the MSE vary if we increase the model's complexity. For polynomial degrees between 1 and 20, repeatedly fit the model using a simple linear regression model and cross-validation, and plot the evolution of the mean train and test MSE in function of the polynomial degree. What do you observe ? Use the function 'PolynomialRegression' as the 'model' argument of the class 'cross_validate', and make sure to understand what the 'make_pipeline' method does.**\n",
|
|
"\n",
|
|
"**Moreover, can you identify which degree would be the best choice to model the relationship between the variables 'Weigth' and 'Height' ?**"
|
|
],
|
|
"metadata": {
|
|
"id": "aZm20YCJjb2j",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"def PolynomialRegression(degree=2, **kwargs):\n",
|
|
" return make_pipeline(PolynomialFeatures(degree),\n",
|
|
" LinearRegression(**kwargs))"
|
|
],
|
|
"metadata": {
|
|
"id": "QwZ91HOXjpEd",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 9,
|
|
"outputs": []
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"degree = [i for i in range(1,20)]\n",
|
|
"for i in degree:\n",
|
|
" pf = PolynomialFeatures(i)\n",
|
|
" x_quad = pf.fit_transform(x)\n",
|
|
" cv_results = cross_validate(model, x_quad, y, cv=10, scoring= 'neg_mean_squared_error')\n",
|
|
" print('Degree : '+str(i))\n",
|
|
" print(np.mean(-cv_results['test_score']))\n"
|
|
],
|
|
"metadata": {
|
|
"id": "ZpaKtxAwk_Xf",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 10,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Degree : 1\n",
|
|
"13.16128508606439\n",
|
|
"Degree : 2\n",
|
|
"7.8976789878676\n",
|
|
"Degree : 3\n",
|
|
"8.449587815628941\n",
|
|
"Degree : 4\n",
|
|
"8.761140882830777\n",
|
|
"Degree : 5\n",
|
|
"9.209181297075522\n",
|
|
"Degree : 6\n",
|
|
"9.9944500787517\n",
|
|
"Degree : 7\n",
|
|
"31.667574181416388\n",
|
|
"Degree : 8\n",
|
|
"286.7772385900282\n",
|
|
"Degree : 9\n",
|
|
"305.3532104481505\n",
|
|
"Degree : 10\n",
|
|
"1027.2250799743188\n",
|
|
"Degree : 11\n",
|
|
"3601.914162149535\n",
|
|
"Degree : 12\n",
|
|
"12496.404587249188\n",
|
|
"Degree : 13\n",
|
|
"41988.84494267203\n",
|
|
"Degree : 14\n",
|
|
"136516.8400021473\n",
|
|
"Degree : 15\n",
|
|
"433044.73757164006\n",
|
|
"Degree : 16\n",
|
|
"1353010.4991111865\n",
|
|
"Degree : 17\n",
|
|
"4189866.8990477705\n",
|
|
"Degree : 18\n",
|
|
"13134590.2350142\n",
|
|
"Degree : 19\n",
|
|
"38576291.29159396\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**8) Maybe selecting more than one feature might also result in a better fit. Use a forward selection strategy to select to best subset of 3 predictors for predicting the target variable 'Height'. See the class 'SequentialFeaturesSelector()' of the scikit-learn library. Use the scorer below as a scoring function.**\n",
|
|
"\n",
|
|
"**Report the selected variables, and fit a multiple linear regression model to predict the variable 'Height' using a 10 folds cross-validation. Report the mean MSE across each folds.** "
|
|
],
|
|
"metadata": {
|
|
"id": "zdqIb5hmlgf3",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [
|
|
"X = df.loc[:, df.columns != 'Height']\n",
|
|
"\n",
|
|
"sfs = SequentialFeatureSelector(model, n_features_to_select=3, scoring= 'neg_mean_squared_error')\n",
|
|
"fit = sfs.fit(X,y)\n",
|
|
"\n",
|
|
"#print(sfs.transform(X))\n",
|
|
"\n",
|
|
"cv_results = cross_validate(model, sfs.transform(X), y, cv=10, scoring= 'neg_mean_squared_error')\n",
|
|
"\n",
|
|
"print(sorted(cv_results.keys()))\n",
|
|
"print(-cv_results['test_score'])\n",
|
|
"print(np.mean(-cv_results['test_score']))\n",
|
|
"print(np.var(-cv_results['test_score']))"
|
|
],
|
|
"metadata": {
|
|
"id": "-vRZNmkel8j3",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": 33,
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
" Width Species_Bream Species_Parkki\n",
|
|
"0 4.0200 1 0\n",
|
|
"1 4.3056 1 0\n",
|
|
"2 4.6961 1 0\n",
|
|
"3 4.4555 1 0\n",
|
|
"4 5.1340 1 0\n",
|
|
".. ... ... ...\n",
|
|
"154 1.3936 0 0\n",
|
|
"155 1.2690 0 0\n",
|
|
"156 1.2558 0 0\n",
|
|
"157 2.0672 0 0\n",
|
|
"158 1.8792 0 0\n",
|
|
"\n",
|
|
"[159 rows x 3 columns]\n",
|
|
"['fit_time', 'score_time', 'test_score']\n",
|
|
"[2.28289699 8.07556819 0.77142363 0.86577968 0.64646845 0.2546358\n",
|
|
" 0.41165865 0.65943568 1.03734828 0.29821018]\n",
|
|
"1.5303425527389796\n",
|
|
"5.060370573803941\n"
|
|
]
|
|
}
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"source": [
|
|
"**9) Usually, models have more than one hyper-parameter that be tuned in order to find the model that best captures the relationship between our input and target variables. For instance, in the case of a simple linear regression using a polynomial transformation on the input variables, the hyper-parameter space would be the polymial's degree, and whether or not to fit the intercept. Inspecting each combination of of hyper-parameters and selecting the combination that results in the best model is called Grid Search. However, manually inspecting each combination could be a very tedious task, but fortunatly, scikit-learn provides a class GridSearchCV() that implements this protocol for you.**\n",
|
|
"\n",
|
|
"**Select 'Weight' as predictor and 'Height' as target variable, and perform a Grid Search on the hyper-parameters space of a simple linear regression using a polynomial transformation of input variable. Search for degrees varying between 1 and 20, and whether or not the intercept should be fit. Report the best subset of hyper-parameters, as well as the test MSE for a model fitted using this subset. Plot the obtained regression curve**"
|
|
],
|
|
"metadata": {
|
|
"id": "s0qZ8rGFpBKp",
|
|
"pycharm": {
|
|
"name": "#%% md\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"outputs": [],
|
|
"source": [],
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
}
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"source": [],
|
|
"metadata": {
|
|
"id": "BJJ0uaGlrtAk",
|
|
"pycharm": {
|
|
"name": "#%%\n"
|
|
}
|
|
},
|
|
"execution_count": null,
|
|
"outputs": []
|
|
}
|
|
]
|
|
} |