diff --git a/doc/src/week37/programs/codeexamplesscaling.do.txt b/doc/src/week37/programs/codeexamplesscaling.do.txt index 85b35c90e..1060014d4 100644 --- a/doc/src/week37/programs/codeexamplesscaling.do.txt +++ b/doc/src/week37/programs/codeexamplesscaling.do.txt @@ -310,21 +310,23 @@ for p in range(d): #Split data in train and test X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) -print(X_train) # Scale data by subtracting mean value using scikit-learn -from sklearn.preprocessing import StandardScaler scaler = StandardScaler() scaler.fit(X_train) X_train_scaled = scaler.transform(X_train) -print(X_train_scaled) +X_test_scaled = scaler.transform(X_test) + +y_train_scaled = y_train - np.mean(y_train) +y_test_scaled = y_test - np.mean(y_test) + #Calculate beta OLS = LinearRegression() -OLS.fit(X_train,y_train_scaled) -ypredictOLS = OLS.predict(X_test) +OLS.fit(X_train_scaled,y_train_scaled) +ypredictOLS = OLS.predict(X_test_scaled) RegRidge = linear_model.Ridge(Lambda) -RegRidge.fit(X_train_scaled,y_train) -ypredictRidge = RegRidge.predict(X_test) +RegRidge.fit(X_train_scaled,y_train_scaled) +ypredictRidge = RegRidge.predict(X_test_scaled) print(OLS.coef_) print(RegRidge.coef_) print(OLS.intercept_) diff --git a/doc/src/week37/programs/codeexamplesscaling.ipynb b/doc/src/week37/programs/codeexamplesscaling.ipynb index 4fb54d1f8..ace108d34 100644 --- a/doc/src/week37/programs/codeexamplesscaling.ipynb +++ b/doc/src/week37/programs/codeexamplesscaling.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "a6d58753", + "id": "8dd69296", "metadata": {}, "source": [ " 22\u001b[0m X_train_scaled \u001b[38;5;241m=\u001b[39m \u001b[43mscaler\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtransform\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX_train\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;66;03m#Calculate beta\u001b[39;00m\n\u001b[1;32m 26\u001b[0m OLS \u001b[38;5;241m=\u001b[39m LinearRegression()\n", - "File \u001b[0;32m~/miniforge3/envs/myenv/lib/python3.9/site-packages/sklearn/preprocessing/_data.py:970\u001b[0m, in \u001b[0;36mStandardScaler.transform\u001b[0;34m(self, X, copy)\u001b[0m\n\u001b[1;32m 955\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mtransform\u001b[39m(\u001b[38;5;28mself\u001b[39m, X, copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m 956\u001b[0m \u001b[38;5;124;03m\"\"\"Perform standardization by centering and scaling.\u001b[39;00m\n\u001b[1;32m 957\u001b[0m \n\u001b[1;32m 958\u001b[0m \u001b[38;5;124;03m Parameters\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 968\u001b[0m \u001b[38;5;124;03m Transformed array.\u001b[39;00m\n\u001b[1;32m 969\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 970\u001b[0m \u001b[43mcheck_is_fitted\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 972\u001b[0m copy \u001b[38;5;241m=\u001b[39m copy \u001b[38;5;28;01mif\u001b[39;00m copy \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcopy\n\u001b[1;32m 973\u001b[0m X \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_validate_data(\n\u001b[1;32m 974\u001b[0m X,\n\u001b[1;32m 975\u001b[0m reset\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 980\u001b[0m force_all_finite\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mallow-nan\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m 981\u001b[0m )\n", - "File \u001b[0;32m~/miniforge3/envs/myenv/lib/python3.9/site-packages/sklearn/utils/validation.py:1222\u001b[0m, in \u001b[0;36mcheck_is_fitted\u001b[0;34m(estimator, attributes, msg, all_or_any)\u001b[0m\n\u001b[1;32m 1217\u001b[0m fitted \u001b[38;5;241m=\u001b[39m [\n\u001b[1;32m 1218\u001b[0m v \u001b[38;5;28;01mfor\u001b[39;00m v \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mvars\u001b[39m(estimator) \u001b[38;5;28;01mif\u001b[39;00m v\u001b[38;5;241m.\u001b[39mendswith(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m_\u001b[39m\u001b[38;5;124m\"\u001b[39m) \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m v\u001b[38;5;241m.\u001b[39mstartswith(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m__\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 1219\u001b[0m ]\n\u001b[1;32m 1221\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m fitted:\n\u001b[0;32m-> 1222\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m NotFittedError(msg \u001b[38;5;241m%\u001b[39m {\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mname\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;28mtype\u001b[39m(estimator)\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m})\n", - "\u001b[0;31mNotFittedError\u001b[0m: This StandardScaler instance is not fitted yet. Call 'fit' with appropriate arguments before using this estimator." + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.5 5. ]\n", + "[0.50026997 4.99920486]\n", + "1.9999999999999893\n", + "2.00251547683159\n", + " \n", + "test MSE of OLS\n", + "2.524354896707238e-28\n", + " \n", + "test MSE of Ridge\n", + "1.5906811889393548e-05\n" ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABY30lEQVR4nO3deVxU5eIG8OfMAMMOosCAgOIu4oZobuUWuJRZllq23rrllmVWdsv7S7yVpZWZWbbcm3kzMq3MvBVJLliZsYkbigsgyiIqyj7DLOf3x8Ak4jKDM3Nmeb6fj58bw+Hw8DqXeXzPO+8RRFEUQURERGQjMqkDEBERkWth+SAiIiKbYvkgIiIim2L5ICIiIpti+SAiIiKbYvkgIiIim2L5ICIiIpti+SAiIiKbcpM6wOX0ej1KSkrg5+cHQRCkjkNEREQmEEUR1dXVCA8Ph0x27bkNuysfJSUliIyMlDoGERERtcKpU6cQERFxzWPMKh9JSUlYvHhxs8dCQ0NRVlYGwNB6Fi9ejI8//hgXLlzATTfdhPfffx+9evUy+Xv4+fkZw/v7+5sT77o0Gg22bt2KxMREuLu7W/TczoZjZTqOlek4VubheJmOY2U6a41VVVUVIiMjja/j12L2zEevXr3wyy+/GD+Wy+XG/162bBmWL1+Ozz77DN26dcOrr76KhIQE5OXlmRQGgPFSi7+/v1XKh7e3N/z9/fnkvA6Olek4VqbjWJmH42U6jpXprD1WpiyZMHvBqZubG5RKpfFPcHAwAMOsx4oVK7Bw4UJMnjwZsbGxWLt2Lerq6pCcnGx+eiIiInJKZpePY8eOITw8HNHR0bj33nuRn58PACgoKEBZWRkSExONxyoUCowYMQK7d++2XGIiIiJyaGZddrnpppvw3//+F926dcOZM2fw6quvYujQoTh06JBx3UdoaGizrwkNDcXJkyevek61Wg21Wm38uKqqCoBhWkij0ZgT77qazmfp8zojjpXpOFam41iZh+NlOo6V6aw1VuacTxBFUWztN6qtrUXnzp2xYMECDB48GMOGDUNJSQnCwsKMxzz++OM4deoUUlJSrniOKy1iBYDk5GR4e3u3NhoRERHZUF1dHaZPn47Kysrrrtm8obfa+vj4oHfv3jh27BjuvPNOAEBZWVmz8lFeXt5iNuRSL774IubPn2/8uGm1bGJiolUWnKampiIhIYELkq6DY2U6jpXpOFbm4XiZjmNlOmuNVdOVC1PcUPlQq9U4fPgwbr75ZkRHR0OpVCI1NRX9+/cHADQ0NCAtLQ1Lly696jkUCgUUCkWLx93d3a32BLLmuZ0Nx8p0HCvTcazMw/EyHcfKdJYeK3POZVb5eO655zBx4kRERUWhvLwcr776KqqqqvDwww9DEATMmzcPS5YsQdeuXdG1a1csWbIE3t7emD59utk/BBERETkns8rH6dOncd999+HcuXMIDg7G4MGDsWfPHnTo0AEAsGDBAtTX12P27NnGTca2bt1q8h4fRERE5PzMKh/r16+/5ucFQUBSUhKSkpJuJBMRERFZgU4vIr2gAgCQXlCBwV1CIJfZ/j5qvKstERGRC0g5WIrhS7fj0bUZAIBH12Zg+NLtSDlYavMsLB9EREROLuVgKWaty0ZpparZ42WVKsxal23zAsLyQURE5MR0ehGLt+TiSpt6NT22eEsudPpWb/tlNpYPIiIiJ5ZeUNFsxkOjB3afEdC0xagIoLRSZVwLYgssH0RERE6svLr5pRa9KODrgpYv/5cfZ00sH0RERE4sxM+zxWO9g0QIwvWPsxaWDyIiIic2KDoIYQHNi8Xg4L/WdwgAwgI8MSg6yGaZWD6IiIicmFwmYNHEmEseEdE90FA+miY/Fk2Msel+HywfRERETm5cbBhiww03a5ULQFPPUAZ4YvUDcRgXG3aNr7a8G7qxHBEREdm/M1Uq5JYa7jq7dHJvoCQHnz48kDucEhERkXV8m10MvQjEd2iDSf3CARjWgkhRPACWDyIiIqcmiiI2Zp4CAEyNj5Q4jQHLBxERkRPLLrqA/HO18HKXY0If267tuBqWDyIiIie2MfM0AGBC7zD4KuxjqSfLBxERkZOqa9Dif/sNN42bEh8hcZq/sHwQERE5qZSDZahRaxEV5I2bbLiJ2PWwfBARETmpDY0LTacMiIBw+X7qEmL5ICIickJF5+uwJ78CggDcPcB+LrkALB9ERERO6etsw0LT4V3aITzQS+I0zbF8EBERORm9XsQ3WYbycY+dzXoALB9ERERO54/88yi+WA8/TzeM7aWUOk4LLlM+dHoR6QUVAID0ggro9OJ1voKIiMgxNS00ndQvHJ7uconTtOQS5SPlYCmGL92OR9dmAAAeXZuB4Uu3I+VgqcTJiIiILKuyXoOUg2UAgCkD7GM79cs5fflIOViKWeuyUVqpavZ4WaUKs9Zls4AQEZFT+d/+Eqi1enQL9UWfiACp41yRU5cPnV7E4i25uNIFlqbHFm/J5SUYIiJyGk3bqU8ZEGlXe3tcyqnLR3pBRbMZD60e+OjwXz+yCKC0UmVcC0JEROTIjp2pRs6pi5DLBNzZv73Uca7KqctHeXXzSy06Eci9KMPlEx2XH0dEROSINja+vXZ0jxAE+ykkTnN1Tl0+Qvw8m30sa5x90onXPo6IiMjRaHR6fJtdDMCwnbo9c+ryMSg6CGEBnmi64iVv/I+mmQ8BQFiAJwbZ0c12iIiIWiMt7yzO1ajRztcDo3qESB3nmpy6fMhlAhZNjAFgKBoCgLYKsfG/DBZNjIFcZp8LcoiIiEy1Mcuwt8ed/drDXW7fL+/2nc4CxsWGYfUDcVAGeEIQgMEhegCAu1zA6gfiMC42TOKEREREN+Z8jRrbDpcDAKbE2+feHpdykzqALYyLDUNCjBJ7jpfjxL49+Ok0oNGJ6BLiK3U0IiKiG7ZpbzG0ehF9IwLQXekndZzrcvqZjyZymYBB0UEI8ABGdQ8GAKxPPyVxKiIiohsjiiK+brqJnAPMegAuVD4uNTXesAr4m+zTUGt1EqchIiJqvYPFVThSVg0PNxnu6BMudRyTuGT5uKVLWyj9PXGhToOth85IHYeIiKjVmhaaju2lRIC3u8RpTOOS5cNNLjPOfqzPKJI4DRERUeuoNDpszikBYP97e1zKJcsHAEwdGAlBAH4/fh4nz9dKHYeIiMhsqblnUFmvQViAJ4Z1aSd1HJO5bPmIaOONm7saFp5+lcGFp0RE5HiatlO/Z0CEQ+1Z5bLlAwDuG2hYFbwx6zQ0Or3EaYiIiExXWlmPX4+dBWAoH47EpcvHmJ6haOfrgbPVamw/Ui51HCIiIpN9m10MUTTcSqRDWx+p45jFpcuHh5sMdze2RV56ISIiRyGKIjZmGl63HGmhaROXLh8AMK1xQ5adeeUorayXOA0REdH1ZRReQOH5Onh7yDGht+PdJsTly0enYF/cFB0EvQhsyDgtdRwiIqLrapr1uL1PGHwUjnenFJcvHwBw36AoAMCGzFPQ6UWJ0xAREV1drVqLHw6UAnCMm8hdCcsHgHGxSgR4uaP44l8rh4mIiOzRjwdKUdegQ3Q7H8R3aCN1nFZh+QDg6S7HXf3bA+DN5oiIyL5dureHIDjO3h6XYvlo1HTp5ZfDZ3C2Wi1xGiIiopYKz9UivaACMgGYHNde6jitxvLRqLvSD/2jAqHV/3VrYiIiInvS9Pp0c9dghAV4SZym9Vg+LnHfQMPsx1cZRRBFLjwlIiL7odOL+CbbUD6mxDve3h6XYvm4xO19w+CrcEPh+Tr8kX9e6jhERERGvx8/h9JKFQK83HFrz1Cp49wQlo9LeHu44Y5+4QC48JSIiOxL00LTSf3C4ekulzjNjWH5uEzTpZeUg2W4UNsgcRoiIiKgsk6Dnw+VAQCmDHDMvT0uxfJxmdj2/ogJ80eDTo9Ne4uljkNERITv9xWjQatHD6UfYtv7Sx3nhrF8XEYQBNw3yNAq13PhKRER2YGmSy5T4iMddm+PS7F8XMGk/u3h6S7D0TM1yC66KHUcIiJyYXll1dh/uhJuMgF3Nq5LdHQsH1fg7+mO23o3LTwtkjgNERG5sqabyI3pGYK2vgqJ01gGy8dVNF16+d/+UlSrNBKnISIiV6S5ZP2hMyw0bcLycRUDOrRBlxBf1Gt02JxTInUcIiJyQduPlON8bQOC/RQY2T1Y6jgWw/JxFYIg4N6Bfy08JSIisrWNmYaFppP7t4eb3Hlesp3nJ7GCyXER8JDLcLC4CgeLK6WOQ0RELuRstRo78soBOP526pdj+biGIB8PjI1VAgC+5MJTIiKyoe/2FkOnF9EvMhBdQvykjmNRLB/XcV/jpZfNOSWoa9BKnIaIiFyBKIrY0PguF2eb9QBusHy8/vrrEAQB8+bNMz4miiKSkpIQHh4OLy8vjBw5EocOHbrRnJIZ3KktOrT1Ro1ai//tL5U6DhERuYB9pytxrLwGCjcZJvZ1jr09LtXq8pGRkYGPP/4Yffr0afb4smXLsHz5cqxatQoZGRlQKpVISEhAdXX1DYeVgkwmYFrTwlNeeiEiIhto2ttjfKwS/p7uEqexvFaVj5qaGtx///345JNP0KZNG+PjoihixYoVWLhwISZPnozY2FisXbsWdXV1SE5OtlhoW7snLgJymYDsoos4esYxSxQRETkGlUaH7/cZtniYEu88e3tcyq01XzRnzhzcdtttuPXWW/Hqq68aHy8oKEBZWRkSExONjykUCowYMQK7d+/GjBkzWpxLrVZDrVYbP66qqgIAaDQaaDSW3dyr6XzmnreNlxyjuwcj9XA5kvcUYuGEHhbNZY9aO1auiGNlOo6VeThepnOmsfpxfymqVVq0D/REfKS/3bwWmnpeU5hdPtavX4/s7GxkZGS0+FxZmeF2v6Ghoc0eDw0NxcmTJ694vtdffx2LFy9u8fjWrVvh7e1tbjyTpKammv010aIAQI4NGScRq8+Hu4ss1W3NWLkqjpXpOFbm4XiZzhnG6sNcGQAZYn3rkJLyk9W+j6XHqq6uzuRjzSofp06dwtNPP42tW7fC09Pzqsddfsc9URSvehe+F198EfPnzzd+XFVVhcjISCQmJsLf37K3DdZoNEhNTUVCQgLc3c27hjZWL2LL8l9RWqmCLKo/JvQJs2g2e3MjY+VqOFam41iZh+NlOmcZq5KL9Ti251cAwAtTb0FkG8v/I9xaY9V05cIUZpWPrKwslJeXY8CAAcbHdDoddu3ahVWrViEvLw+AYQYkLOyvF+fy8vIWsyFNFAoFFIqWN8pxd3e32hOoNed2BzA1PhLvbjuGjVklmDwgyirZ7I01/x6cDcfKdBwr83C8TOfoY7V5fyFEERjSqS06hQRY9XtZeqzMOZdZFw/GjBmDAwcOICcnx/gnPj4e999/P3JyctCpUycolcpmUzkNDQ1IS0vD0KFDzflWdmnqwEgIAvBH/nkUnquVOg4RETkRvV7E11mG7dSdcW+PS5k18+Hn54fY2Nhmj/n4+KBt27bGx+fNm4clS5aga9eu6Nq1K5YsWQJvb29Mnz7dcqkl0j7QCyO6BWNn3lmszziFf4x3/oWnRERkG+mFFSiqqIOvwg3jY5370r7Fl00uWLAA8+bNw+zZsxEfH4/i4mJs3boVfn7OsTXsvQMNl1u+zjoNjU4vcRoiInIWTTua3t4nDF4econTWFer3mp7qZ07dzb7WBAEJCUlISkp6UZPbZfG9AxBO18FztWose3wGYxz8nZKRETWV6PW4qcDhneMOuveHpdykTeMWo67XGa8Fvdl+imJ0xARkTP4YX8J6jU6dAr2QVxUoNRxrI7loxWmNbbSXcfO4vQF09/XTEREdCUbMxsXmg6IvOrWFM6E5aMVOrbzwZBObSGKwIbGJwwREVFr5J+tQebJC5AJwOS49lLHsQmWj1a6d5Bh9mNj5ino9KLEaYiIyFFtbHx77YhuwQj1v/oGns6E5aOVxvZSItDbHaWVKuw6elbqOERE5IB0ehHfZhvKx1QXWGjahOWjlTzd5Zjcv2nhaZHEaYiIyBHtOnYWZ6rUaOPtjjE9r7wTuDNi+bgB9zVeetl2pBzlVSqJ0xARkaP5unHd4KR+7eHh5jovya7zk1pB11A/DOjQBjq9aLxmR0REZIoLtQ1IzT0DwPm3U78cy8cNunegYfbjq4xT0HPhKRERmWhzTjEadHrEhPmjV7h1byJnb1g+btBtfcLgp3BDUUUd/sg/L3UcIiJyEE0z5lNdbNYDYPm4Yd4ebpjUPxwAF54SEZFpckuqcKikCh5yGSb1c429PS7F8mEBTTeb23roDCpqGyROQ0RE9m5jluH2HLfGhKCNj4fEaWyP5cMCYtsHoHf7ADTo9Mb3axMREV1Jg1aP7/YWAzBsp+6KWD4sZFrjwtMv04sgilx4SkREV7bt8BlcqNMgxE+Bm7u2kzqOJFg+LGRSv3B4uctx4mwtMk9ekDoOERHZqaaFpncPiICb3DVfhl3zp7YCP0933N4nDAAXnhIR0ZWVV6mwM68cADBlgOu9y6UJy4cF3TvIsPD0xwOlqKzXSJyGiIjszbd7i6EXgQEd2qBTsK/UcSTD8mFBcVGB6BbqC5VGj+9ziqWOQ0REdkQURWzMNLzLxZVnPQCWD4sSBMH4ttsv009x4SkRERllF13EibO18HSX4bbGy/SuiuXDwibHGW4OlFtahQPFlVLHISIiO/F1494eE3qHwc/TXeI00mL5sLBAbw+Mj1UCMMx+EBER1TfosGVfKQDX3dvjUiwfVtB06eX7nGLUqrUSpyEiIqmlHCpFjVqLyCAv3BQdJHUcybF8WMHgTkGIbueD2gYd/re/ROo4REQksY2Zhr097omLhEwmSJxGeiwfViAIwiU7nvLSCxGRKztVUYfdJ85DEIC7B7jeTeSuhOXDSu6Oi4CbTEDOqYs4UlYldRwiIpLI1407mg7r3A4RbbwlTmMfWD6sJNhPgVt7hgIA1nP2g4jIJen1orF8TIl37b09LsXyYUX3DjJcevk2+zRUGp3EaYiIyNb25J9H8cV6+Hm6YWwvpdRx7AbLhxXd3DUY7QO9UKXS4qeDpVLHISIiG2u6idzEvuHwdJdLnMZ+sHxYkVwmYGq8YfaDl16IiFxLlUqDHw807e3BSy6XYvmwsqkDIyATgD8LKpB/tkbqOEREZCP/21cKtVaPriG+6BcZKHUcu8LyYWVhAV4Y2T0EAPBVBmc/iIhcxcbG7dSnxEdAELi3x6VYPmzg3sY9P77OOo0GrV7iNEREZG3Hy6uxt+gi5DIBd/bn3h6XY/mwgdE9QhDip8D52gb8cviM1HGIiMjKmhaajuoejBA/T4nT2B+WDxtwk8uM7+/+Mr1I4jRERGRNWp0e32YXAwDu4U3krojlw0amxRtuNvfb8XM4VVEncRoiIrKWtKNncbZajSAfD4zuESJ1HLvE8mEjUW29MbxLO4gisCGTC0+JiJxV003k7urfHh5ufJm9Eo6KDTXdbG5D5ilodVx4SkTkbCpqG7DtiGFtH7dTvzqWDxtK7BWKNt7uOFOlxs68s1LHISIiC/tubzE0OhG92wegh9Jf6jh2i+XDhhRuctwdZ2jC6zO48JSIyJmIomi8rM5Zj2tj+bCxppvNbT9SjrJKlcRpiIjIUg6VVOFIWTU85DLc0Tdc6jh2jeXDxrqE+GFgxzbQi8BGLjwlInIaTb/TE3uFItDbQ+I09o3lQwL3DjS87farzFPQ60WJ0xAR0Y1Sa3XYvK8EADAlnnt7XA/LhwQm9A6Dn6cbTl+ox+8nzkkdh4iIbtAvueW4WKdBWIAnhndpJ3Ucu8fyIQEvDznuatzrf306L70QETm6poWmk+PaQy7jTeSuh+VDIk2XXrbmluF8jVriNERE1FpllSr8esywfQK3UzcNy4dEYsL90TciABqdiG+yT0sdh4iIWumb7NPQi8CgjkGIbucjdRyHwPIhoXsHGWY/1mecgihy4SkRkaMRRRFfN97B9h7u7WEylg8JTewbDm8POfLP1iK9oELqOEREZKaskxdQcK4W3h5y3NY7TOo4DoPlQ0K+CjdM7GPYiGZ9BheeEhE5mqaFphN6h8FH4SZxGsfB8iGxph1PfzxQiso6jcRpiIjIVHUNWvywvxQAMGUAL7mYg+VDYv0iA9FD6Qe1Vo9Ne7nwlIjIUfx4oAy1DTp0bOuNQdFBUsdxKCwfEhMEAfcONMx+cOEpEZHjaNpO/Z4BERAE7u1hDpYPO3BX/wgo3GQ4UlaNnFMXpY5DRETXcfJ8Lf4sqIAgAJPjeMnFXCwfdiDA2x0TGldJf8WFp0REdq/p7bXDu7RDeKCXxGkcD8uHnWi69PL9vhLUqLUSpyEioqvR6UV801g+eBO51mH5sBODooPQKdgHdQ06bGm8MyIREdmf3SfOoaRSBX9PNyTGhEodxyGxfNiJZgtP04skTkNERFezMdMw6zGpX3t4usslTuOYWD7syN1xEXCXC9h3uhK5JVVSxyEiostU1mmQcqgMADCF26m3GsuHHWnrq0BijBIAsD6Dsx9ERPZCpxfxx4nzePWHXDRo9egW6ove7QOkjuWwWD7szLTGSy+b9hajvkEncRoiIko5WIrhS7fjvk/2YGPjQtPSShV+bpwBIfOxfNiZ4V3aIaKNF6pVWvx4oFTqOERELi3lYClmrctGaaWq2ePVKi1mrctGykH+nm4Nlg87I5MJmBbftOMpL70QEUlFpxexeEsurrXv9OItudDpuTO1ucwqH6tXr0afPn3g7+8Pf39/DBkyBD/99JPx86IoIikpCeHh4fDy8sLIkSNx6NAhi4d2dlPiIyETgIzCCzheXi11HCIil5ReUNFixuNSIgyXX9ILKmwXykmYVT4iIiLwxhtvIDMzE5mZmRg9ejQmTZpkLBjLli3D8uXLsWrVKmRkZECpVCIhIQHV1XwBNYcywBOje4QAANanc8dTIiIplFdfvXi05jj6i1nlY+LEiZgwYQK6deuGbt264bXXXoOvry/27NkDURSxYsUKLFy4EJMnT0ZsbCzWrl2Luro6JCcnWyu/07p3YBQA4Jvs01BrufCUiMjWQvw8LXoc/cWttV+o0+mwceNG1NbWYsiQISgoKEBZWRkSExONxygUCowYMQK7d+/GjBkzrngetVoNtVpt/LiqyrC/hUajgUajaW28K2o6n6XPaw3DOgUi1E+BM9VqpOwvwYTeSpt+f0caK6lxrEzHsTIPx8t01hir/hF+6NBGgZKLKmhEAYAIDxnQdANbAUCovyf6R/g51N+RtZ5X5pxPEM28h/uBAwcwZMgQqFQq+Pr6Ijk5GRMmTMDu3bsxbNgwFBcXIzw83Hj8E088gZMnT+Lnn3++4vmSkpKwePHiFo8nJyfD29vbnGhO54ciGbYWy9AtQI85MXqp4xARuRydCCzZK8c5tYA7onQY056LS6+mrq4O06dPR2VlJfz9/a95rNkzH927d0dOTg4uXryIb775Bg8//DDS0tKMnxeaKmEjURRbPHapF198EfPnzzd+XFVVhcjISCQmJl43vLk0Gg1SU1ORkJAAd3d3i57bGvpcqEfqO7/iaKUMsYNvQVSQ7cqYo42VlDhWpuNYmYfjZTprjdXX2cU4t+cQZAKQclqGn4sNjyv9PfGP8T1wa0/Hu7eLtcaq6cqFKcwuHx4eHujSpQsAID4+HhkZGXj33XfxwgsvAADKysoQFhZmPL68vByhoVf/y1EoFFAoFC0ed3d3t9r/2ax5bkuKDnHH8C7t8Ouxc/g2pxTPj+1h8wyOMlb2gGNlOo6VeTheprPkWGl0eqxOKwAAPD+2B/pFBqK8WoUQP08Mig6CXHb1f1g7Aks/r8w51w3v8yGKItRqNaKjo6FUKpGammr8XENDA9LS0jB06NAb/TYu675BhoWnGzNPQ6vjpRciIlvZlF2Mooo6tPP1wMNDO2BI57aY1K89hnRu6/DFQ2pmzXy89NJLGD9+PCIjI1FdXY3169dj586dSElJgSAImDdvHpYsWYKuXbuia9euWLJkCby9vTF9+nRr5Xd6t/YMRVsfD5RXq7H9SDkSe9l24SkRkSvS6PR4b8cxAMCMWzrD26PV78+gKzBrNM+cOYMHH3wQpaWlCAgIQJ8+fZCSkoKEhAQAwIIFC1BfX4/Zs2fjwoULuOmmm7B161b4+flZJbwr8HCT4e4BEfh4Vz7WZ5xi+SAisoFvsk7jVEU92vkq8MDgDlLHcTpmlY///Oc/1/y8IAhISkpCUlLSjWSiy0wbGImPd+VjZ145SivrERbgJXUkIiKn1aDV473txwEAM0d0gpeHXOJEzof3dnEAnYN9MSg6CHoR2JBxWuo4RERO7Zvs0yi+aJj1uP8mznpYA8uHg7hvkOFmcxsyT/EmRkREVtKg1WNV46zHrJGdOethJSwfDmJ8bBj8Pd1QfLEeH+/Kx+acYvxx4jyLCBGRBW3MOoXii/UI8VPg/puipI7jtLh810F4ussxoEMb7Mg7i6UpR4yPhwV4YtHEGIyLDbvGVxMR0fWotTq8f8msh6c7Zz2shTMfDiLlYCl25J1t8XhZpQqz1mUj5WCpBKmIiJzHhszTKKlUIcRPYdxjiayD5cMB6PQiFm/JveLnmi66LN6Sy0swREStpNbq8MEOw6zHbM56WB3LhwNIL6hAaaXqqp8XAZRWqpBeUGG7UERETmRDximUVqqg9PfEvZz1sDqWDwdQXn314tGa44iI6C8qjQ7v7zgBAJg9irMetsDy4QBC/DwtehwREf3lq4xTKKtSISzAE9MGRkodxyWwfDiAQdFBCAvwxNVuYyTA8K6XQdFBtoxFROTwVBodPtjZuNZjVBco3DjrYQssHw5ALhOwaGIMAFyxgIgAFk2M4V0WiYjMtD69CGeq1AgP8MTU+Aip47gMlg8HMS42DKsfiIMyoOWllanxEdzng4jITIZZj6a1Hpz1sCVuMuZAxsWGISFGifSCCpRXq7D/dCX+81sBth0uR41aC18F/zqJiEyV/GcRyqvVaB/ohanxXOthS5z5cDBymYAhndtiUr/2+Mf4HujY1hvnaxvw6W8FUkcjInIYKo0Oq9MMsx5zRnWBhxtfDm2Jo+3A3OUyzE/sDgD4ZFc+KmobJE5EROQY1u05ibONsx73DOBaD1tj+XBwt/cOQ0yYP6rVWqxuXLFNRERXV9+gw4dp+QCAuaM56yEFjriDk8kEPD/OMPux9o+TKK2slzgREZF9++LPkzhXo0ZEGy/czVkPSbB8OIGR3YIxKDoIDVo9Vm47JnUcIiK7VdegxYeNaz3mju4CdzlfBqXAUXcCgiDghcbZjw2Zp5F/tkbiRERE9mndnpM4V9OAqCBvTI7jrIdUWD6cxIAOQRjTIwQ6vYi3U49KHYeIyO7UNWjxUeNajyc56yEpjrwTeW5sdwgC8MP+UhwsrpQ6DhGRXfnvHydxvrYBHdp6Y3L/9lLHcWksH06kZ5g/JvUNBwAs+zlP4jRERPajVq3Fx7saZz1GdYEbZz0kxdF3MvMTusNNJmDX0bP448R5qeMQEdmF//5xEhW1DejY1ht3cdZDciwfTiaqrTfuGxQFAFj28xGIoihxIiIiadWotfh4V9M7XLpy1sMO8G/ACc0d3QVe7nLsLbqIXw6XSx2HiEhSa3cX4kKdBtHtfDCpX7jUcQgsH04pxN8TfxvWEQDw5s9HoNNz9oOIXFO1SoNPfjWs9XhqDNd62Av+LTipGbd0hr+nG46eqcHmnGKp4xARSeK/f5zExToNOrXzwR19udbDXrB8OKkAb3fMGtkFALA89SgatHqJExER2Va1SmN8h8tTY7pCLhMkTkRNWD6c2CNDOyLET4HTF+rxZXqR1HGIiGzqs98LUVmvQedgH0zsy7Ue9oTlw4l5ecjx1JiuAID3th9HrVorcSIiItuoarbWg7Me9oblw8lNGxiJDm29ca5GjTW/F0gdh4jIJtb8VogqlRZdQnxxex/Oetgblg8n5y6XYX5CNwDAR2n5uFDbIHEiIiLrqqzX4D+/GWY9nuash11i+XABE/uEo4fSD9Xqv24lTUTkrNb8XoAqlRZdQ3wxoXeY1HHoClg+XIBMJmDBuO4AgM92F6KsUiVxIiIi66iq1+A/vxkuMT99K2c97BXLh4sY1T0EAzu2gVqrx7vbjkkdh4jIKtbsPolqlRbdQ/0wIZazHvaK5cNFCIKABeN6AAA2ZJ5CwblaiRMREVlWnRb47A/DtgJP39oVMs562C2WDxcysGMQRvcIgU4vYnnqUanjEBFZ1I4SGWrUWvRQ+mFcL6XUcegaWD5czHOJhrUfW/aV4GBxpcRpiIgs42KdBmllhpmOp8dw1sPesXy4mJhwf9zRuNPfW1vzJE5DRGQZn+4uhFonoEeoL8Zy1sPusXy4oPkJ3eAmE7Az7yz+zD8vdRwiohtyobYB/21c6zF3dGfOejgAlg8X1LGdD6YNjAQALPs5D6IoSpyIiKj1Pvk1H7UNOrT3FpHQM0TqOGQClg8X9dSYrvB0lyHr5AVsP1IudRwiolapqG3A2t2FAIDxkXoIAmc9HAHLh4sK9ffEI0OjAQBv/pwHvZ6zH0TkeJpmPWLC/BDbhr/HHAXLhwubNaIz/D3dcKSsGt/vK5E6DhGRWc7XqI2zHk+N6gxOejgOlg8XFuDtjhkjOgMA3k7NQ4NWL3EiIiLTffxrPuoadOjdPgCjewRLHYfMwPLh4v42rCPa+SpwqqIeX2UUSR2HiMgk52rU+O/ukwCAebd25VoPB8Py4eK8Pdzw9JguAIB3tx1HXYNW4kRERNf38a581Gt06BsRgNE9+A4XR8PyQZg2MAqRQV44V6PGmt8LpY5DRHRN52rU+O8fhQCAebd246yHA2L5IHi4yfBsgmHb9Y/STqCyTiNxIiKiq/so7QRUGj36RgZiZHeu9XBELB8EALijbzh6KP1QpdJiddoJqeMQEV1RebUKn+/hWg9Hx/JBAACZTDDedO6z3QU4U6WSOBERUUsfpeVDpdGjX2QgRnbjrIejYvkgozE9QzCgQxuoNHq8vzNf6jhERM2UV6mwrnHW45kErvVwZCwfZCQIAhaMNcx+bMwqxjlOfhCRHfkwLR9qrR5xUYG4pWs7qePQDWD5oGZu6tQWI7sHQ6sX8eMpPj2IyD6UV6nwxZ9Naz046+Ho+OpCLTSt/cg6J8Ph0mqJ0xARAR/sPAG1Vo8BHdrgZs56ODyWD2ohtn0AbuutBAAs/+WYxGmIyNWVVaqQnG7YgfkZzno4BZYPuqJ5YzpDBhE7j55DRmGF1HGIyIWt3nkcDVo9BnZsg2Fd2kodhyyA5YOuqGNbHwwOMdyeeulPRyCKvFU1EdleaWU9vkw/BYCzHs6E5YOuamyEHgo3GTJPXsCOvHKp4xCRC1q98wQadHoMig7CkM6c9XAWLB90VYEK4MHBUQCAN38+Cr2esx9EZDslF+uxvnHWg7uZOheWD7qmGTdHw8/TDYdLq7Blf4nUcYjIhXyw8zgadHrcFB2EoZ35DhdnwvJB1xTo7Y4Zt3QCACxPPQqNTi9xIiJyBcUX6/FVRuNaj4RuEqchSzOrfLz++usYOHAg/Pz8EBISgjvvvBN5eXnNjhFFEUlJSQgPD4eXlxdGjhyJQ4cOWTQ02dbfhkWjna8CJ8/XGX8ZEBFZ0/s7jkOjEzGkU1sM7sS1Hs7GrPKRlpaGOXPmYM+ePUhNTYVWq0ViYiJqa2uNxyxbtgzLly/HqlWrkJGRAaVSiYSEBFRXc7MqR+WjcMPc0V0AACu3HUN9g07iRETkzE5fqMPGzL/WepDzMat8pKSk4JFHHkGvXr3Qt29frFmzBkVFRcjKygJgmPVYsWIFFi5ciMmTJyM2NhZr165FXV0dkpOTrfIDkG3cNygKEW28UF6txme7C6WOQ0RO7P0dJ6DRiRjauS1u4qyHU3K7kS+urKwEAAQFBQEACgoKUFZWhsTEROMxCoUCI0aMwO7duzFjxowW51Cr1VCr1caPq6qqAAAajQYajeZG4rXQdD5Ln9cZXT5WAoCnR3fG898cxOqdxzElLgwBXu4SJrQffF6ZjmNlHlccr9MX6o2zHnNHdTL5Z3fFsWota42VOecTxFbuHiWKIiZNmoQLFy7g119/BQDs3r0bw4YNQ3FxMcLDw43HPvHEEzh58iR+/vnnFudJSkrC4sWLWzyenJwMb2/v1kQjK9GLwNJ9cpTVC7i1vR4To7j4lIgsa/0JGf4ol6FbgB5zYvg7xpHU1dVh+vTpqKyshL+//zWPbfXMx5NPPon9+/fjt99+a/G5y9+LLYriVd+f/eKLL2L+/PnGj6uqqhAZGYnExMTrhjeXRqNBamoqEhIS4O7Of7Vfy9XGyqtTOWYm5+D3cjcsvv9mhPgpJExpH/i8Mh3HyjyuNl6nLtTh2T9/ByDilamDERcVaPLXutpY3QhrjVXTlQtTtKp8zJ07F99//z127dqFiIgI4+NKpeFmZGVlZQgLCzM+Xl5ejtDQ0CueS6FQQKFo+QLm7u5utSeQNc/tbC4fq7G9wxEXVYjsoov4cFchXrkzVsJ09oXPK9NxrMzjKuP10a6T0OpF3Ny1HW7qHNyqc7jKWFmCpcfKnHOZteBUFEU8+eST+Pbbb7F9+3ZER0c3+3x0dDSUSiVSU1ONjzU0NCAtLQ1Dhw4151uRnRIEAQvG9QAAfJlehKLzdRInIiJnUHS+Dl9nnwYAzLuV+3o4O7PKx5w5c7Bu3TokJyfDz88PZWVlKCsrQ319PQDDC9O8efOwZMkSbNq0CQcPHsQjjzwCb29vTJ8+3So/ANne4E5tcUu3YGj1Ipan5l3/C4iIruO97ceg04u4pVswBnRoI3UcsjKzysfq1atRWVmJkSNHIiwszPjnq6++Mh6zYMECzJs3D7Nnz0Z8fDyKi4uxdetW+Pn5WTw8SWfB2O4AgM37SnCkzPTrfERElys8V4tv9xYDAJ7hvh4uwezLLlf688gjjxiPEQQBSUlJKC0thUqlQlpaGmJjuS7A2cS2D8BtfcIgisBbP3P2g4ha773tx6HTixjZPRj9ozjr4Qp4bxdqtWcTukEuE/DL4XJkFlZIHYeIHFDhuVp8l2OY9eBaD9fB8kGt1inYF1PjDe92WpaSh1ZuGUNELmxl41qP0T1C0C8yUOo4ZCMsH3RDnhrTFR5uMqQXVmDn0bNSxyEiB5J/tgbfNa71eHoM13q4EpYPuiFhAV54eEgHAMCbKXnQ6zn7QUSmeW/7cehFYEyPEPTlrIdLYfmgGzZ7ZBf4KdyQW1qFHw6USh2HiBzAibM12My1Hi6L5YNuWBsfDzx+SycAwNtb86DR8X4MRHRtK7cdg14Ebu0Zit4RAVLHIRtj+SCLeGx4NNr6eKDwfB02Zp6WOg4R2bHj5TX4fl8JAGAe9/VwSSwfZBE+Cjc8OboLAODdbUeh0ugkTkRE9mrltmMQRSAhJhSx7Tnr4YpYPshipt8UhfaBXjhTpcba3YVSxyEiO3TsTDW27Oesh6tj+SCLUbjJ8UyCYeHYBztPoLJeI3EiIrI37zbOeoztFYpe4Zz1cFUsH2RRd/Vvj64hvqis1+CTXflSxyEiO3L0TLXxHXF8h4trY/kgi5LLBDzXeNO5//xWgLPVaokTEZG9aJr1GB+rRM8wf6njkIRYPsjiEmNC0S8yEPUaHVZtPyZ1HCKSkE4v4o8T57F653H8uN8w6/E013q4PJYPsjhBELBgnGH2Izm9CKcq6iRORERSSDlYiuFLt+O+T/ZgaUoeRACebjIUnquVOhpJjOWDrGJo53a4uWs7aHQi3kk9KnUcIrKxlIOlmLUuG6WVqmaPq7R6zFqXjZSD3A3ZlbF8kNU837j2Y1NOMfLKqiVOQ0S2otOLWLwlF9e609PiLbnQ8V5QLovlg6ymT0QgJvRWQhSBt7bmSR2HiGwkvaCixYzHpUQApZUqpBdU2C4U2RWWD7Kq+QndIROA1NwzyDp5Qeo4RGQD5dVXLx6tOY6cD8sHWVWXEF9MGRAJAFiWcgSiyGlWImcX4udp0ePI+bB8kNU9fWtXeLjJ8GdBBXYdOyd1HCKyskHRQfD3dLvq5wUAYQGeGBQdZLtQZFdYPsjqwgO98NDgDgCAN38+Aj0XmRE5tcOlVahruPLNJYXG/100MQZymXDFY8j5sXyQTcwe1QW+CjccLK7Cj3yLHZHTqlZpMCc5G1q9iL4RAVD6K5p9XhngidUPxGFcbJhECckeXH1ejMiCgnw88PjNnfDOL0exfOtRjOulhJuc3ZfImYiiiBe/PYCT5+vQPtALax8dBD9Pd6QXVKC8WoUQP8OlFs54EH/7k808dnM02vp4IP9cLb7OOi11HCKysOT0IvxvfyncZALem94fgd4ekMsEDOncFpP6tceQzm1ZPAgAywfZkK/CDbNHdQEArPjlGNLyzmJzTjH+OHGemw0RObjckios3pILAFgwrjviotpInIjsGS+7kE3df1MU3t9xHGVVKjy8Jt34eFiAJxZNjOF1YCIHVKPW4snkbDRo9RjdIwR/H95J6khk5zjzQTa1M68cFbUNLR4vq1Txfg9EDkgURfxz0wHkn6tFWIAn3p7SFzJeWqHrYPkgm2m638OVNF104f0eiBzLxszT+C6nBHKZgPfu6482Ph5SRyIHwPJBNsP7PRA5l6NnqvHy9wcBAM8mdkN8R24aRqZh+SCb4f0eiJxHXYMWs7/Ihkqjx4huwZh5S2epI5EDYfkgm+H9Hoicx8ubD+F4eQ1C/RVYPpXrPMg8LB9kM4OigxAW4Ilr/Yri/R6I7N83WafxddZpyATg3Xv7o62v4vpfRHQJlg+yGblMwKKJMQBw1QJy36BIbkJEZMeOl1fjn98Z1nnMu7UbBndqK3EickQsH2RT42LDsPqBOCgDml9a8XQ3PBXX7j6J4ov1UkQjouuob9Bhzhd7Ua/RYViXtpjTuGkgkbm4yRjZ3LjYMCTEKJvd76F3+wBM/egP5JZWYcbnmfh65lB4usuljkpEl1i85RDyzlSjna8CK6b15ywltRpnPkgSl9/vwdfTDR8/NABBPh44WFyFF789AFHkfh9E9mJzTjHWZ5yCIADv3tsPwX5c50Gtx/JBdiOijTdWTTf8a2rT3mJ8+nuh1JGICED+2Rq89O0BAMDc0V0xrEs7iRORo2P5ILsytHM7LJzQEwCw5MfD+P34OYkTEbk2lUaHOcl7Udugw+BOQXh6TFepI5ETYPkgu/O3YR0xOa49dHoRTyZn41RFndSRiFzWqz/k4nBpFdr6eODde7nOgyyD5YPsjiAIWHJXb/SJCMCFOg2e+DwL9Q06qWMRuZz/7S/Buj1FAIB3pvVDqD83ACTLYPkgu+TpLseHDwxAO18PHC6twoJv9nMBKpENnTxfi398Y1jnMXtkZ9zSLVjiRORMWD7IboUHeuGD+wfATSZgy74SfLwrX+pIRC5BrdVhTnI2atRaDOzYBvMTukkdiZwMywfZtUHRQcZdUZemHEHa0bMSJyJyfq//eAQHi6vQxtsdK+/rDzc5XyrIsviMIrv3wOAOmBYfCb0IzE3OxsnztVJHInJaKQdL8dnuQgDA8qn9EBbgJW0gckosH2T3BEHAv+7shf5RgahSafHEf7NQq9ZKHYvI6ZyqqMPzX+8HAMy4pRNG9QiROBE5K5YPcggKN8MC1GA/BfLOVOO5jfu4AJXIghq0ejyZnI1qlRZxUYF4bmx3qSORE2P5IIcR6u+JDx+Ig7tcwE8Hy/DBzhNSRyJyGstSjmDf6UoEeLnjvelxcOc6D7IiPrvIoQzoEIR/TYoFALy1NQ/bj5yROBGR4/sl9wz+/VsBAOCtKX3RPpDrPMi6WD7I4dw3KAr33xQFUQSeXp+D/LM1UkcicljFF+vx7MZ9AIBHh0UjISZU4kTkClg+yCEtmtgL8R3aoFqlxROfZ6FapZE6EpHD0ej0mJucjcp6DfpGBOAf43tIHYlcBMsHOSQPNxk+eCAOof4KHC+vwfwN+6DXcwEqkTne2pqH7KKL8PN0w6rpcfBw40sC2QafaeSwQvw88dGD8fCQy5CaewYrtx+TOhKRw9hxpBwfpRl2DX7znj6IDPKWOBG5EpYPcmj9IgPx6l2GBagrfjmGrYfKJE5EZP9KK+sxf0MOAODhIR0wLjZM2kDkclg+yOFNjY/EI0M7AgDmb9iH4+XV0gYismNanR5PfbkXF+o0iG3vj5du6yl1JHJBLB/kFBbe1hM3RQehRq3F4//NQmU9F6ASXck7vxxFRuEF+CrcsOq+OCjc5FJHIhfE8kFOwV0uw/v3xyE8wBMF52oxb/1e6LgAlaiZXUfPGjfne+Pu3ujYzkfiROSqWD7IabTzVeCjB+OhcJNhR95ZvJN6VOpIRHajvEqFZ77KgSgC998Uhdv7hEsdiVwYywc5ld4RAXjj7t4AgFU7juOnA6USJyKSnk4v4un1OThf24AeSj/83+0xUkciF8fyQU7nrv4ReGx4NADg2Y37kFfGBajk2lZuO4Y/8s/D20OO9++Pg6c713mQtFg+yCm9OL4HhnZui7oGHR7/byYu1jVIHYlIEruPnzPugbPkrt7oHOwrcSIilg9yUm5yGVZNj0NEGy8UVdRh7pdcgEqu52y1Gk83rvOYFh+JO/u3lzoSEQCWD3JiQT4e+PjBeHi6y/DrsXNY9vMRqSMR2YxOL+KZr3JwtlqN7qF+SLqjl9SRiIzMLh+7du3CxIkTER4eDkEQ8N133zX7vCiKSEpKQnh4OLy8vDBy5EgcOnTIUnmJzBIT7o9l9/QFAHyUlo8t+0okTkRkGx/sOI7fjp+Dl7sc79/fH14eXOdB9sPs8lFbW4u+ffti1apVV/z8smXLsHz5cqxatQoZGRlQKpVISEhAdTUX/ZE07ugbjhkjOgEAnv96H3JLqiRORGRde/LP451fDG81f+XOWHQJ8ZM4EVFzZpeP8ePH49VXX8XkyZNbfE4URaxYsQILFy7E5MmTERsbi7Vr16Kurg7JyckWCUzUGgvG9sDNXdtBpdHjic8zUVHLBajknM7XqPH0+r3Qi8DdcRG4Z0CE1JGIWnCz5MkKCgpQVlaGxMRE42MKhQIjRozA7t27MWPGjBZfo1aroVarjR9XVRn+VarRaKDRWHaL7KbzWfq8zsgZx2r5Pb1x90d7UFRRjzlfZOHTh+LgJr/xZU/OOFbWwrEyj7njpdeLeOarvThTpUandj54+bZuLjPWfG6ZzlpjZc75BFEUW/0WAEEQsGnTJtx5550AgN27d2PYsGEoLi5GePhfu+c98cQTOHnyJH7++ecW50hKSsLixYtbPJ6cnAxvb97imSyrpA5454AcDXoBI8P0uKujXupIRBbzS7GALUVyuAsi5vfWIZy7p5MN1dXVYfr06aisrIS/v/81j7XozEcTQRCafSyKYovHmrz44ouYP3++8eOqqipERkYiMTHxuuHNpdFokJqaioSEBLi7u1v03M7GmccqoucZzF2/DztLZbh9aB9M6ndj20w781hZGsfKPOaMV9bJC/jxz0wAIpLu6IWp8a51uYXPLdNZa6yarlyYwqLlQ6lUAgDKysoQFhZmfLy8vByhoaFX/BqFQgGFQtHicXd3d6s9gax5bmfjjGM1sV8E8s7UYtWO41i4ORfdwwLROyLghs/rjGNlLRwr81xvvC7UNuCZjQeg04uY1C8c0wd3vOo/+Jwdn1ums/RYmXMui+7zER0dDaVSidTUVONjDQ0NSEtLw9ChQy35rYhuyDMJ3TC6RwjUWj1mfJ6JczXq638RkR0SRRHPbdyH0koVOrXzwWt39XbZ4kGOw+zyUVNTg5ycHOTk5AAwLDLNyclBUVERBEHAvHnzsGTJEmzatAkHDx7EI488Am9vb0yfPt3S2YlaTS4T8M60fujUzgcllSrM/iIbGh3Xf5Dj+fevBdh2pBweboZdfX0VVrmaTmRRZpePzMxM9O/fH/379wcAzJ8/H/3798fLL78MAFiwYAHmzZuH2bNnIz4+HsXFxdi6dSv8/Pg+c7IvAV7u+PihAfBVuCG9oAKv/XBY6khEZskuuoClKYade1++PQYx4ZZdJ0dkLWZX5JEjR+Jab5ARBAFJSUlISkq6kVxENtElxA/Lp/bFE59n4bPdhYgJ98fU+EipYxFdV2WdBnOT90KrF3FbnzDcf1OU1JGITMZ7u5DLS+ylxLxbuwIA/rnpIPYWXZA4EdG1iaKI57/eh+KL9ejQ1htvTOY6D3IsLB9EAJ4a3RUJMaFo0Okxc10WyqtVUkciuqrPdhdia+4ZeMhleH96HPw8+e4OciwsH0QAZDIBy6f2RedgH5ypUmP2umw0aLkAlezP/tMXseRHw/qklyb0QGz7G3+bOJGtsXwQNfLzdMcnD8XDT+GGzJMXsHgL78ZM9qVKpcGTyXuh0YkY10uJh4d2lDoSUauwfBBdolOwL969rx8EAfjizyIk/1kkdSQiAIZ1Hv/4Zj+KKuoQ0cYLS+/pw3Ue5LBYPoguM7pHKJ5L7A4AWPT9QWSdrJA4ERGwbs9J/HigDO5yAaumxyHAi+s8yHGxfBBdweyRnTGhtxIanYiZ67JxpooLUMm2dHoR6QWG4rsh8zRe+Z9hnccL43qgX2SghMmIbhzLB9EVCIKAN+/pi+6hfjhbrcaMz7Og1uqkjkUuIuVgKYYv3Y5H12ZApQX+ufkQGnR69IkIwGPDo6WOR3TDWD6IrsJH4YaPHxoAf0835Jy6iP/77uA1N9gjsoSUg6WYtS4bpZUqiCLwVb4MIgxrO/afrsTPh8okTkh041g+iK6hQ1sfvDc9DjLBMPW9bs9JqSORE9PpRSzekoumiqsXgezzMqDxEQHA4i250OlZgsmxsXwQXceIbsF4YVwPAIZf/H/mn5c4ETmr9IIKlFb+tb5I29gx5I1vahEBlFaqjGtBiByVw97+UKfTQaPRmPU1Go0Gbm5uUKlU0Ol4/f5arDVW7u7ukMvlFjufrTxxSyccLKnCln0lmP1FNrbMHY7wQC+pY5GTabmzroCYQD2OVwnQidc6jsixOFz5EEURZWVluHjxYqu+VqlU4tSpU3x//HVYc6wCAwOhVCod6u9AEAQsvbs3jpfX4HBpFWZ8noWNM4fA093xihTZr2qVttnHckHE37rp8c+s5s+zED9PW8YisjiHKx9NxSMkJATe3t5mvYDp9XrU1NTA19cXMhmvOF2LNcZKFEXU1dWhvLwcABAWFmaR89qKt4cbPn5wAO5Y9RsOFFfipU0H8PaUvlLHIieg1enx/o4TWLn9WLPH3WSAxyW9QwCgDPDEoOgg2wYksjCHKh86nc5YPNq2bWv21+v1ejQ0NMDT05Pl4zqsNVZeXoZLFeXl5QgJCXG4SzCRQd5YNT0OD32ajm+zixET5o+eoT4ADNfrB3cJgVzmODM6JL38szV4ZsM+7Dt1EQAwoEMgsk5exOXPoqaPF02M4XOMHJ5DvQI3rfHw9vaWOAndiKa/P3PX7NiLYV3a4cXxhgWor/5wGI98lgEAeHRtBoYv3Y6Ug6VSxiMHIYoi1u05idtW/oZ9py7Cz9MN797bD1/PHIoPH4iDMqD5pRVlgCdWPxCHcbGONWNIdCUONfPRxJHWClBLzvD31z7wrxcGjR4437j+r6xShVnrsvkiQddUXqXCgm/2Y2feWQDA0M5t8daUvsZFzONiw5AQo8Se4+U4d3gPPn14IGfVyKk41MwHkT3Q6UX8q3GrawMB/8mTQxRh3J+BezHQ1fx0oBRjV+zCzryz8HCT4eXbY7DusZtavHtKLhOMazsGRQexeJBTcciZDyIpXb4XAyCiuE5AU/W4dC+GIZ3NX5tEzqlKpUHS94fwbXYxAKBXuD9WTOuHrqF+Eicjsj3OfNjII488AkEQIAgC3N3dERoaioSEBHz66afQ6/Umn+ezzz5DYGCg9YLSdV2+x4K7DAjzFoHLlghyLwZqsif/PMav+BXfZhdDJgBzRnXGptnDWDzIZbls+dDpRfxx4jw25xTjjxPnbTJFPm7cOJSWlqKwsBA//fQTRo0ahaeffhq33347tFrt9U9AduHyPRZkAvB8Hx3kQvPn0Hd7i3G+Rm3LaGRnVBodlvx4GPd9sgfFF+sRFeSNDTOG4PmxPeDh5rK/folcs3ykHCzD8KXbcd8ne/D0+hzc98kem7xLQaFQQKlUon379oiLi8NLL72EzZs346effsJnn30GAFi+fDl69+4NHx8fREZGYvbs2aipqQEA7Ny5E3/7299QWVlpnEVJSkoCAKxbtw7x8fHw8/ODUqnE9OnTjftpkGUNig5CWIBns3kOuWDYk+FSO/LOYtRbO/H5H4Vc/+GCDpdW4c73f8fHu/IhisC9AyPx49M3I74j9+ggcrnysS3vPOYk773smv1f71Kw9dskR48ejb59++Lbb78FAMhkMqxcuRIHDx7E2rVrsX37dixYsAAAMHToUKxYsQL+/v4oLS1FaWkpnnvuOQBAQ0MDXnnlFezbtw/fffcdCgoK8Mgjj9j0Z3EVcpmARRNjAFx+ocXwsQDgucRu6BnmjyqVFv+3+RDuWPUbsk5esHVUkoBOL+LDtBO4Y9VvOFJWjXa+HvjkoXi8cXcf+Cq4zI4IcLEFpzq9iGW/5ONK/wZtumK/eEsuEmKUNl1Z3qNHD+zfvx8AMG/ePOPj0dHReOWVVzBr1ix88MEH8PDwQEBAAARBgFKpbHaORx991PjfnTp1wsqVKzFo0CDjLqVkWeNiw7D6gTgs3pKLipp64+PKAE8smhiDcbFhmDmiM774swhvbc3DoZIq3L16N+4ZEIF/jO+Bdr4KCdOTtZyqqMOzG/YhvdBw47dbe4bijbt78++b6DIuVT4yCitwprrhqp+X6l0Koiga977YsWMHlixZgtzcXFRVVUGr1UKlUqG2thY+Pj5XPcfevXuRlJSEnJwcVFRUGBexFhUVISYmxiY/h6u53l4MbnIZHh7aEbf1CcPSn45gY9ZpfJ11Gj8fKsOzCd3wwOAOcJO73OSjUxJFEV9nncbiLbmoUWvh4yHHoom9MCU+win2tSGyNJf6zVdebdriP1u/S+Hw4cOIjo7GyZMnMWHCBMTGxuKbb75BVlYW3n//fQDX3g20trYWiYmJ8PX1xbp165CRkYFNmzYBMFyOIesxZS+Gdr4KvDmlL76ZNRS9wv1RrdIiaUsubn/vN2QU8tboju58jRoz12Xh+a/3o0atRXyHNvjp6VswdWAkiwfRVbjUzEeIn2lTn7a8Y+T27dtx4MABPPPMM8jMzIRWq8Xbb79tvJ/Khg0bmh3v4eHR4hb3R44cwblz5/DGG28gMjISAJCZmWmbH4BMNqBDG3z/5HAkpxfhrZ/zcKSsGlM+/AOT+7fHPyb04J1KHdC2w2fwwjcHcK5GDXe5gGcSumHGLZ25IRjRdbjUzMfAjkEI9fNosUiwiQAgzIp3jFSr1SgrK0NxcTGys7OxZMkSTJo0CbfffjseeughdO7cGVqtFu+99x7y8/Px+eef48MPP2x2jo4dO6Kmpgbbtm3DuXPnUFdXh6ioKHh4eBi/7vvvv8crr7xilZ+BboxcJuDBwR2w47mRuHdgJAQB+HZvMca8lYb//FYArc70PV9IOrVqLV789gAeW5uJczVqdAv1xXdzhmH2yC4sHkQmcKnyIZcJWHBrJwBXfpcCYN07RqakpCAsLAwdO3bEuHHjsGPHDqxcuRKbN2+GXC5Hv379sHz5cixduhSxsbH44osv8Prrrzc7x9ChQzFz5kxMmzYNwcHBWLZsGYKDg/HZZ59h48aNiImJwRtvvIG33nrLKj8DWUaQjwfeuLsPNs0ehj4RAahWa/HK/3Jx28rf8Gf+eanj0TVknbyACSt/xZfpRQCAvw+PxvdPDkev8ACJkxE5DkEURbvagKCqqgoBAQGorKyEv79/s8+pVCoUFBQgOjoanp7mT1Hr9XpUVVVhd1EdXvnhcLO324Zd8i4F+mus/P39jZeALOVG/x7tjUajwY8//ogJEybA3d3d7K/X6UV8lXEKy34+got1hrU9k/qF46UJPRHq7/jjc6kbHSspaXR6rNx2DO/vOA69CIQHeOKtKX0xtEs7631PBx4vW+NYmc5aY3Wt1+/LudSajybjYpUYGxuG9IIKlFerEOLnyRs3kWTkMgHTb4rC+Fgl3tyahy/Ti7A5pwS/5J7BvFu74ZFhHeHOd8VI6nh5NZ75ah8OFFcCAO7q3x5Jd/RCgBdf5IhawyXLB2D4hc+bfpE9aePjgSV39ca9AyPxf5sPYd+pi3jtx8PYkHkK/5oUy+erBPR6EWv/KMQbPx2BWqtHgJc7XrsrFrf3CZc6GpFD4z+niOxMn4hAbJo1FEvv7o0gHw8cK6/BfZ/swdwv96Kskjers5XSyno89Gk6Fm/JhVqrx81d22HrM7eweBBZgMvOfBDZM5lMwLSBURjbS4m3tx7Fuj9PYsu+Emw7fAZPjemKR4dF88ZkVvT9vhL8c9MBVKm08HSX4aUJPfHg4A7ct4PIQvjbi8iOBXp74JU7Y7HlyeGIiwpEXYMOb/x0BOPf3YXfj5+TOp7TqazT4Kkv9+KpL/eiSqVF34gA/PDUzXhoSEcWDyILYvkgcgCx7QPw9cyhePOePmjr44ETZ2tx/7//xJwvslFysf76J6Dr+u3YOYxdsQvf7yuBXCbg6TFd8fWsoegczHsjEVkaL7sQOQiZTMCU+Egk9lLindSj+O8fhfjhQCm2HynH3DFd8PfhnXgpphVUGsNs0me7CwEA0e18sHxqX/SPaiNtMCInxt9URA4mwMsdSXf0wpa5wxHfoQ3qNTosS8nDuBW7sOvoWanjOZQDpytx+3u/GYvHA4Oj8MNTw1k8iKyMMx9EDqpXeAA2zhyCb7OL8fpPR5B/rhYPfZqOcb2U+L+JMWgf6CV1RLul1enxYdoJrPjlGLR6EcF+Ciy7pw9GdQ+ROhqRS+DMhx0RBAHffffdVT9fWFgIQRCQk5Njs0zX0rFjR6xYsULqGC5NEATcPSAC258bgUeHRUMuE5ByqAxj3t6JVduPQa3VXf8kLqbwXC2mfvQH3tp6FFq9iPGxSvw87xYWDyIbYvmwkUceeQSCIEAQBLi5uSEqKgqzZs3ChQsXjMeUlpZi/PjxEqZsafTo0cbcl/7RarXIyMjAE088YTz2euWJrMff0x0vT4zBD08Nx6COQVBp9Hhr61GMfWcXduSVSx3PLoiiiC/TizBh5a/ILroIP4Ublk/tiw/uj0OQj4fU8YhcCi+72NC4ceOwZs0aaLVa5Obm4tFHH8XFixfx5ZdfAgCUSqXECa/s8ccfx7/+9a9mj7m5uSE4OFiiRHQ1PZT++GrGYGzOKcFrPx5G4fk6/G1NBhJiQvHy7TGIDPKWOqLV6fRii1snnK9V48VvDmDbEUMRuyk6CG9P7YuINs4/HkT2iDMfNqRQKKBUKhEREYHExERMmzYNW7duNX7+8pmD9PR09O/fH56enoiPj8fevXtbnPP7779H165d4eXlhVGjRmHt2rUQBAEXL140HrN7927ccsst8PLyQmRkJJ566inU1taanNvb2xtKpbLZH6D5ZZeOHTsCAO666y4IgmD8mGxPEATc2b89tj87An8fbrgUk5p7BrcuT8O7vxyDSuO8l2JSDpZi+NLtuO+TPXh6fQ7u+2QP4l9Nxei30rDtSDk85DIsnNATXz4+mMWDSEIOXz5EUURdg9bkP/UNOrOOv9afG7khcH5+PlJSUq56R8Ha2lrcfvvt6N69O7KyspCUlITnnnuu2TGFhYW45557cOeddyInJwczZszAwoULmx1z4MABjB07FpMnT8b+/fvx1Vdf4bfffsOTTz7Z6uxXkpGRAQBYs2YNSktLjR+TdPw83fHP22Pw09M3Y3CnIKi1erzzy1EkvrML2w6faXasTi/ijxPnsTmnGH+cOA+d3q5udm2SlIOlmLUuu9ndqgHgQp0GNWot2gd64fu5w/D4LZ0g400kiSTl8Jdd6jU6xLz8syTfO/dfY+HtYfoQ/u9//4Ovry90Oh1UKsMvyOXLl1/x2C+++AI6nQ6ffvopvL290atXL5w+fRqzZs0yHvPhhx+ie/fuePPNNwEA3bt3x8GDB/Haa68Zj3nzzTcxffp0zJs3DwDQtWtXrFy5EiNGjMDq1atNuqX9Bx98gH//+9/Gj2fMmIG333672TFNl2ACAwPt9vKRq+oW6ocvHx+MLftL8doPuSiqqMNjazMxpkcIFk3shdzSSizektvsRTsswBOLJsZgXGyYhMlNp9OLWLwlF9eqTDq9Hl1D/GyWiYiuzuHLhyMZNWoUVq9ejbq6Ovz73//G0aNHMXfu3Csee/jwYfTt2xfe3n9NDQ8ZMqTZMXl5eRg4cGCzxwYNGtTs46ysLBw/fhxffPGF8TFRFKHX61FQUICePXteN/f999/fbEYlMDDwul9D9kUQBNzRNxyje4TgvW3H8J/fCrDtSDnSjp6F9gqzHGWVKsxal43VD8RZvYCIooh6jQ41ai3q1I3/26BDbYMWtc0e06JGbZi5rFXrUKvWGo85W9PQYsajxc9UpUZ6QQXvDkxkBxy+fHi5y5H7r7EmHavX61FdVQ0/fz/IZDd+xcnLXW7W8T4+PujSpQsAYOXKlRg1ahQWL16MV155pcWxplzSEUWxxf0mLv86vV6PGTNm4Kmnnmrx9VFRUSblDggIMOYmx+arcMOLE3piSnwEXt58CLtPnL/icSIAAcDiLblIiFFCfslligat/pIXfkNJuLQgGD7XWA4ay0JToahWaVB6Vo53jv6GugadsWTcwBVMs5RX867ARPbA4cuHIAgmX/rQ6/XQesjh7eFmkfJxoxYtWoTx48dj1qxZCA9vfpvumJgYfP7556ivr4eXl2GzqD179jQ7pkePHvjxxx+bPZaZmdns47i4OBw6dMgm5cHd3R06nfMuZnQmXUL88OSoLlctH4ChgJRWqnDLsh0AYJxl0OhutCkIQG1dy0cFwMfDDd4ecvgq3OCtkMPHww0+ikse83CDr0IOb4UbfDzkjZ9zQ1FFHZb8ePi63znE7/qXGYnI+hy+fDiykSNHolevXliyZAlWrVrV7HPTp0/HwoUL8dhjj+Gf//wnCgsL8dZbbzU7ZsaMGVi+fDleeOEFPPbYY8jJycFnn30GAMYZkRdeeAGDBw/GnDlz8Pjjj8PHxweHDx9Gamoq3nvvPYv+PB07dsS2bdswbNgwKBQKtGnDLart2dkatUnHFV/lxnUKN5lZJcHHww2ebsChfXsxYthgBPgojF/no5DDy13e6jvH6vQi1vxegLJK1RXXfQgAlAGGt90SkfRYPiQ2f/58/O1vf8MLL7zQ7HFfX19s2bIFM2fORP/+/RETE4OlS5fi7rvvNh4THR2Nr7/+Gs8++yzeffddDBkyBAsXLsSsWbOgUCgAAH369EFaWhoWLlyIm2++GaIoonPnzpg2bZrFf5a3334b8+fPxyeffIL27dujsLDQ4t+DLMfUWYB/3tYTAzsGGUuCt4ehULjJzZ891Gg0EItEDOzY5qrv9GoNuUzAookxmLUuGwLQrIA01ZlFE2OaXT4iIumwfNhI04zE5aZPn47p06cDaLleY/DgwS22Ur/8mDvuuAN33HGH8ePXXnsNERERzd7FMnDgwGb7iZhj+/btV71EdXm5mDhxIiZOnNiq70O2Nyg6CGEBntedLfhb47bt9m5cbBhWPxDX4p07Sgd75w6RK2D5cHAffPABBg4ciLZt2+L333/Hm2++afE9PMg5OeNswbjYMCTEKFvscOpIPwORK2D5cHDHjh3Dq6++ioqKCkRFReHZZ5/Fiy++aNLX/vrrr9e8l8zp06ctFZPslDPOFshlAt9OS2TnWD4c3DvvvIN33nmnVV8bHx9/1Tvk6vX6G0hFjoSzBURkaywfLszLy+uqb8HV6/WoqqqycSKSCmcLiMiWpN/sgoiIiFyKQ5YPXhJwbPz7IyJybQ512cXDwwMymQwlJSUIDg6Gh4eHWZsS6fV6NDQ0QKVS2cUOp/bMGmMliiIaGhpw9uxZyGQyeHh4WOS8RETkWByqfMhkMkRHR6O0tBQlJSVmf70oisbtylu7k6KrsOZYeXt7IyoqigWQiMhFOVT5AAyzH1FRUdBqtWbfR0Sj0WDXrl245ZZbLLq7ojOy1ljJ5XK4ubmx/BERuTCHKx+A4b4l7u7uZr8oyuVyaLVaeHp6snxcB8eKiIishfPeREREZFMsH0RERGRTLB9ERERkU3a35qPprq3W2F1To9Ggrq4OVVVVXMdwHRwr03GsTMexMg/Hy3QcK9NZa6yaXrcvv/v6ldhd+aiurgYAREZGSpyEiIiIzFVdXY2AgIBrHiOIplQUG9Lr9SgpKYGfn5/F345ZVVWFyMhInDp1Cv7+/hY9t7PhWJmOY2U6jpV5OF6m41iZzlpjJYoiqqurER4eft19nOxu5kMmkyEiIsKq38Pf359PThNxrEzHsTIdx8o8HC/TcaxMZ42xut6MRxMuOCUiIiKbYvkgIiIim3Kp8qFQKLBo0SIoFAqpo9g9jpXpOFam41iZh+NlOo6V6exhrOxuwSkRERE5N5ea+SAiIiLpsXwQERGRTbF8EBERkU2xfBAREZFNuWz5uOOOOxAVFQVPT0+EhYXhwQcfRElJidSx7E5hYSEee+wxREdHw8vLC507d8aiRYvQ0NAgdTS79dprr2Ho0KHw9vZGYGCg1HHsygcffIDo6Gh4enpiwIAB+PXXX6WOZJd27dqFiRMnIjw8HIIg4LvvvpM6kl16/fXXMXDgQPj5+SEkJAR33nkn8vLypI5lt1avXo0+ffoYNxcbMmQIfvrpJ0myuGz5GDVqFDZs2IC8vDx88803OHHiBO655x6pY9mdI0eOQK/X46OPPsKhQ4fwzjvv4MMPP8RLL70kdTS71dDQgClTpmDWrFlSR7ErX331FebNm4eFCxdi7969uPnmmzF+/HgUFRVJHc3u1NbWom/fvli1apXUUexaWloa5syZgz179iA1NRVarRaJiYmora2VOppdioiIwBtvvIHMzExkZmZi9OjRmDRpEg4dOmT7MCKJoiiKmzdvFgVBEBsaGqSOYveWLVsmRkdHSx3D7q1Zs0YMCAiQOobdGDRokDhz5sxmj/Xo0UP8xz/+IVEixwBA3LRpk9QxHEJ5ebkIQExLS5M6isNo06aN+O9//9vm39dlZz4uVVFRgS+++AJDhw7lrZhNUFlZiaCgIKljkANpaGhAVlYWEhMTmz2emJiI3bt3S5SKnE1lZSUA8PeTCXQ6HdavX4/a2loMGTLE5t/fpcvHCy+8AB8fH7Rt2xZFRUXYvHmz1JHs3okTJ/Dee+9h5syZUkchB3Lu3DnodDqEhoY2ezw0NBRlZWUSpSJnIooi5s+fj+HDhyM2NlbqOHbrwIED8PX1hUKhwMyZM7Fp0ybExMTYPIdTlY+kpCQIgnDNP5mZmcbjn3/+eezduxdbt26FXC7HQw89BNFFNnw1d6wAoKSkBOPGjcOUKVPw97//XaLk0mjNeFFLgiA0+1gUxRaPEbXGk08+if379+PLL7+UOopd6969O3JycrBnzx7MmjULDz/8MHJzc22ew83m39GKnnzySdx7773XPKZjx47G/27Xrh3atWuHbt26oWfPnoiMjMSePXskmYKyNXPHqqSkBKNGjcKQIUPw8ccfWzmd/TF3vKi5du3aQS6Xt5jlKC8vbzEbQmSuuXPn4vvvv8euXbsQEREhdRy75uHhgS5dugAA4uPjkZGRgXfffRcfffSRTXM4VfloKhOt0TTjoVarLRnJbpkzVsXFxRg1ahQGDBiANWvWQCZzqgkzk9zIc4sMv/AGDBiA1NRU3HXXXcbHU1NTMWnSJAmTkSMTRRFz587Fpk2bsHPnTkRHR0sdyeGIoijJ655TlQ9TpaenIz09HcOHD0ebNm2Qn5+Pl19+GZ07d3aJWQ9zlJSUYOTIkYiKisJbb72Fs2fPGj+nVColTGa/ioqKUFFRgaKiIuh0OuTk5AAAunTpAl9fX2nDSWj+/Pl48MEHER8fb5xBKyoq4vqhK6ipqcHx48eNHxcUFCAnJwdBQUGIioqSMJl9mTNnDpKTk7F582b4+fkZZ9YCAgLg5eUlcTr789JLL2H8+PGIjIxEdXU11q9fj507dyIlJcX2YWz+/ho7sH//fnHUqFFiUFCQqFAoxI4dO4ozZ84UT58+LXU0u7NmzRoRwBX/0JU9/PDDVxyvHTt2SB1Ncu+//77YoUMH0cPDQ4yLi+NbIq9ix44dV3wOPfzww1JHsytX+920Zs0aqaPZpUcffdT4/7/g4GBxzJgx4tatWyXJIoiii6ywJCIiIrvgehfviYiISFIsH0RERGRTLB9ERERkUywfREREZFMsH0RERGRTLB9ERERkUywfREREZFMsH0RERGRTLB9ERERkUywfREREZFMsH0RERGRTLB9ERERkU/8Pi+CA51rwAuoAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ "from sklearn import linear_model\n", "np.random.seed(2018)\n", - "n = 100\n", + "n = 10\n", "d = 2\n", "Lambda = 0.01\n", "\n", "# Make data set.\n", "x = np.linspace(-3, 3, n)\n", - "y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n)\n", + "y = 2.0 + 0.5*x + 5.0*(x**2)#+ np.random.randn(n)\n", "\n", "# Design matrix X does not include the intercept. \n", "X = np.zeros((n, d))\n", @@ -594,25 +607,34 @@ "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", "\n", "# Scale data by subtracting mean value using scikit-learn\n", - "from sklearn.preprocessing import StandardScaler\n", - "scaler = StandardScaler()\n", + "scaler = StandardScaler(with_std=False)\n", + "scaler.fit(X_train)\n", + "X_train_mean = np.mean(X_train,axis=0)\n", + "#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)\n", + "y_scaler = np.mean(y_train)\n", "X_train_scaled = scaler.transform(X_train)\n", + "X_test_scaled = scaler.transform(X_test)\n", + "y_train_scaled = y_train - np.mean(y_train)\n", + "y_test_scaled = y_test-np.mean(y_test)\n", "\n", "#Calculate beta\n", "OLS = LinearRegression()\n", - "OLS.fit(X_train_scaled,y_train)\n", - "ypredictOLS = OLS.predict(X_test)\n", - "RegRidge = linear_model.Ridge(Lambda)\n", - "RegRidge.fit(X_train_scaled,y_train)\n", - "ypredictRidge = RegRidge.predict(X_test)\n", - "print(OLS.coef_)\n", - "print(RegRidge.coef_)\n", - "print(OLS.intercept_)\n", - "interceptRidge = RegRidge.intercept_\n", - "print(RegRidge.intercept_)\n", - "#predict value without intercept\n", - "ytilde_test_Ridge = X_test @ RegRidge.coef_+ RegRidge.intercept_\n", - "ytilde_test_OLS = X_test @ OLS.coef_+ OLS.intercept_\n", + "betaOLS=OLS.fit(X_train_scaled,y_train_scaled)\n", + "ypredictOLS = OLS.predict(X_test_scaled)\n", + "linear_model.Ridge(Lambda)\n", + "RegRidge.fit(X_train_scaled,y_train_scaled)\n", + "ypredictRidge = RegRidge.predict(X_test_scaled)\n", + "betaOLS = OLS.coef_\n", + "betaRidge = RegRidge.coef_\n", + "print(betaOLS)\n", + "print(betaRidge)\n", + "interceptOLS = np.mean(y_train) - X_train_mean @ betaOLS\n", + "interceptRidge = y_scaler - X_train_mean @ betaRidge\n", + "print(interceptOLS)\n", + "print(interceptRidge)\n", + "#predict value \n", + "ytilde_test_Ridge = X_test_scaled @ betaRidge+y_scaler\n", + "ytilde_test_OLS = X_test_scaled @ betaOLS+y_scaler\n", "\n", "#Calculate MSE\n", "print(\" \")\n", @@ -622,7 +644,7 @@ "print(\"test MSE of Ridge\")\n", "print(MSE(y_test,ytilde_test_Ridge))\n", "plt.scatter(x,y,label='Data')\n", - "plt.plot(x, X @ RegRidge.coef_ + RegRidge.intercept_ , label=\"Ridge_Fit\")\n", + "plt.plot(x, X @ betaRidge+interceptRidge, label=\"Ridge_Fit\")\n", "plt.grid()\n", "plt.legend()\n", "plt.show()" @@ -631,7 +653,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b99be971", + "id": "c6130d15", "metadata": {}, "outputs": [], "source": []