lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
be6621e222b048f158ed40e9f1e45d4a2a4b04dd
FYS-STK4155/doc/LectureNotes/exercisesweek35.ipynb
T
lars be6621e222 Finish ex 37
2025-10-08 08:08:24 +02:00

1177 lines
88 KiB
Plaintext
Raw Blame History

{
"cells": [
{
"cell_type": "markdown",
"id": "b4005770",
"metadata": {},
"source": [
"# Exercises week 35\n",
"\n",
"## Deriving and Implementing Ordinary Least Squares"
]
},
{
"cell_type": "markdown",
"id": "2ca1b589",
"metadata": {},
"source": [
"This week you will be deriving the analytical expressions for linear regression, building up the model from scratch. This will include taking several derivatives of products of vectors and matrices. Such derivatives are central to the optimization of many machine learning models. Although we will often use automatic differentiation in actual calculations, to be able to have analytical expressions is extremely helpful in case we have simpler derivatives as well as when we analyze various properties (like second derivatives) of the chosen cost functions.\n",
"\n",
"Vectors are always written as boldfaced lower case letters and matrices as upper case boldfaced letters. You will find useful the notes from week 35 on derivatives of vectors and matrices. See also the textbook of Faisal at al, chapter 5 and in particular sections 5.3-5.5 at <https://github.com/CompPhysics/MachineLearning/blob/master/doc/Textbooks/MathMLbook.pdf>"
]
},
{
"cell_type": "markdown",
"id": "41e92bf9",
"metadata": {},
"source": [
"### Learning goals\n",
"\n",
"After completing these exercises, you will know how to\n",
"- Take the derivatives of simple products between vectors and matrices\n",
"- Implement OLS using the analytical expressions\n",
"- Create a feature matrix from a set of data\n",
"- Create a feature matrix for a polynomial model\n",
"- Evaluate the MSE score of various model on training and test data, and comparing their performance\n",
"\n",
"### Deliverables\n",
"\n",
"Complete the following exercises while working in a jupyter notebook. Then, in canvas, include\n",
"- The jupyter notebook with the exercises completed\n",
"- An exported PDF of the notebook (https://code.visualstudio.com/docs/datascience/jupyter-notebooks#_export-your-jupyter-notebook)"
]
},
{
"cell_type": "markdown",
"id": "f7a9209d",
"metadata": {},
"source": [
"## How to take derivatives of Matrix-Vector expressions"
]
},
{
"cell_type": "markdown",
"id": "45f3712e",
"metadata": {},
"source": [
"In these exercises it is always useful to write out with summation indices the various quantities. Take also a look at the weekly slides from week 35 and the various examples included there.\n",
"\n",
"As an example, consider the function\n",
"\n",
"$$\n",
"f(\\boldsymbol{x}) =\\boldsymbol{A}\\boldsymbol{x},\n",
"$$\n",
"\n",
"which reads for a specific component $f_i$ (we define the matrix $\\boldsymbol{A}$ to have dimension $n\\times n$ and the vector $\\boldsymbol{x}$ to have length $n$)\n",
"\n",
"$$\n",
"f_i =\\sum_{j=0}^{n-1}a_{ij}x_j,\n",
"$$\n",
"\n",
"which leads to\n",
"\n",
"$$\n",
"\\frac{\\partial f_i}{\\partial x_j}= a_{ij},\n",
"$$\n",
"\n",
"and written out in terms of the vector $\\boldsymbol{x}$ we have\n",
"\n",
"$$\n",
"\\frac{\\partial f(\\boldsymbol{x})}{\\partial \\boldsymbol{x}}= \\boldsymbol{A}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "5fa8a4e6",
"metadata": {},
"source": [
"## Exercise 1 - Finding the derivative of Matrix-Vector expressions"
]
},
{
"cell_type": "markdown",
"id": "df7a2270",
"metadata": {},
"source": [
"**a)** Consider the expression\n",
"\n",
"$$\n",
"\\frac{\\partial (\\boldsymbol{a}^T\\boldsymbol{x})}{\\partial \\boldsymbol{x}},\n",
"$$\n",
"\n",
"Where $\\boldsymbol{a}$ and $\\boldsymbol{x}$ are column-vectors with length $n$.\n",
"\n",
"What is the *shape* of the expression we are taking the derivative of?\n",
"\n",
"What is the *shape* of the thing we are taking the derivative with respect to?\n",
"\n",
"What is the *shape* of the result of the expression?"
]
},
{
"cell_type": "markdown",
"id": "5fa87526",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"\n",
"- We're taking the derivative of a expression of shape $1\\times1$, as the product of row-vector $a^T$ and column-vector $x$ is of shape $1\\times 1$.\n",
"- We're taking the derivative with respect to a column vector of length n.\n",
"- The resulting shape is a row vector of length n.\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "c0396734",
"metadata": {},
"source": [
"**b)** Show that\n",
"\n",
"$$\n",
"\\frac{\\partial (\\boldsymbol{a}^T\\boldsymbol{x})}{\\partial \\boldsymbol{x}} = \\boldsymbol{a}^T,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "8d3e59ee",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"\n",
"We define $\\alpha \\equiv a^T x = \\sum_j a_j x_j$. Then\n",
"\n",
"$$\n",
"\\frac{\\partial \\alpha}{\\partial x} = \\frac{\\partial}{\\partial x_j} (a_j x_j) = a_j = a^T\n",
"$$\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "dc39d541",
"metadata": {},
"source": [
"**c)** Show that\n",
"\n",
"$$\n",
"\\frac{\\partial (\\boldsymbol{a}^T\\boldsymbol{A}\\boldsymbol{a})}{\\partial \\boldsymbol{a}} = \\boldsymbol{a}^T(\\boldsymbol{A}+\\boldsymbol{A}^T),\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "fc826cc4",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"\n",
"If we define $\\phi \\equiv a^T A a$ which is evidently a scalar (thus $\\phi = \\phi^T$), we can write it as:\n",
"$$\n",
"\\phi = \\sum_{i, j} a_i A_{ij} a_j\n",
"$$\n",
"Using the scalar property we can rewrite the problem as\n",
"$$\n",
"\\frac{\\partial \\phi}{\\partial a} = \\frac{\\partial}{\\partial a} a^T A a = \\frac{\\partial \\phi^T}{\\partial a} = \\frac{\\partial}{\\partial a}a^T A^T a\n",
"$$\n",
"Now evaluating\n",
"$$\n",
"\\frac{\\partial \\phi}{\\partial a_k} = \\sum_{i,j} (\\frac{\\partial a_i}{\\partial a_k} A_{ij} a_j + a_j A_{ij} \\frac{\\partial a_j}{\\partial a_k}) = \\sum_{i,j} (\\delta_{ik} A_{ij} a_j + a_j A_{ij} \\delta_{jk}) = \\sum_{j} A_{kj} a_j + \\sum_{i} A_{ik} a_i = a^T (A + A^T)\n",
"$$\n",
"with the Kronecker-Delta $\\delta_{xy} = \\begin{cases} 0 & x \\neq y \\\\ 1 & x = y \\end{cases}$. As the derivative of two components can be written as $\\frac{\\partial a_x}{\\partial a_y} = \\delta_{xy}$.\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "498d13ec",
"metadata": {},
"source": [
"## Exercise 2 - Deriving the expression for OLS"
]
},
{
"cell_type": "markdown",
"id": "f3f771de",
"metadata": {},
"source": [
"The ordinary least squares method finds the parameters $\\boldsymbol{\\theta}$ which minimizes the squared error between our model $\\boldsymbol{X\\theta}$ and the true values $\\boldsymbol{y}$.\n",
"\n",
"To find the parameters $\\boldsymbol{\\theta}$ which minimizes this error, we take the derivative of the squared error expression with respect to $\\boldsymbol{\\theta}$, and set it equal to 0."
]
},
{
"cell_type": "markdown",
"id": "49690237",
"metadata": {},
"source": [
"**a)** Very briefly explain why the approach above finds the parameters $\\boldsymbol{\\theta}$ which minimizes this error."
]
},
{
"cell_type": "markdown",
"id": "b7cccc9d",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"\n",
"$$\n",
"\\vert\\vert\\boldsymbol{y} - \\boldsymbol{X\\theta}\\vert\\vert^2\n",
"$$\n",
"\n",
"which we can rewrite in matrix-vector form as\n",
"\n",
"$$\n",
"\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\theta}\\right)^T\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\theta}\\right)\n",
"$$.\n",
"If we take the derivative and set it to 0 we find extrema in the squared error term. Since the squared error is positivly definit this point will be the minimum. We minimized the squared error.\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "8fbecf74",
"metadata": {},
"source": [
"**b)** If $\\boldsymbol{X}$ is invertible, what is the expression for the optimal parameters $\\boldsymbol{\\theta}$? (**Hint:** Don't compute any derivatives, but solve $\\boldsymbol{X\\theta}=\\boldsymbol{y}$ for $\\boldsymbol{\\theta}$)"
]
},
{
"cell_type": "markdown",
"id": "0a4e9afe",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"\n",
"$$\n",
"\\theta = X^{-1} y\n",
"$$\n",
"since then $$X \\theta = X X^{-1} y = y$$\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "f37af8f0",
"metadata": {},
"source": [
"**c)** Show that\n",
"\n",
"$$\n",
"\\frac{\\partial \\left(\\boldsymbol{x}-\\boldsymbol{A}\\boldsymbol{s}\\right)^T\\left(\\boldsymbol{x}-\\boldsymbol{A}\\boldsymbol{s}\\right)}{\\partial \\boldsymbol{s}} = -2\\left(\\boldsymbol{x}-\\boldsymbol{A}\\boldsymbol{s}\\right)^T\\boldsymbol{A},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "a7e35dca",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"\\frac{\\partial}{\\partial s} (x-As)^T(x-As) &= \\frac{\\partial (x-As)^T}{\\partial s} (x-As) + (x-As)^T \\frac{\\partial (x-As)}{\\partial s} \\\\\n",
"&= - A^T (x-As) - (x-As)^T A \\\\\n",
"&= -2 (x-As)^T A\n",
"\\end{aligned}\n",
"$$\n",
"We can make the last step due to the fact, that $(x-As)^T(x-As) \\equiv \\gamma$ is a scalar and thus $\\gamma = \\gamma^T$. Therefore the term for the derivative of $\\gamma^T$ has to be equivalent and we can derive the last equivalence.\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "869fca4d",
"metadata": {},
"source": [
"**d)** Using the expression from **c)**, but substituting back in $\\boldsymbol{\\theta}$, $\\boldsymbol{y}$ and $\\boldsymbol{X}$, find the expression for the optimal parameters $\\boldsymbol{\\theta}$ in the case that $\\boldsymbol{X}$ is not invertible, but $\\boldsymbol{X^T X}$ is, which is most often the case."
]
},
{
"cell_type": "markdown",
"id": "aee5a457",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"\n",
"$$\n",
"\\boldsymbol{\\hat{\\theta}_{OLS}} = (X^T X)^{-1} X^T y\n",
"$$\n",
"\n",
"since then\n",
"$X\\theta = X (X^TX)^{-1} X^T y = (X X^{-1}) ((X^T)^{-1} X^T) y = y$\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "57ca3d74",
"metadata": {},
"source": [
"## Exercise 3 - Creating feature matrix and implementing OLS using the analytical expression"
]
},
{
"cell_type": "markdown",
"id": "5dc179f7",
"metadata": {},
"source": [
"With the expression for $\\boldsymbol{\\hat{\\theta}_{OLS}}$, you now have what you need to implement OLS regression with your input data and target data $\\boldsymbol{y}$. But before you can do that, you need to set up you input data as a feature matrix $\\boldsymbol{X}$.\n",
"\n",
"In a feature matrix, each row is a datapoint and each column is a feature of that data. If you want to predict someones spending based on their income and number of children, for instance, you would create a row for each person in your dataset, with the montly income and the number of children as columns.\n",
"\n",
"We typically also include an intercept in our models. The intercept is a value that is added to our prediction regardless of the value of the other features. The intercept tries to account for constant effects in our data that are not dependant on anything else. In our current example, the intercept could account for living expenses which are typical regardless of income or childcare expenses.\n",
"\n",
"We calculate the optimal intercept by including a feature with the constant value of 1 in our model, which is then multplied by some parameter $\\theta_0$ from the OLS method into the optimal intercept value (which will be $\\theta_0$). In practice, we include the intercept in our model by adding a column of ones to the start of our feature matrix."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "e5ff2a69",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "a3cf2792",
"metadata": {},
"outputs": [],
"source": [
"n = 20\n",
"income = np.array([116., 161., 167., 118., 172., 163., 179., 173., 162., 116., 101., 176., 178., 172., 143., 135., 160., 101., 149., 125.])\n",
"children = np.array([5, 3, 0, 4, 5, 3, 0, 4, 4, 3, 3, 5, 1, 0, 2, 3, 2, 1, 5, 4])\n",
"spending = np.array([152., 141., 102., 136., 161., 129., 99., 159., 160., 107., 98., 164., 121., 93., 112., 127., 117., 69., 156., 131.])\n"
]
},
{
"cell_type": "markdown",
"id": "5da61481",
"metadata": {},
"source": [
"**a)** Create a feature matrix $\\boldsymbol{X}$ for the features income and children, including an intercept column of ones at the start."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "5ad87a65",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1., 116., 5.],\n",
" [ 1., 161., 3.],\n",
" [ 1., 167., 0.],\n",
" [ 1., 118., 4.],\n",
" [ 1., 172., 5.],\n",
" [ 1., 163., 3.],\n",
" [ 1., 179., 0.],\n",
" [ 1., 173., 4.],\n",
" [ 1., 162., 4.],\n",
" [ 1., 116., 3.],\n",
" [ 1., 101., 3.],\n",
" [ 1., 176., 5.],\n",
" [ 1., 178., 1.],\n",
" [ 1., 172., 0.],\n",
" [ 1., 143., 2.],\n",
" [ 1., 135., 3.],\n",
" [ 1., 160., 2.],\n",
" [ 1., 101., 1.],\n",
" [ 1., 149., 5.],\n",
" [ 1., 125., 4.]])"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"X = np.stack((np.ones(n), income, children)).T\n",
"display(X)"
]
},
{
"cell_type": "markdown",
"id": "e0ddfac2",
"metadata": {},
"source": [
"**b)** Use the expression from **3d)** to find the optimal parameters $\\boldsymbol{\\hat{\\beta}_{OLS}}$ for predicting spending based on these features. Create a function for this operation, as you are going to need to use it a lot."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8f3f68aa",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 9.12808583, 0.5119025 , 14.60743095])"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def OLS_parameters(X, y):\n",
" return np.linalg.inv(X.T @ X) @ X.T @ y\n",
"\n",
"beta = OLS_parameters(X, spending)\n",
"display(beta)"
]
},
{
"cell_type": "markdown",
"id": "0cb6da80",
"metadata": {},
"source": [
"## Exercise 4 - Fitting a polynomial"
]
},
{
"cell_type": "markdown",
"id": "71015064",
"metadata": {},
"source": [
"In this course, we typically do linear regression using polynomials, though in real world applications it is also very common to make linear models based on measured features like you did in the previous exercise.\n",
"\n",
"When fitting a polynomial with linear regression, we make each polynomial degree($x, x^2, x^3, ..., x^p$) its own feature."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "d7476c84",
"metadata": {},
"outputs": [],
"source": [
"n = 100\n",
"x = np.linspace(-3, 3, n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + np.random.normal(n)"
]
},
{
"cell_type": "markdown",
"id": "8321451b",
"metadata": {},
"source": [
"**a)** Create a feature matrix $\\boldsymbol{X}$ for the features $x, x^2, x^3, x^4, x^5$, including an intercept column of ones at the start. Make this into a function, as you will do this a lot over the next weeks."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "91496e40",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1.00000000e+00, -3.00000000e+00, 9.00000000e+00,\n",
" -2.70000000e+01, 8.10000000e+01, -2.43000000e+02],\n",
" [ 1.00000000e+00, -2.93939394e+00, 8.64003673e+00,\n",
" -2.53964716e+01, 7.46502347e+01, -2.19426447e+02],\n",
" [ 1.00000000e+00, -2.87878788e+00, 8.28741965e+00,\n",
" -2.38577232e+01, 6.86813245e+01, -1.97718964e+02],\n",
" [ 1.00000000e+00, -2.81818182e+00, 7.94214876e+00,\n",
" -2.23824192e+01, 6.30777269e+01, -1.77764503e+02],\n",
" [ 1.00000000e+00, -2.75757576e+00, 7.60422406e+00,\n",
" -2.09692239e+01, 5.78242235e+01, -1.59454677e+02],\n",
" [ 1.00000000e+00, -2.69696970e+00, 7.27364555e+00,\n",
" -1.96168016e+01, 5.29059195e+01, -1.42685662e+02],\n",
" [ 1.00000000e+00, -2.63636364e+00, 6.95041322e+00,\n",
" -1.83238167e+01, 4.83082440e+01, -1.27358098e+02],\n",
" [ 1.00000000e+00, -2.57575758e+00, 6.63452709e+00,\n",
" -1.70889334e+01, 4.40169497e+01, -1.13376992e+02],\n",
" [ 1.00000000e+00, -2.51515152e+00, 6.32598714e+00,\n",
" -1.59108162e+01, 4.00181133e+01, -1.00651618e+02],\n",
" [ 1.00000000e+00, -2.45454545e+00, 6.02479339e+00,\n",
" -1.47881292e+01, 3.62981354e+01, -8.90954232e+01],\n",
" [ 1.00000000e+00, -2.39393939e+00, 5.73094582e+00,\n",
" -1.37195370e+01, 3.28437400e+01, -7.86259231e+01],\n",
" [ 1.00000000e+00, -2.33333333e+00, 5.44444444e+00,\n",
" -1.27037037e+01, 2.96419753e+01, -6.91646091e+01],\n",
" [ 1.00000000e+00, -2.27272727e+00, 5.16528926e+00,\n",
" -1.17392938e+01, 2.66802131e+01, -6.06368480e+01],\n",
" [ 1.00000000e+00, -2.21212121e+00, 4.89348026e+00,\n",
" -1.08249715e+01, 2.39461490e+01, -5.29717842e+01],\n",
" [ 1.00000000e+00, -2.15151515e+00, 4.62901745e+00,\n",
" -9.95940117e+00, 2.14278025e+01, -4.61022418e+01],\n",
" [ 1.00000000e+00, -2.09090909e+00, 4.37190083e+00,\n",
" -9.14124718e+00, 1.91135168e+01, -3.99646261e+01],\n",
" [ 1.00000000e+00, -2.03030303e+00, 4.12213039e+00,\n",
" -8.36917383e+00, 1.69919590e+01, -3.44988258e+01],\n",
" [ 1.00000000e+00, -1.96969697e+00, 3.87970615e+00,\n",
" -7.64184545e+00, 1.50521198e+01, -2.96481148e+01],\n",
" [ 1.00000000e+00, -1.90909091e+00, 3.64462810e+00,\n",
" -6.95792637e+00, 1.32833140e+01, -2.53590540e+01],\n",
" [ 1.00000000e+00, -1.84848485e+00, 3.41689624e+00,\n",
" -6.31608092e+00, 1.16751799e+01, -2.15813931e+01],\n",
" [ 1.00000000e+00, -1.78787879e+00, 3.19651056e+00,\n",
" -5.71497343e+00, 1.02176798e+01, -1.82679729e+01],\n",
" [ 1.00000000e+00, -1.72727273e+00, 2.98347107e+00,\n",
" -5.15326822e+00, 8.90109965e+00, -1.53746267e+01],\n",
" [ 1.00000000e+00, -1.66666667e+00, 2.77777778e+00,\n",
" -4.62962963e+00, 7.71604938e+00, -1.28600823e+01],\n",
" [ 1.00000000e+00, -1.60606061e+00, 2.57943067e+00,\n",
" -4.14272199e+00, 6.65346258e+00, -1.06858641e+01],\n",
" [ 1.00000000e+00, -1.54545455e+00, 2.38842975e+00,\n",
" -3.69120962e+00, 5.70459668e+00, -8.81619487e+00],\n",
" [ 1.00000000e+00, -1.48484848e+00, 2.20477502e+00,\n",
" -3.27375685e+00, 4.86103290e+00, -7.21789734e+00],\n",
" [ 1.00000000e+00, -1.42424242e+00, 2.02846648e+00,\n",
" -2.88902802e+00, 4.11467627e+00, -5.86029651e+00],\n",
" [ 1.00000000e+00, -1.36363636e+00, 1.85950413e+00,\n",
" -2.53568745e+00, 3.45775562e+00, -4.71512130e+00],\n",
" [ 1.00000000e+00, -1.30303030e+00, 1.69788797e+00,\n",
" -2.21239948e+00, 2.88282356e+00, -3.75640646e+00],\n",
" [ 1.00000000e+00, -1.24242424e+00, 1.54361800e+00,\n",
" -1.91782842e+00, 2.38275652e+00, -2.96039447e+00],\n",
" [ 1.00000000e+00, -1.18181818e+00, 1.39669421e+00,\n",
" -1.65063862e+00, 1.95075473e+00, -2.30543741e+00],\n",
" [ 1.00000000e+00, -1.12121212e+00, 1.25711662e+00,\n",
" -1.40949439e+00, 1.58034220e+00, -1.77189883e+00],\n",
" [ 1.00000000e+00, -1.06060606e+00, 1.12488522e+00,\n",
" -1.19306008e+00, 1.26536675e+00, -1.34205564e+00],\n",
" [ 1.00000000e+00, -1.00000000e+00, 1.00000000e+00,\n",
" -1.00000000e+00, 1.00000000e+00, -1.00000000e+00],\n",
" [ 1.00000000e+00, -9.39393939e-01, 8.82460973e-01,\n",
" -8.28978490e-01, 7.78737370e-01, -7.31541165e-01],\n",
" [ 1.00000000e+00, -8.78787879e-01, 7.72268136e-01,\n",
" -6.78659877e-01, 5.96398074e-01, -5.24107398e-01],\n",
" [ 1.00000000e+00, -8.18181818e-01, 6.69421488e-01,\n",
" -5.47708490e-01, 4.48125128e-01, -3.66647832e-01],\n",
" [ 1.00000000e+00, -7.57575758e-01, 5.73921028e-01,\n",
" -4.34788658e-01, 3.29385347e-01, -2.49534354e-01],\n",
" [ 1.00000000e+00, -6.96969697e-01, 4.85766758e-01,\n",
" -3.38564710e-01, 2.35969344e-01, -1.64463482e-01],\n",
" [ 1.00000000e+00, -6.36363636e-01, 4.04958678e-01,\n",
" -2.57700977e-01, 1.63991531e-01, -1.04358247e-01],\n",
" [ 1.00000000e+00, -5.75757576e-01, 3.31496786e-01,\n",
" -1.90861786e-01, 1.09890119e-01, -6.32700686e-02],\n",
" [ 1.00000000e+00, -5.15151515e-01, 2.65381084e-01,\n",
" -1.36711467e-01, 7.04271195e-02, -3.62806373e-02],\n",
" [ 1.00000000e+00, -4.54545455e-01, 2.06611570e-01,\n",
" -9.39143501e-02, 4.26883410e-02, -1.94037913e-02],\n",
" [ 1.00000000e+00, -3.93939394e-01, 1.55188246e-01,\n",
" -6.11347636e-02, 2.40833917e-02, -9.48739674e-03],\n",
" [ 1.00000000e+00, -3.33333333e-01, 1.11111111e-01,\n",
" -3.70370370e-02, 1.23456790e-02, -4.11522634e-03],\n",
" [ 1.00000000e+00, -2.72727273e-01, 7.43801653e-02,\n",
" -2.02854996e-02, 5.53240899e-03, -1.50883882e-03],\n",
" [ 1.00000000e+00, -2.12121212e-01, 4.49954086e-02,\n",
" -9.54448062e-03, 2.02458680e-03, -4.29457806e-04],\n",
" [ 1.00000000e+00, -1.51515152e-01, 2.29568411e-02,\n",
" -3.47830926e-03, 5.27016555e-04, -7.98509932e-05],\n",
" [ 1.00000000e+00, -9.09090909e-02, 8.26446281e-03,\n",
" -7.51314801e-04, 6.83013455e-05, -6.20921323e-06],\n",
" [ 1.00000000e+00, -3.03030303e-02, 9.18273646e-04,\n",
" -2.78264741e-05, 8.43226488e-07, -2.55523178e-08],\n",
" [ 1.00000000e+00, 3.03030303e-02, 9.18273646e-04,\n",
" 2.78264741e-05, 8.43226488e-07, 2.55523178e-08],\n",
" [ 1.00000000e+00, 9.09090909e-02, 8.26446281e-03,\n",
" 7.51314801e-04, 6.83013455e-05, 6.20921323e-06],\n",
" [ 1.00000000e+00, 1.51515152e-01, 2.29568411e-02,\n",
" 3.47830926e-03, 5.27016555e-04, 7.98509932e-05],\n",
" [ 1.00000000e+00, 2.12121212e-01, 4.49954086e-02,\n",
" 9.54448062e-03, 2.02458680e-03, 4.29457806e-04],\n",
" [ 1.00000000e+00, 2.72727273e-01, 7.43801653e-02,\n",
" 2.02854996e-02, 5.53240899e-03, 1.50883882e-03],\n",
" [ 1.00000000e+00, 3.33333333e-01, 1.11111111e-01,\n",
" 3.70370370e-02, 1.23456790e-02, 4.11522634e-03],\n",
" [ 1.00000000e+00, 3.93939394e-01, 1.55188246e-01,\n",
" 6.11347636e-02, 2.40833917e-02, 9.48739674e-03],\n",
" [ 1.00000000e+00, 4.54545455e-01, 2.06611570e-01,\n",
" 9.39143501e-02, 4.26883410e-02, 1.94037913e-02],\n",
" [ 1.00000000e+00, 5.15151515e-01, 2.65381084e-01,\n",
" 1.36711467e-01, 7.04271195e-02, 3.62806373e-02],\n",
" [ 1.00000000e+00, 5.75757576e-01, 3.31496786e-01,\n",
" 1.90861786e-01, 1.09890119e-01, 6.32700686e-02],\n",
" [ 1.00000000e+00, 6.36363636e-01, 4.04958678e-01,\n",
" 2.57700977e-01, 1.63991531e-01, 1.04358247e-01],\n",
" [ 1.00000000e+00, 6.96969697e-01, 4.85766758e-01,\n",
" 3.38564710e-01, 2.35969344e-01, 1.64463482e-01],\n",
" [ 1.00000000e+00, 7.57575758e-01, 5.73921028e-01,\n",
" 4.34788658e-01, 3.29385347e-01, 2.49534354e-01],\n",
" [ 1.00000000e+00, 8.18181818e-01, 6.69421488e-01,\n",
" 5.47708490e-01, 4.48125128e-01, 3.66647832e-01],\n",
" [ 1.00000000e+00, 8.78787879e-01, 7.72268136e-01,\n",
" 6.78659877e-01, 5.96398074e-01, 5.24107398e-01],\n",
" [ 1.00000000e+00, 9.39393939e-01, 8.82460973e-01,\n",
" 8.28978490e-01, 7.78737370e-01, 7.31541165e-01],\n",
" [ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,\n",
" 1.00000000e+00, 1.00000000e+00, 1.00000000e+00],\n",
" [ 1.00000000e+00, 1.06060606e+00, 1.12488522e+00,\n",
" 1.19306008e+00, 1.26536675e+00, 1.34205564e+00],\n",
" [ 1.00000000e+00, 1.12121212e+00, 1.25711662e+00,\n",
" 1.40949439e+00, 1.58034220e+00, 1.77189883e+00],\n",
" [ 1.00000000e+00, 1.18181818e+00, 1.39669421e+00,\n",
" 1.65063862e+00, 1.95075473e+00, 2.30543741e+00],\n",
" [ 1.00000000e+00, 1.24242424e+00, 1.54361800e+00,\n",
" 1.91782842e+00, 2.38275652e+00, 2.96039447e+00],\n",
" [ 1.00000000e+00, 1.30303030e+00, 1.69788797e+00,\n",
" 2.21239948e+00, 2.88282356e+00, 3.75640646e+00],\n",
" [ 1.00000000e+00, 1.36363636e+00, 1.85950413e+00,\n",
" 2.53568745e+00, 3.45775562e+00, 4.71512130e+00],\n",
" [ 1.00000000e+00, 1.42424242e+00, 2.02846648e+00,\n",
" 2.88902802e+00, 4.11467627e+00, 5.86029651e+00],\n",
" [ 1.00000000e+00, 1.48484848e+00, 2.20477502e+00,\n",
" 3.27375685e+00, 4.86103290e+00, 7.21789734e+00],\n",
" [ 1.00000000e+00, 1.54545455e+00, 2.38842975e+00,\n",
" 3.69120962e+00, 5.70459668e+00, 8.81619487e+00],\n",
" [ 1.00000000e+00, 1.60606061e+00, 2.57943067e+00,\n",
" 4.14272199e+00, 6.65346258e+00, 1.06858641e+01],\n",
" [ 1.00000000e+00, 1.66666667e+00, 2.77777778e+00,\n",
" 4.62962963e+00, 7.71604938e+00, 1.28600823e+01],\n",
" [ 1.00000000e+00, 1.72727273e+00, 2.98347107e+00,\n",
" 5.15326822e+00, 8.90109965e+00, 1.53746267e+01],\n",
" [ 1.00000000e+00, 1.78787879e+00, 3.19651056e+00,\n",
" 5.71497343e+00, 1.02176798e+01, 1.82679729e+01],\n",
" [ 1.00000000e+00, 1.84848485e+00, 3.41689624e+00,\n",
" 6.31608092e+00, 1.16751799e+01, 2.15813931e+01],\n",
" [ 1.00000000e+00, 1.90909091e+00, 3.64462810e+00,\n",
" 6.95792637e+00, 1.32833140e+01, 2.53590540e+01],\n",
" [ 1.00000000e+00, 1.96969697e+00, 3.87970615e+00,\n",
" 7.64184545e+00, 1.50521198e+01, 2.96481148e+01],\n",
" [ 1.00000000e+00, 2.03030303e+00, 4.12213039e+00,\n",
" 8.36917383e+00, 1.69919590e+01, 3.44988258e+01],\n",
" [ 1.00000000e+00, 2.09090909e+00, 4.37190083e+00,\n",
" 9.14124718e+00, 1.91135168e+01, 3.99646261e+01],\n",
" [ 1.00000000e+00, 2.15151515e+00, 4.62901745e+00,\n",
" 9.95940117e+00, 2.14278025e+01, 4.61022418e+01],\n",
" [ 1.00000000e+00, 2.21212121e+00, 4.89348026e+00,\n",
" 1.08249715e+01, 2.39461490e+01, 5.29717842e+01],\n",
" [ 1.00000000e+00, 2.27272727e+00, 5.16528926e+00,\n",
" 1.17392938e+01, 2.66802131e+01, 6.06368480e+01],\n",
" [ 1.00000000e+00, 2.33333333e+00, 5.44444444e+00,\n",
" 1.27037037e+01, 2.96419753e+01, 6.91646091e+01],\n",
" [ 1.00000000e+00, 2.39393939e+00, 5.73094582e+00,\n",
" 1.37195370e+01, 3.28437400e+01, 7.86259231e+01],\n",
" [ 1.00000000e+00, 2.45454545e+00, 6.02479339e+00,\n",
" 1.47881292e+01, 3.62981354e+01, 8.90954232e+01],\n",
" [ 1.00000000e+00, 2.51515152e+00, 6.32598714e+00,\n",
" 1.59108162e+01, 4.00181133e+01, 1.00651618e+02],\n",
" [ 1.00000000e+00, 2.57575758e+00, 6.63452709e+00,\n",
" 1.70889334e+01, 4.40169497e+01, 1.13376992e+02],\n",
" [ 1.00000000e+00, 2.63636364e+00, 6.95041322e+00,\n",
" 1.83238167e+01, 4.83082440e+01, 1.27358098e+02],\n",
" [ 1.00000000e+00, 2.69696970e+00, 7.27364555e+00,\n",
" 1.96168016e+01, 5.29059195e+01, 1.42685662e+02],\n",
" [ 1.00000000e+00, 2.75757576e+00, 7.60422406e+00,\n",
" 2.09692239e+01, 5.78242235e+01, 1.59454677e+02],\n",
" [ 1.00000000e+00, 2.81818182e+00, 7.94214876e+00,\n",
" 2.23824192e+01, 6.30777269e+01, 1.77764503e+02],\n",
" [ 1.00000000e+00, 2.87878788e+00, 8.28741965e+00,\n",
" 2.38577232e+01, 6.86813245e+01, 1.97718964e+02],\n",
" [ 1.00000000e+00, 2.93939394e+00, 8.64003673e+00,\n",
" 2.53964716e+01, 7.46502347e+01, 2.19426447e+02],\n",
" [ 1.00000000e+00, 3.00000000e+00, 9.00000000e+00,\n",
" 2.70000000e+01, 8.10000000e+01, 2.43000000e+02]])"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def polynomial_features(x, p):\n",
" n = len(x)\n",
" X = np.power(x[:, np.newaxis], np.arange(p + 1))\n",
" return X\n",
"\n",
"X = polynomial_features(x, 5)\n",
"display(X)"
]
},
{
"cell_type": "markdown",
"id": "b84b1e31",
"metadata": {},
"source": [
"**b)** Use the expression from **3d)** to find the optimal parameters $\\boldsymbol{\\hat{\\beta}_{OLS}}$ for predicting $\\boldsymbol{y}$ based on these features. If you have done everything right so far, this code will not need changing."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "034f502c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.92452576, 0.27464654, -0.02326439, 0.05342623, -0.0034652 ,\n",
" -0.0087781 ])"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"beta = OLS_parameters(X, y)\n",
"display(beta)"
]
},
{
"cell_type": "markdown",
"id": "d703f788",
"metadata": {},
"source": [
"**c)** Like in exercise 4 last week, split your feature matrix and target data into a training split and test split."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "29171358",
"metadata": {},
"outputs": [],
"source": [
"from sklearn.model_selection import train_test_split\n",
"\n",
"x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)"
]
},
{
"cell_type": "markdown",
"id": "a0e3509f",
"metadata": {},
"source": [
"**d)** Train your model on the training data(find the parameters which best fit) and compute the MSE on both the training and test data."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "1e346f4c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Training MSE: 0.0133\n",
"Testing MSE: 0.0162\n"
]
}
],
"source": [
"beta = OLS_parameters(polynomial_features(x_train, 5), y_train)\n",
"\n",
"from sklearn.metrics import mean_squared_error\n",
"\n",
"def evaluate_model(beta, X, y):\n",
" y_pred = X @ beta\n",
" mse = mean_squared_error(y, y_pred)\n",
" return mse\n",
"\n",
"mse_train = evaluate_model(beta, polynomial_features(x_train, 5), y_train)\n",
"mse_test = evaluate_model(beta, polynomial_features(x_test, 5), y_test)\n",
"print(f\"Training MSE: {mse_train:.4f}\")\n",
"print(f\"Testing MSE: {mse_test:.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "7e431889",
"metadata": {},
"source": [
"**e)** Do the same for each polynomial degree from 2 to 10, and plot the MSE on both the training and test data as a function of polynomial degree. The aim is to reproduce Figure 2.11 of [Hastie et al](https://github.com/CompPhysics/MLErasmus/blob/master/doc/Textbooks/elementsstat.pdf). Feel free to read the discussions leading to figure 2.11 of Hastie et al. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "ceb57457",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.microsoft.datawrangler.viewer.v0+json": {
"columns": [
{
"name": "index",
"rawType": "int64",
"type": "integer"
},
{
"name": "degree",
"rawType": "int64",
"type": "integer"
},
{
"name": "mse_train",
"rawType": "float64",
"type": "float"
},
{
"name": "mse_test",
"rawType": "float64",
"type": "float"
}
],
"ref": "66cfaf7d-1e2b-4459-a642-f9fc0decfff7",
"rows": [
[
"0",
"2",
"0.05164086416666284",
"0.016256195807635355"
],
[
"1",
"3",
"0.02206575873360362",
"0.0203940940253871"
],
[
"2",
"4",
"0.020904375395772917",
"0.023035341621416943"
],
[
"3",
"5",
"0.013347385151109586",
"0.01615480202876137"
],
[
"4",
"6",
"0.00955995213624267",
"0.008151629748962193"
],
[
"5",
"7",
"0.00574771913202137",
"0.006163149518752208"
],
[
"6",
"8",
"0.0010834910069786221",
"0.0008846272939455341"
],
[
"7",
"9",
"0.000958469368562645",
"0.001134083627637028"
],
[
"8",
"10",
"7.94173975842872e-05",
"9.913320977908047e-05"
]
],
"shape": {
"columns": 3,
"rows": 9
}
},
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>degree</th>\n",
" <th>mse_train</th>\n",
" <th>mse_test</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>2</td>\n",
" <td>0.051641</td>\n",
" <td>0.016256</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>3</td>\n",
" <td>0.022066</td>\n",
" <td>0.020394</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>4</td>\n",
" <td>0.020904</td>\n",
" <td>0.023035</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>5</td>\n",
" <td>0.013347</td>\n",
" <td>0.016155</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>6</td>\n",
" <td>0.009560</td>\n",
" <td>0.008152</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>7</td>\n",
" <td>0.005748</td>\n",
" <td>0.006163</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>8</td>\n",
" <td>0.001083</td>\n",
" <td>0.000885</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>9</td>\n",
" <td>0.000958</td>\n",
" <td>0.001134</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>10</td>\n",
" <td>0.000079</td>\n",
" <td>0.000099</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" degree mse_train mse_test\n",
"0 2 0.051641 0.016256\n",
"1 3 0.022066 0.020394\n",
"2 4 0.020904 0.023035\n",
"3 5 0.013347 0.016155\n",
"4 6 0.009560 0.008152\n",
"5 7 0.005748 0.006163\n",
"6 8 0.001083 0.000885\n",
"7 9 0.000958 0.001134\n",
"8 10 0.000079 0.000099"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import pandas as pd\n",
"\n",
"results = []\n",
"\n",
"for degree in range(2, 11):\n",
" beta = OLS_parameters(polynomial_features(x_train, degree), y_train)\n",
" mse_train = evaluate_model(beta, polynomial_features(x_train, degree), y_train)\n",
" mse_test = evaluate_model(beta, polynomial_features(x_test, degree), y_test)\n",
" results.append({\"degree\": degree, \"mse_train\": mse_train, \"mse_test\": mse_test})\n",
"\n",
"df_results = pd.DataFrame(results)\n",
"display(df_results)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "c3d81800",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fa9522bfe00>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAgwVJREFUeJzt3XdcleX/x/HXOey9pyKIorhyD9wDw5mj0sxyoPbLTC3TUr/lyNKytKVlmbssNc3KzHLnIBdqLsCNA0QcICDz3L8/Dpw8gcaBA4fxeT4e5yHnPte5z/sAej5e9zVUiqIoCCGEEEJUImpTBxBCCCGEKG1SAAkhhBCi0pECSAghhBCVjhRAQgghhKh0pAASQgghRKUjBZAQQgghKh0pgIQQQghR6ZibOkBZpNFouH79Og4ODqhUKlPHEUIIIUQhKIrCvXv38PX1Ra1+dB+PFEAFuH79On5+fqaOIYQQQogiuHLlClWrVn1kGymACuDg4AC530BHR0dTxxFCCCFEISQnJ+Pn56f7HH8UKYAKkHfZy9HRUQogIYQQopwpzPAVGQQthBBCiEpHCiAhhBBCVDpSAAkhhBCi0pExQEIIIcqUnJwcsrKyTB1DlEEWFhaYmZkZ5VxSAAkhhCgTFEUhPj6eu3fvmjqKKMOcnZ3x9vYu9jp9UgAJIYQoE/KKH09PT2xtbWUhWqFHURTS0tJISEgAwMfHp1jnkwJICCGEyeXk5OiKHzc3N1PHEWWUjY0NAAkJCXh6ehbrcpgMghZCCGFyeWN+bG1tTR1FlHF5vyPFHScmBZAQQogyQy57if9irN8RKYCEEEIIUelIASSEEEKISkcKICGEEKKMCQgI4OOPPy50+127dqFSqWQJAQNIAVTKLiamcv3ufVPHEEIIYQQqleqRtxkzZhTpvIcOHeKFF14odPvWrVsTFxeHk5NTkV6vsPIKLRcXF9LT0/UeO3TokO59P2jx4sU0bNgQe3t7nJ2dady4MXPmzNE9PmPGjAK/d8HBwSX6XmQafCmatek0S/Ze5P86BDKlex1TxxFCCFFMcXFxuq/XrFnDtGnTiI6O1h2zt7fXfa0oCjk5OZib//dHr4eHh0E5LC0t8fb2Nug5xeHg4MCPP/7IoEGDdMeWLFlCtWrViI2N1R1bunQpr7zyCp9++ikdOnQgIyODv//+m5MnT+qdr169emzbtk3vWGG+T8UhPUCl6LGq2sp8Z1SCqaMIIUSZpygKaZnZJrkpilKojN7e3rqbk5MTKpVKdz8qKgoHBwd+++03mjZtipWVFXv37uX8+fP06dMHLy8v7O3tad68eb4P/39fAlOpVHz99df069cPW1tbgoKC+Pnnn3WP//sS2PLly3F2dub333+nTp062Nvb061bN72CLTs7m3HjxuHs7IybmxtvvPEGQ4cOpW/fvv/5vocOHcrSpUt19+/fv8/333/P0KFD9dr9/PPPDBgwgBEjRlCzZk3q1avHoEGDePfdd/XamZub630vvb29cXd3L9TPoKikB6gUdajlgVoFMTdSuHonjaoust6FEEI8zP2sHOpO+90kr3367TBsLY3zETl58mQ+/PBDAgMDcXFx4cqVK/To0YN3330XKysrVq5cSe/evYmOjqZatWoPPc/MmTOZO3cuH3zwAZ999hmDBw/m8uXLuLq6Ftg+LS2NDz/8kFWrVqFWq3nuueeYOHEi3377LQDvv/8+3377LcuWLaNOnTp88sknbNy4kU6dOv3ne3r++ef54IMPiI2NpVq1aqxfv56AgACaNGmi187b25vdu3dz+fJl/P39Df7elSTpASpFzraWNPV3AekFEkKISuPtt9+ma9eu1KhRA1dXVxo2bMj//d//Ub9+fYKCgpg1axY1atTQ69EpyLBhwxg0aBA1a9Zk9uzZpKSkcPDgwYe2z8rKYtGiRTRr1owmTZrw8ssvs337dt3jn332GVOmTKFfv34EBwezYMECnJ2dC/WePD096d69O8uXL4fcS13h4eH52k2fPh1nZ2cCAgKoXbs2w4YNY+3atWg0Gr12J06cwN7eXu/24osvFipLUUkPUCnrFOzJoUt32BGVwPMhAaaOI4QQZZaNhRmn3w4z2WsbS7NmzfTup6SkMGPGDH799Vfi4uLIzs7m/v37emNnCvLYY4/pvrazs8PR0VG3L1ZBbG1tqVGjhu6+j4+Prn1SUhI3btygRYsWusfNzMxo2rRpvuLkYcLDwxk/fjzPPfccERERrFu3jj179ui18fHxISIigpMnT/Lnn3+yf/9+hg4dytdff82WLVtQq7X9MLVr185XADo6OhYqR1FJAVTKOgd7MndLNPvP3yI9KwdrI/4lE0KIikSlUhntMpQp2dnZ6d2fOHEiW7du5cMPP6RmzZrY2Njw1FNPkZmZ+cjzWFhY6N1XqVSPLFYKal/YsU2F0b17d1544QVGjBhB7969H7mHW/369alfvz4vvfQSL774Iu3atWP37t26y22WlpbUrFnTaNkKQy6BlbLaXg74OlmTka0h4vwtU8cRQghRyvbt28ewYcPo168fDRo0wNvbm0uXLpVqBicnJ7y8vDh06JDuWE5ODpGRkYU+h7m5OUOGDGHXrl0FXv56mLp16wKQmppqYGrjKv+ldTmjUqnoGOzJ6gOx7IhKoFOwp6kjCSGEKEVBQUFs2LCB3r17o1KpeOuttwp92cmYxo4dy5w5c6hZsybBwcF89tln3Llzx6C9tmbNmsWkSZMe2vszevRofH196dy5M1WrViUuLo533nkHDw8PQkJCdO2ys7OJj4/Xe65KpcLLy6sY7/DRykQP0MKFCwkICMDa2pqWLVs+clAXwLp16wgODsba2poGDRqwefNmvceHDRuWb0Glbt26lfC7KLzOtbVFz46oBKN2RwohhCj75s+fj4uLC61bt6Z3796EhYXlmz1VGt544w0GDRrEkCFDCAkJwd7enrCwMKytrQt9DktLS9zd3R9aNIWGhvLXX3/x9NNPU6tWLZ588kmsra3Zvn27XtF06tQpfHx89G4lPWtMpZj4E3jNmjUMGTKERYsW0bJlSz7++GPWrVtHdHQ0np75e0f2799P+/btmTNnDr169WL16tW8//77REZGUr9+fcgtgG7cuMGyZct0z7OyssLFxaVQmZKTk3FyciIpKalEBmGlZWbT6O2tZGZr+OPV9tTycjD6awghRHmSnp7OxYsXqV69ukEfwMJ4NBoNderUYcCAAcyaNcvUcR7qUb8rhnx+m7wHaP78+YwaNYrhw4dTt25dFi1ahK2trd4CSw/65JNP6NatG5MmTaJOnTrMmjWLJk2asGDBAr12VlZWegsqFbb4KQ22luaEBGor3x0yHV4IIYQJXL58mcWLFxMTE8OJEycYPXo0Fy9e5NlnnzV1tFJh0gIoMzOTI0eOEBoa+k8gtZrQ0FAiIiIKfE5ERIRee4CwsLB87Xft2oWnpye1a9dm9OjR3Lr18AHHGRkZJCcn691KWufgfy6DCSGEEKVNrVazfPlymjdvTps2bThx4gTbtm2jTp3KsVWTSQdBJyYmkpOTk2+Qk5eXF1FRUQU+Jz4+vsD2Dw6e6tatG/3796d69eqcP3+eqVOn0r17dyIiIjAzyz/tfM6cOcycOdNo76swOgd7Mv3nUxy5fIektCycbC0K8SwhhBDCOPz8/Ni3b5+pY5iMyS+BlYRnnnmGJ554ggYNGtC3b182bdrEoUOH2LVrV4Htp0yZQlJSku525cqVEs/o52pLTU97cjQKe87dLPHXE0IIIcQ/TFoAubu7Y2Zmxo0bN/SO37hx46G72np7exvUHiAwMBB3d3fOnTtX4ONWVlY4Ojrq3UqDXAYTQgghTMOkBZClpSVNmzbV25tEo9Gwfft2vfUBHhQSEqLXHmDr1q0PbQ9w9epVbt26hY+PjxHTF1/H2h4A7I6+iUYj0+GFEEKI0mLyS2ATJkxg8eLFrFixgjNnzjB69GhSU1MZPnw4AEOGDGHKlCm69uPHj2fLli3MmzePqKgoZsyYweHDh3n55Zchd4+VSZMm8ddff3Hp0iW2b99Onz59qFmzJmFhptlT5mGaB7jiYGXOrdRMjl+9a+o4QgghRKVh8pWgBw4cyM2bN5k2bRrx8fE0atSILVu26AY6x8bG6jZLA2jdujWrV6/mzTffZOrUqQQFBbFx40bdGkBmZmb8/fffrFixgrt37+Lr68vjjz/OrFmzsLKyMtn7LIiFmZp2tdzZfCKenVEJNK5WdqbqCyGEEBWZyRdCLItKeiHEB607fIVJP/xN/SqObBrbrkRfSwghyipZCLFwZsyYwcaNGzl27Jipo5hMhVkIsbLrmLstxslrySQkp5s6jhBCCAP8e9ulf99mzJhRrHNv3LhR79jEiRPzjYMtCTNmzHjoNlIffPCBdl/Ljh11x9LS0pgyZQo1atTA2toaDw8POnTowE8//aRr07FjxwK/Ry+++GKJv5+CmPwSWGXn4WBFw6pOHL+axK7omwxo7mfqSEIIIQopLi5O9/WaNWuYNm0a0dHRumP29vZGfT17e3ujn/NhfHx82LlzJ1evXqVq1aq640uXLqVatWp6bV988UUOHDjAZ599Rt26dbl16xb79+/PtwjxqFGjePvtt/WO2dralvA7KZj0AJUBnWQ6vBBClEsPbrnk5OSESqXSO/b9999Tp04drK2tCQ4O5vPPP9c9NzMzk5dffhkfHx+sra3x9/dnzpw5AAQEBADQr18/VCqV7v6MGTNo1KiR7hzDhg2jb9++fPjhh/j4+ODm5saYMWPIysrStYmLi6Nnz57Y2NhQvXp1Vq9eTUBAAB9//PEj35unpyePP/44K1as0B3bv38/iYmJ9OzZU6/tzz//zNSpU+nRowcBAQE0bdqUsWPHEh4ertfO1tZW7/vj7e1dakvP/Jv0AJUBnYM9+XjbWfaeSyQzW4OludSlQgiBokBWmmle28IWHrLDeWF9++23TJs2jQULFtC4cWOOHj3KqFGjsLOzY+jQoXz66af8/PPPrF27lmrVqnHlyhXdQryHDh3C09OTZcuW0a1btwJ3Mcizc+dOXW/NuXPnGDhwII0aNWLUqFGQO5s6MTGRXbt2YWFhwYQJE0hIKNx/uMPDw3n99df53//+B7m9P4MHD87Xztvbm82bN9O/f38cHMrHBt9SAJUB9X2dcLe3IjElg0OXbtOmprupIwkhhOllpcFsX9O89tTrYGlXrFNMnz6defPm0b9/fwCqV6/O6dOn+fLLLxk6dCixsbEEBQXRtm1bVCoV/v7+uud6eGjXiXN2dn7kQr8ALi4uLFiwADMzM4KDg+nZsyfbt29n1KhRREVFsW3bNg4dOkSzZs0A+PrrrwkKCirUe+jVqxcvvvgif/75J02bNmXt2rXs3bs334blX331FYMHD8bNzY2GDRvStm1bnnrqKdq0aaPX7vPPP+frr7/WO/bll18WWFSVNOlqKAPUapVuUUS5DCaEEOVfamoq58+fZ8SIEbpxO/b29rzzzjucP38eci9fHTt2jNq1azNu3Dj++OOPIr1WvXr19HqIfHx8dD080dHRmJub06RJE93jNWvWxMWlcMuuWFhY8Nxzz7Fs2TLWrVtHrVq1eOyxx/K1a9++PRcuXGD79u089dRTnDp1inbt2jFr1iy9doMHD+bYsWN6tyeeeKJI77u4pAeojOgc7MkPR66yMyqBt3rVNXUcIYQwPQtbbU+MqV67GFJSUgBYvHgxLVu21Hssr1hp0qQJFy9e5LfffmPbtm0MGDCA0NBQfvjhB8OiWuhvpq1SqdBoNMXK/6Dw8HBatmzJyZMn843p+XeOdu3a0a5dO9544w3eeecd3n77bd544w0sLS0BcHJyombNmkbLVhxSAJURbYPcMVeruJCYyqXEVALci9f1KoQQ5Z5KVezLUKbi5eWFr68vFy5ceOTlHUdHRwYOHMjAgQN56qmn6NatG7dv38bV1RULCwtycnKKlaN27dpkZ2dz9OhRmjZtCsC5c+e4c+dOoc9Rr1496tWrx99//82zzz5b6OfVrVuX7Oxs0tPTdQVQWSIFUBnhaG1B8wBXIi7cYkdUAuFtq5s6khBCiGKYOXMm48aNw8nJiW7dupGRkcHhw4e5c+cOEyZMYP78+fj4+NC4cWPUajXr1q3D29sbZ2dnyJ0Jtn37dtq0aYOVlVWhL1s9KDg4mNDQUF544QW++OILLCwseO2117CxsUFlwCDvHTt2kJWVpcv2bx07dmTQoEE0a9YMNzc3Tp8+zdSpU+nUqZPeLK+0tDTi4+P1nlvU91ZcMgaoDMnbHX5ntIwDEkKI8m7kyJF8/fXXLFu2jAYNGtChQweWL19O9era/+A6ODgwd+5cmjVrRvPmzbl06RKbN2/Wbf80b948tm7dip+fH40bNy5yjpUrV+Ll5UX79u3p168fo0aNwsHBwaAVt+3s7B5a/ACEhYWxYsUKHn/8cerUqcPYsWMJCwtj7dq1eu0WL16Mj4+P3m3QoEFFfm/FIVthFKA0t8J40LmEFELn78bSTM3RaV2xs5IOOiFE5SBbYZSeq1ev4ufnx7Zt2+jSpYup4xhMtsKogGp42OHnakNmjoZ95xJNHUcIIUQFsGPHDn7++WcuXrzI/v37eeaZZwgICKB9+/amjmZSUgCVISqVis615TKYEEII48nKymLq1KnUq1ePfv364eHhoVsUsTKTayxlTKdgT1ZEXGZn1E0URTFokJoQQgjxb2FhYYSFhZk6RpkjPUBlTKtAN2wszIhPTud0XLKp4wghhBAVkhRAZYy1hRltaroBsFNWhRZCVDIyL0f8F2P9jkgBVAZ10k2Hv2nqKEIIUSryxqOkpZlo81NRbuT9jhR3DJOMASqDOuUOhD4ae4c7qZm42JW9FTSFEMKYzMzMcHZ21u1hZWtrK2MghR5FUUhLSyMhIQFnZ2e9/c+KQgqgMsjX2YZgbwei4u+xO+YmfRtXMXUkIYQocXm7nucVQUIUxNnZWfe7UhxSAJVRnYI9iYq/x46oBCmAhBCVgkqlwsfHB09PT7KyskwdR5RBFhYWxe75ySMFUBnVOdiTL3adZ3fMTbJzNJibyXAtIUTlYGZmZrQPOSEeRj5Vy6jGfs442ViQdD+Lo1fumjqOEEIIUaFIAVRGmZup6VDLA4AdMh1eCCGEMCqDC6AVK1bw66+/6u6//vrrODs707p1ay5fvmzsfJWabnd4KYCEEEIIozK4AJo9ezY2NjYAREREsHDhQubOnYu7uzuvvvpqSWSstDrU8kCtgqj4e1y/e9/UcYQQQogKw+AC6MqVK9SsWROAjRs38uSTT/LCCy8wZ84c9uzZUxIZKy0XO0saV3MB2RxVCCGEMCqDCyB7e3tu3boFwB9//EHXrl0BsLa25v596aUwtk61teOA5DKYEEIIYTwGF0Bdu3Zl5MiRjBw5kpiYGHr06AHAqVOnCAgIKImMlVrethj7zt0iPSvH1HGEEEKICsHgAmjhwoWEhIRw8+ZN1q9fj5ubduPOI0eOMGjQoJLIWKnV9XHE29Ga+1k5/HXhlqnjCCGEEBWCSpGtd/NJTk7GycmJpKQkHB0dTR2HKRv+5ruDVxga4s/MPvVNHUcIIYQokwz5/C7SOkB79uzhueeeo3Xr1ly7dg2AVatWsXfv3qIlFo+UtznqjugEpF4VQgghis/gAmj9+vWEhYVhY2NDZGQkGRkZACQlJTF79uySyFjptanpjqWZmiu373P+Zqqp4wghhBDlnsEF0DvvvMOiRYtYvHgxFhYWuuNt2rQhMjLS2PkEYGdlTstAV5DZYEIIIYRRGFwARUdH0759+3zHnZycuHtX9qwqKXmrQsu2GEIIIUTxGVwAeXt7c+7cuXzH9+7dS2BgoLFyiX/JGwd06NJtktOzTB1HCCGEKNcMLoBGjRrF+PHjOXDgACqViuvXr/Ptt98yceJERo8eXTIpBQHudgS625GtUdh7NtHUcYQQQohyzdzQJ0yePBmNRkOXLl1IS0ujffv2WFlZMXHiRMaOHVsyKQXkLop4Ye9FdkQl0KOBj6njCCGEEOVWkdcByszM5Ny5c6SkpFC3bl3s7e2Nn85Eyto6QHn2nUtk8NcHcLe35ODUUNRqlakjCSGEEGWGIZ/fBvcA5bG0tKRu3bpFfbooguYBrthbmZOYksmJa0k09HM2dSQhhBCiXDK4AOrUqRMq1cN7Hnbs2FHcTOIhLM3VtK3pzpZT8eyMTpACSAghhCgigwdBN2rUiIYNG+pudevWJTMzk8jISBo0aFAyKYVO3nR4WQ9ICCGEKDqDe4A++uijAo/PmDGDlJQUY2QSj9CxtgcAx68mcfNeBh4OVqaOJIQQQpQ7RdoLrCDPPfccS5cuNdbpxEN4OlpTv4p2YNeuaOkFEkIIIYrCaAVQREQE1tbWxjqdeITOuYsi7pQCSAghhCgSgy+B9e/fX+++oijExcVx+PBh3nrrLWNmEw/RKdiTT3ecY09MIlk5GizMjFbHCiGEEJWCwQWQk5OT3n21Wk3t2rV5++23efzxx42ZTTxEw6rOuNlZcis1k0OXbtO6hrupIwkhhBDlisEF0LJly0omiSg0tVpFh9oebIi8xq7om1IACSGEEAaSayfllOwOL4QQQhRdoXqAXFxcHrn44YNu375d3EyiENoFeWCmVnEuIYUrt9Pwc7U1dSQhhBCi3ChUAfTxxx+XfBJhECcbC5r6u3Dw4m12RCUwtHWAqSMJIYQQ5UahCqChQ4eWfBJhsM7BnlIACSGEEEVQrDFA6enpJCcn691E6ckbBxRx4RZpmdmmjiOEEEKUGwYXQKmpqbz88st4enpiZ2eHi4uL3k2UniBPe6o425CZrWH/uVumjiOEEEKUGwYXQK+//jo7duzgiy++wMrKiq+//pqZM2fi6+vLypUrSyalKJBKpfpnNpisCi2EEEIUmsEF0C+//MLnn3/Ok08+ibm5Oe3atePNN99k9uzZfPvtt0UKsXDhQgICArC2tqZly5YcPHjwke3XrVtHcHAw1tbWNGjQgM2bNz+07YsvvohKpaqwA7nzCqBdUQkoimLqOEIIIUS5YHABdPv2bQIDAwFwdHTUTXtv27Ytf/75p8EB1qxZw4QJE5g+fTqRkZE0bNiQsLAwEhIK7tHYv38/gwYNYsSIERw9epS+ffvSt29fTp48ma/tjz/+yF9//YWvr6/BucqLkBpuWFuouZ6UTvSNe6aOI4QQQpQLBhdAgYGBXLx4EYDg4GDWrl0LuT1Dzs7OBgeYP38+o0aNYvjw4dStW5dFixZha2v70J3lP/nkE7p168akSZOoU6cOs2bNokmTJixYsECv3bVr1xg7dizffvstFhYWBucqL6wtzHQrQcuiiEIIIUThGFwADR8+nOPHjwMwefJkFi5ciLW1Na+++iqTJk0y6FyZmZkcOXKE0NDQfwKp1YSGhhIREVHgcyIiIvTaA4SFhem112g0PP/880yaNIl69er9Z46MjIxyPZutU20PAHZKASSEEEIUSqH3Aps4cSIjR47k1Vdf1R0LDQ0lKiqKI0eOULNmTR577DGDXjwxMZGcnBy8vLz0jnt5eREVFVXgc+Lj4wtsHx8fr7v//vvvY25uzrhx4wqVY86cOcycOdOg7GVJp2BP+OkURy7f4W5aJs62lqaOJIQQQpRphe4B+umnn6hXrx6tW7dm6dKlpKamAuDv70///v0NLn5KypEjR/jkk09Yvnx5obfvmDJlCklJSbrblStXSjynMVV1saWWlz0aBXbH3DR1HCGEEKLMK3QBdPbsWXbu3EmtWrUYP3483t7ehIeHs3///iK/uLu7O2ZmZty4cUPv+I0bN/D29i7wOd7e3o9sv2fPHhISEqhWrRrm5uaYm5tz+fJlXnvtNQICCl4t2crKCkdHR71bedMpdzaYXAYTQggh/ptBY4Dat2/P8uXLiY+P55NPPuHs2bO0bduWOnXq8OGHH+YrTP6LpaUlTZs2Zfv27bpjGo2G7du3ExISUuBzQkJC9NoDbN26Vdf++eef5++//+bYsWO6m6+vL5MmTeL33383KF950rm2tgDaHXOTHI1MhxdCCCEepUhbYdjZ2REeHs6ePXuIiYmhf//+zJkzh2rVqhl8rgkTJrB48WJWrFjBmTNnGD16NKmpqQwfPhyAIUOGMGXKFF378ePHs2XLFubNm0dUVBQzZszg8OHDvPzyywC4ublRv359vZuFhQXe3t7Url27KG+3XGjq74KjtTl30rI4duWuqeMIIYQQZVqhB0EXJDU1lT179rB7927u3LlTpAJj4MCB3Lx5k2nTphEfH0+jRo3YsmWLbqBzbGwsavU/dVrr1q1ZvXo1b775JlOnTiUoKIiNGzdSv3794ryVcs/cTE37Wh5s+juOnVEJNPWXbUmEEEKIh1EpRVg+eO/evSxdupQffvgBRVF4+umnGTFiBG3atCmZlKUsOTkZJycnkpKSytV4oPVHrvLauuPU9XFk8/h2po4jhBBClCpDPr8L3QMUFxfHihUrWL58OTExMbRq1Yr58+fzzDPPYG9vb4zcopg61vZApYLTccnEJ6Xj7WRt6khCCCFEmVToAsjPzw83Nzeef/55RowYQZ06dUo2mTCYm70VDas6c+zKXXZGJzCoheFjsoQQQojKoNAF0Nq1a3niiScwNy/WsCFRwjoHe3Lsyl12REkBJIQQQjxMoWeB9e/fX4qfciBvd/h95xLJyM4xdRwhhBCiTCrSNHhRdtXzdcTTwYq0zBwOXrxt6jhCCCFEmSQFUAWjUqnolLsoouwOL4QQQhRMCqAKSLbFEEIIIR6tyAXQuXPn+P3337l//z4ARVhOSJSQtkHuWJipuHQrjQs3U0wdRwghhChzDC6Abt26RWhoKLVq1aJHjx7ExcUBMGLECF577bWSyCgMZG9lTovqriCXwYQQQogCGVwAvfrqq5ibmxMbG4utra3u+MCBA9myZYux84kiyhsHtDNaCiAhhBDi3wwugP744w/ef/99qlatqnc8KCiIy5cvGzObKIa86fAHL94mJSPb1HGEEEKIMsXgAig1NVWv5yfP7du3sbKyMlYuUUyBHvYEuNmSlaOw92yiqeMIIYQQZYrBBVC7du1YuXKl7r5KpUKj0TB37lw6depk7HyiGGQ2mBBCCFEwg5d2njt3Ll26dOHw4cNkZmby+uuvc+rUKW7fvs2+fftKJqUoks7Bnizbd4md0QkoioJKpTJ1JCGEEKJMMLgHqH79+sTExNC2bVv69OlDamoq/fv35+jRo9SoUaNkUooiaVHdFVtLMxLuZXDqerKp4wghhBBlRpE293JycuJ///uf8dMIo7IyN6NNTXe2nr7BjqgE6ldxMnUkIYQQokwwuAD6888/H/l4+/bti5NHGFnnYE9dATSuS5Cp4wghhBBlgsEFUMeOHfMde3BsSU6O7EBeluStB3T86l1upWTgZi8z9YQQQgiDxwDduXNH75aQkMCWLVto3rw5f/zxR8mkFEXm7WRNXR9HFAV2x9w0dRwhhBCiTDC4B8jJKf84kq5du2JpacmECRM4cuSIsbIJI+kc7MnpuGR2RCXQv0nVQjxDCCGEqNiMthu8l5cX0dHRxjqdMKK89YD+jLlJdo7G1HGEEEIIkzO4B+jvv//Wu68oCnFxcbz33ns0atTImNmEkTTyc8bF1oI7aVkcuXyHloFupo4khBBCmJTBBVCjRo1QqVQoiqJ3vFWrVixdutSY2YSRmKlVdKjlwcZj19kRnSAFkBBCiErP4ALo4sWLevfVajUeHh5YW1sbM5cwsk7Bnmw8dp2dUQlM6V7H1HGEEEIIkzJoDFBWVhbh4eFkZmbi7++Pv78/fn5+UvyUAx1qeaBWQcyNFK7eSTN1HCGEEMKkDCqALCws8o0BEuWDs60lTf1dQDZHFUIIIQyfBfbcc8+xZMmSkkkjSpRud/hoWQ9ICCFE5WbwGKDs7GyWLl3Ktm3baNq0KXZ2dnqPz58/35j5hBF1DvZk7pZo9p9PJD0rB2sLM1NHEkIIIUyi0AWQmZkZcXFxnDx5kiZNmgAQExOj1+bBLTFE2VPbywFfJ2uuJ6UTcf6WrkdICCGEqGwKXQDlTXvfuXNnSeYRJUilUtEp2JNvD8SyIypBCiAhhBCVltFWghblQ97mqDuiEvKt5SSEEEJUFgaNAfr666+xt7d/ZJtx48YVN5MoQa1rumFpruba3fucTUihlpeDqSMJIYQQpc6gAmjRokWYmT184KxKpZICqIyztTQnJNCN3TE32RGVIAWQEEKISsmgAujw4cN4esq4kfKuc7Anu2NusjMqgRc71DB1HCGEEKLUFXoMkMzwqjg65w5+Pnz5Dkn3s0wdRwghhCh1hS6AZMBsxeHnaktNT3tyNAp7zsqiiEIIISqfQhdA06dP/88B0KL8yOsF2iHbYgghhKiEDCqAbG1tSzaNKDUda3sAsDv6JhqN9O4JIYSoXGQdoEqqeYArDlbm3ErN5PjVu6aOI4QQQpQqKYAqKQszNe1quYPsDi+EEKISkgKoEtOtCh0tBZAQQojKRQqgSqxjbgF08loyCcnppo4jhBBClBqjFUBTp04lPDzcWKcTpcDDwYqGVZ0A2BUt0+GFEEJUHkYrgK5evcrFixeNdTpRSjrJdHghhBCVkNEKoJUrV7Jz505jnU6Ukrz1gPaeSyQzW2PqOEIIIUSpMLgAWrlyJRkZGfmOZ2ZmsnLlSmPlEqWkvq8T7vZWpGRkc+jSbVPHEUIIIUqFwQXQ8OHDSUpKynf83r17DB8+3Fi5RClRq1W6RRHlMpgQQojKwuACSFGUAjdGvXr1Kk5OTsbKJUpR3mUwWQ9ICCFEZWFe2IaNGzdGpVKhUqno0qUL5ub/PDUnJ4eLFy/SrVu3ksopSlDbIHfM1SouJKZyKTGVAHc7U0cSQgghSlShC6C+ffsCcOzYMcLCwvQ2RrW0tCQgIIAnn3yyZFKKEuVobUHzAFciLtxiZ3QCw92rmzqSEEIIUaIKXQBNnz4dgICAAAYOHIi1tXVJ5hKlrHOwJxEXbrEjKoHhbaQAEkIIUbEZPAZo6NChUvxUQHnrAR24cJvUjGxTxxFCCCFKlNHWARo6dCidO3c21ulEKavhYUc1V1syczTsO5do6jhCCCFEiTKoAFIUhdjYWNLT8+8bVaVKFfz9/Y2ZTZQilUpFp9zp8Dtlc1QhhBAVnMEFUM2aNbly5Uq+x2bPns2yZcuKFGLhwoUEBARgbW1Ny5YtOXjw4CPbr1u3juDgYKytrWnQoAGbN2/We3zGjBkEBwdjZ2eHi4sLoaGhHDhwoEjZKpNOuunwN1EUxdRxhBBCiBJjUAGkVqsJCgri1q1bRguwZs0aJkyYwPTp04mMjKRhw4aEhYWRkFBwL8T+/fsZNGgQI0aM4OjRo/Tt25e+ffty8uRJXZtatWqxYMECTpw4wd69ewkICODxxx/n5k3Z8PNRWgW6YWNhRnxyOqfjkk0dRwghhCgxKsXA/+r/8ssvzJ07ly+++IL69esXO0DLli1p3rw5CxYsAECj0eDn58fYsWOZPHlyvvYDBw4kNTWVTZs26Y61atWKRo0asWjRogJfIzk5GScnJ7Zt20aXLl3+M1Ne+6SkJBwdHYv1/sqbkSsOse1MApPCajOmU01TxxFCCCEKzZDPb4MHQQ8ZMoSDBw/SsGFDbGxscHV11bsZIjMzkyNHjhAaGvpPILWa0NBQIiIiCnxORESEXnuAsLCwh7bPzMzkq6++wsnJiYYNGxbYJiMjg+TkZL1bZSW7wwshhKgMCr0OUJ6PP/7YaC+emJhITk4OXl5eese9vLyIiooq8Dnx8fEFto+Pj9c7tmnTJp555hnS0tLw8fFh69atuLu7F3jOOXPmMHPmzGK/n4qgU21tAXQ09g53UjNxsbM0dSQhhBDC6AwugIYOHVoySYysU6dOHDt2jMTERBYvXsyAAQM4cOAAnp6e+dpOmTKFCRMm6O4nJyfj5+dXyonLBl9nG4K9HYiKv8fumJv0bVzF1JGEEEIIozPaOkBF4e7ujpmZGTdu3NA7fuPGDby9vQt8jre3d6Ha29nZUbNmTVq1asWSJUswNzdnyZIlBZ7TysoKR0dHvVtlJpfBhBBCVHRGK4BCQ0MJDAw06DmWlpY0bdqU7du3645pNBq2b99OSEhIgc8JCQnRaw+wdevWh7Z/8LwZGRkG5aus8naH3x1zk+wcjanjCCGEEEZn8CWwh+nXrx+JiYavIDxhwgSGDh1Ks2bNaNGiBR9//DGpqakMHz4ccgddV6lShTlz5gAwfvx4OnTowLx58+jZsyfff/89hw8f5quvvgIgNTWVd999lyeeeAIfHx8SExNZuHAh165d4+mnnzbW263QGvs542RjQdL9LI5euUvzAMMGtwshhBBlncEF0LRp0+jUqRMhISF6e4KNGTOmSAEGDhzIzZs3mTZtGvHx8TRq1IgtW7boBjrHxsaiVv/TUdW6dWtWr17Nm2++ydSpUwkKCmLjxo26KflmZmZERUWxYsUKEhMTcXNzo3nz5uzZs4d69eoVKWNlY26mpkMtD34+fp2dUQlSAAkhhKhwDF4HqGvXrkRERJCdnU3z5s3p0KEDHTt2pE2bNtjY2JRc0lJUmdcByrPx6DVeWXOMYG8HtrzS3tRxhBBCiP9UousAbd26lbt377J9+3Z69OjB4cOH6d+/P87OzrRt27Y4uUUZ0qGWB2oVRMXf4/rd+6aOI4QQQhhVkcYAmZub06ZNGzw8PHB1dcXBwYGNGzc+dO0eUf642FnSuJoLRy7fYWd0AoNbyka3QgghKg6De4C++uornn32WapUqULr1q3ZsmULbdu25fDhw7LXVgXTWbc5qkyHF0IIUbEY3AP04osv4uHhwWuvvcZLL72Evb19ySQTJtextgcf/B7NvnO3SM/KwdrCzNSRhBBCCKMwuAdow4YNDB48mO+//x4PDw9at27N1KlT+eOPP0hLSyuZlMIk6vo44u1ozf2sHP66cMvUcYQQQgijMbgHqG/fvvTt2xeApKQk9uzZw7p16+jVqxdqtZr09PSSyClMQKVS0SnYg+8OXmFnVAIda+ffRkQIIYQoj4o0CPrWrVvs3r2bXbt2sWvXLk6dOoWLiwvt2rUzfkJhUp1qe2oLoOibzFAUVCqVqSMJIYQQxWZwAdSgQQPOnDmDi4sL7du3Z9SoUXTo0IHHHnusZBIKk2pT0x1LMzWxt9M4fzOVmp4y5ksIIUT5V6RB0B06dNCtvCwqNjsrc1oGurLnbCI7oxKkABJCCFEhGDwIesyYMVL8VDKdZXd4IYQQFYzRdoMXFVen3MHPhy7dJjk9y9RxhBBCiGKTAkj8pwB3OwLd7cjWKOw9m2jqOEIIIUSxSQEkCqWTXAYTQghRgRhUAGVnZ/P2229z9erVkkskyqS8cUC7om+i0SimjiOEEEIUi0EFkLm5OR988AHZ2dkll0iUSc0DXLG3MicxJYOT15NMHUcIIYQoFoMvgXXu3Jndu3eXTBpRZlmaq2lb0x3kMpgQQogKwOB1gLp3787kyZM5ceIETZs2xc7OTu/xJ554wpj5RBnSOdiTLafi2RmVwCuhtUwdRwghhCgylaIoBg3oUKsf3mmkUqnIyckxRi6TSk5OxsnJiaSkJBwdHU0dp8xISE6nxeztABz6XygeDlamjiSEEELoGPL5bfAlMI1G89BbRSh+xMN5OlpTv4r2F2pXtFwGE0IIUX7JNHhhkM65iyLulAJICCFEOVakAmj37t307t2bmjVrUrNmTZ544gn27Nlj/HSizMlbD2hPTCJZORpTxxFCCCGKxOAC6JtvviE0NBRbW1vGjRvHuHHjsLGxoUuXLqxevbpkUooyo2FVZ9zsLLmXkc3hS3dMHUcIIYQoEoMHQdepU4cXXniBV199Ve/4/PnzWbx4MWfOnDF2xlIng6AfbcLaY2yIvMYL7QOZ2qOOqeMIIYQQUNKDoC9cuEDv3r3zHX/iiSe4ePGioacT5ZDsDi+EEKK8M7gA8vPzY/v27fmOb9u2DT8/P2PlEmVYuyAPzNQqziWkcOV2mqnjCCGEEAYzeCHE1157jXHjxnHs2DFat24NwL59+1i+fDmffPJJSWQUZYyTjQXN/F04cPE2O6ISGNo6wNSRhBBCCIMYXACNHj0ab29v5s2bx9q1ayF3XNCaNWvo06dPSWQUZVCnYE8pgIQQQpRbBhVA2dnZzJ49m/DwcPbu3VtyqUSZ1znYk/d+iyLiwi3SMrOxtTS4lhZCCCFMxuDd4OfOnSu7wQuCPO2p4mxDZraGiPO3TB1HCCGEMIjBg6C7dOkiu8ELVCqVzAYTQghRbslu8KLIOgd7suqvy+yMSkBRFFQqlakjCSGEEIUiu8EXQBZCLJz0rBwavf0H6VkatrzSjmBv+V4JIYQwHdkNXpQKawszWtdwB7kMJoQQopwxqADKysrC3NyckydPllwiUa50qu0BwE4pgIQQQpQjBhVAFhYWVKtWTXp6hE7e7vBHLt/hblqmqeMIIYQQhWLwJbD//e9/TJ06ldu3b5dMIlGuVHWxpZaXPRoFdsfcNHUcIYQQolAMngW2YMECzp07h6+vL/7+/vlmgUVGRhoznygHOgV7EnMjhV3RN+nTqIqp4wghhBD/yeACqG/fviWTRJRbnWt78uXuC+yKTiBHo2CmroDT4TUauB4J0Zsh+Tp0nAIu/qZOJYQQoogMLoCmT59eMklEudXU3wVHa3PupGVx7Mpdmvq7mDqScWTdhwu7tUVPzBZIufHPYzG/w1NLoUYnUyYUQghRRIUeA3Tw4MFHDn7OyMjQbY4qKhdzMzXta1WQ2WApNyFyFXz3LLxfHb4bCJErtMWPpQPU6we+jeH+bfimP+z7BAxbSksIIUQZUOiFEM3MzIiLi8PTUzvrx9HRkWPHjhEYGAjAjRs38PX1rRAzxGQhRMNtiLzKhLXHqevjyObx7Uwdp/AUBRJjtL080b/BlYPAA38lHKtCcA+o3R3824K5JWSlw6+vwbFvtG3q9YM+C8HS7qEvI4QQouQZ8vld6Etg/66TCqqbDFxUWlQgHWp5oFLB6bhk4pPS8XayNnWkh8vJhisHcouezXD7gv7jPo2gdm7R490A/r3Fh4U19FkAVRrDb5Ph1I9wMxoGfgNuNUr1rQghhCgag8cAPYrsBVV5udlb0bCqM8eu3GVndAKDWlQzdSR9Gffg3HZtL8/Z3+H+nX8eM7OE6u21BU+t7uBUiJlsKhU0Hwle9WHtEEg4DYs7wZNLIKhrib4VIYQQxWfUAkhUbp2DPTl25S4/HbtGgypOVHOzxdHawnSBkq5BzG/aoufin5DzwEKNNi5Qq5u26KnRGawcivYa1VrBC7u1RdDVg/Dt09D5f9BuYv6eIyGEEGWGQQXQ6dOniY+Ph9zLXVFRUaSkpACQmJhYMglFudE52JP5W2P468Jten22FwBnWwuqudrqbv5utvjlfu3jZGPcKfOKAvEn/rm0FXdc/3HXwNxLWz3AryWYGan+d/SBYb/Cljfg8FLY8Q5cPwb9FhW9sBJCCFGiCj0IWq1Wo1KpChznk3dcdoOv3BRF4aNtZ9l79iaxt++TmJLxyPYWZiqqumgLIv/cosgvt0iq5mqLnVUhCpTsTLi0R9vLE/0bJF994EEV+LX4p+hxDyr5XpkjK2DzRG1vk3tteOZb7esKIYQocYZ8fhe6ALp8+XKhXtzfv/wvDicFkHGkZmRz5U4asbfSiL2tvV2+lcaV22lcuZNGVs6jf/Xc7S11vUX+rv/0HAXYZeIR/yfq6M3acT2Z9/55koWt9pJW7e4QFAb2HiX/Rv/t6mFY8zzcuw5WjtDvS+1MMiGEECWqRAqgykQKoJKXo1GIT04nNrcgunw7ldjb94m9lUrs7TTupGXpta+mukFX9RFC1ZE0V0dhrtLoHrtn7sZ1rw7crx6GXZ3OVPVww8bSzATv6gEpCbB2KMTu197v8AZ0mAxqg7ffE0IIUUhSABWTFECml3w/g5tn9kP0Zlyubsc19bze42c0fmzTNGVbThP+VgJR/rWmp6eDlXbckVv+8Uce9lalM2MxJwv+eBMOLNLer9VN2xtk41zyry2EEJWQFEDFJAWQiWSmwcXcrSeit0DqA6tKq8wgoA3U7kF2zTDi1N5cfuDSWuztVN0ltnvp2Y98GRsLs3zjjfKKpaouNliZG7n36Nh3sOkVyE4H1xrwzGrwDDbuawghhJACqLikACpFKQnafbaif4PzOyH7/j+PWTlCzVAI7gk1u2inrv8HRVFIup+lK4Zib+deYsv9Oi7pPppH/MarVODjaP3P2KMHZq35u9nhYmtRtN6j60e144KSroClPfT9HOr2Mfw8QgghHkoKoGKSAqgEKYp21eS8qepXD+tvPeHk988qzP5ttFtPGFFmtoZrd+9re41upeoVSrG300jLfPQsRnsrc+r6OPJ233oEexv4u5GaCD8M165JBNB2AnR+E9QmHq8khBAVhBRAxSQFkJHlZENsRO5U9c1w56L+476N/yl6vOqbbAFBRVG4lZqZr9cobxZbfHK6rq2DlTlfPNeUtkHuhr1ITjZsmw4RC7T3a3SBJ78GW1cjvxshhKh8jF4ANW7cuNDd/pGRkYVPWkZJAWQE6clwPnfriZjfIf3uP4+ZWUFgh9ytJ7qBo68pkxZaelYOsbfTeGvjSQ5cvI25WsWc/g14upmf4Sc78QP89LL2kp+zv3a9IO8GJRFbCCEqDUM+vws1J7dv37706dOHPn36EBYWxvnz57GysqJjx4507NgRa2trzp8/T1hYWJECL1y4kICAAKytrWnZsiUHDx58ZPt169YRHByMtbU1DRo0YPPmzbrHsrKyeOONN2jQoAF2dnb4+voyZMgQrl+/XqRswgBZ6dqFAFf1h7mBsG4Y/L1GW/zYuELDZ2HAKnj9AgxeB83Cy03xA2BtYUYtLwdWjmhBn0a+ZGsUJv3wNx9vizF8I+AGT8HIrdri5+5l+LqrtigSQghRKgy+BDZy5Eh8fHyYNWuW3vHp06dz5coVli5dalCANWvWMGTIEBYtWkTLli35+OOPWbduHdHR0Xh6euZrv3//ftq3b8+cOXPo1asXq1ev5v333ycyMpL69euTlJTEU089xahRo2jYsCF37txh/Pjx5OTkcPjw4UJlkh4gA2WkaLeAiFgAKTf+Oe5aQ7sAYO2e2hWZK9BYF41GYd7WaBbu1E7Pf6ppVWb3a4CluYHr/KTdhvUj4PwO7f2QlyF0pvG26RBCiEqkRMcAOTk5cfjwYYKC9Jf3P3v2LM2aNSMpKcmgsC1btqR58+YsWKAdE6HRaPDz82Ps2LFMnjw5X/uBAweSmprKpk2bdMdatWpFo0aNWLRoUYGvcejQIVq0aMHly5epVu2/dymXAqiQ0m7Dwa+069zk7a7uWBWaj4A6vSvFFhCrD8Ty1k8nydEotK3pzufPNTF8A1hNjnb/sL3ztfert4enloOdW4lkFkKIisrol8AeZGNjw759+/Id37dvH9bW1gadKzMzkyNHjhAaGvpPILWa0NBQIiIiCnxORESEXnuAsLCwh7YHSEpKQqVS4exc8AJ0GRkZJCcn693EI9y7AX+8BR83gF1ztMWPaw14YgGMOwrtJlSK4gfg2ZbV+HpoM2wtzdh7LpEBiyKIS7pfiGc+QG0GodNhwEqwsNPOEvuqg3ZDVSGEECXC4H72V155hdGjRxMZGUmLFi0AOHDgAEuXLuWtt94y6FyJiYnk5OTg5eWld9zLy4uoqKgCnxMfH19g+7xd6v8tPT2dN954g0GDBj20GpwzZw4zZ840KHuldDcW9n0CkasgJ3ejU6/62oKnbt8KdYnLEJ1qe7L2/0IYvvwQUfH36LtwH8uGtaCur4G9h3X7gHst+H4w3D4PS8Og18fQaFBJRRdCiErL4B6gyZMns2LFCo4cOcK4ceMYN24ckZGRLFu2rMBLVqaUlZXFgAEDUBSFL7744qHtpkyZQlJSku525cqVUs1Z5t2MgR9Hw6eN4dDX2uKnanMYtAZe3Av1n6y0xU+e+lWc+PGl1gR52nMjOYMBX0bwZ8xNw0/kWQdG7dBu5JqdDhtfhM2va7fVEEIIYTRFGmk5YMAABgwYUOwXd3d3x8zMjBs3bugdv3HjBt7e3gU+x9vbu1Dt84qfy5cvs2PHjkdeC7SyssLKyqpY76VCijsOe+bB6Z//WawwsCO0ew0C2plsvZ6yqqqLLT+Mbs2Lq44QceEWw5cfYk6/BgxobuA0eRtnGPQ97H4Pdr8PB7+E+BMwYAXY558YIIQQwnBF2pr67t27fP3110ydOpXbt29D7vo/165dM+g8lpaWNG3alO3bt+uOaTQatm/fTkhISIHPCQkJ0WsPsHXrVr32ecXP2bNn2bZtG25uMpjUILF/wTdPwZft4fRP2uKndk8YuR2G/KQdpCvFT4GcbCxYEd6Cfo2rkKNReH3938z/I9rwafJqNXSaCs98B5YO2l3lv+wAV4+UVHQhhKhUDJ4F9vfffxMaGoqTkxOXLl0iOjqawMBA3nzzTWJjY1m5cqVBAdasWcPQoUP58ssvadGiBR9//DFr164lKioKLy8vhgwZQpUqVZgzZw7kToPv0KED7733Hj179uT7779n9uzZumnwWVlZPPXUU0RGRrJp0ya98UKurq5YWv731gqVchaYomgXLtwzHy7nDnJXqbWXt9q+Cl71TJ2wXFEUhflbY/hsxzkA+jeuwntPPmb4NHmAxLPw/bOQGANmltBzHjQZYvzQQghRzpXoNPjQ0FCaNGnC3LlzcXBw4Pjx4wQGBrJ//36effZZLl26ZHDgBQsW8MEHHxAfH0+jRo349NNPadmyJQAdO3YkICCA5cuX69qvW7eON998k0uXLhEUFMTcuXPp0aMHAJcuXaJ69eoFvs7OnTvp2LHjf+apVAWQRgNRm7SXuuJyZx2pLaDRs9D2FXANNHXCcm3NoVim/qidJt+6hhtfPNcUJxsDp8mTu7L2xtHanxVA0+HQ/X0wl0u3QgiRp8TXAYqMjKRGjRp6BdDly5epXbs26enphThL2VYpCqCcbDj5g7bHJzFae8zcBpoN1y7G51TF1AkrjN0xN3npmyOkZuZQy8ueZcNbUMXZxvATaTSwdx7seFd7WbJqC+3UeUefkogthBDlTomuA2RlZVXgOjkxMTF4eHgYejpR2rLS4dAS+KwJ/Ph/2uLHygnaTYRXT0K3OVL8GFmHWh6sfTEEL0crYm6k0G/hPk5eM2zBUMgdF9R+knYbEWsnuHpQu15Q7F8lEVsIISo0gwugJ554grfffpusLO20XJVKRWxsLG+88QZPPvlkSWQUxpCRAvs/g08awq8TtPtP2bpDl2nw6gno8hbYGbizuSi0er5O/PhSG2p7OZBwL4OBX0awMzqhaCcL6gqjdoJnXe3WI8t7apcnMHSgtRBCVGIGXwLL22vr8OHD3Lt3D19fX+Lj4wkJCWHz5s3Y2dmVXNpSUqEugd2/Awe+ggNfPLBdRRVoPU47kNbS1tQJK5Xk9CxGf3OEfeduYaZW8U7f+gxq8d/bsxQoIwV+fhlO/ai93/g56DEPLAxbkV0IISqKEh0DlGffvn0cP36clJQUmjRpkm97ivKsQhRAKQnazUkPLYHMFO0x10DtjK7HngHz/54NJ0pGZraGKRtOsD7yKgBjOtVg4uO1URVlaQFFgf2fwrYZoGjAtwkMXAVOVY0fXAghyrgSK4CysrKwsbHh2LFj1K9f3xhZy6RyXQDdjYV9n8LRVdqVhJHtKsoiRVH4eNtZPtl+FoC+jXx5/6nHsDIv4s/n/A74IVzby2frrl00MaCtcUMLIUQZV2KDoC0sLKhWrRo5OTnFzSiMLfEsbHwpd7uKxdriR7arKLNUKhWvdq3F3Kcew1ytYuOx6wxdepCktCJueVGjM7ywC7wbQFoirHgC/vpCxgUJIcRDGHwJbMmSJWzYsIFVq1bh6upacslMqFz1AMUd105lz1uxGaB6B2g/UbarKCf2nL3J6G8iScnIpqanPcuHN6eqSxHHZmWmwS/j4cRa7f0GA6D3JzLWSwhRKZToGKDGjRtz7tw5srKy8Pf3zzfoOTIysmipy5ByUQDF/qVdvPDsH/8cq91Du09X1WamTCaK4ExcMsOXHSI+OR0PByuWDm1Og6pORTuZosCBRfD7/0DJ0fYKDfwWXPyNHVsIIcoUQz6/Dd4MtW/fvsXJJopDUbRjPfbMh8t7tcdUaqjXXzvGR7arKLfq+Djy45jWDF92iKj4ewz8KoKFzzahU3ARNj9VqaDVaO3Yr3XDtBupftURnloKNTqVRHwhhCh3ijwLrCIrcz1AGg1E/6rt8bl+VHssb7uKNuPBrYapEwojuZeexUvfRrLnbCJqFczqW5/BLYvRc5N0FdY8p/29UakhdIZ2CQS5NCqEqIBKZRp8RVZmCqCcbDi5HvbOh5tR2mOyXUWFl5WjYeqGE6w7op0mP7pjDSY9Xhu1uohFS1Y6/PoaHPtGe79eP+izECzL/5pdQgjxoBItgHJycvjoo49Yu3YtsbGxZGZm6j1++/btoqUuQ0xeAGWlw/HVsPdj7YrNAFaO0OIF7aUNWbG5wlMUhc92nGP+1hgAejf05cOnizFNXlHg8BL47Q3QZGtXkR74jfQeCiEqlBLdC2zmzJnMnz+fgQMHkpSUxIQJE+jfvz9qtZoZM2YUJ7fISIH9C7TbVWx69V/bVZyU7SoqEZVKxbguQXz4dEPM1Sp+OX6d55cc5G5aZiGeXeAJoflIGLoJ7L0g4TQs7gRntxo7uhBClAsG9wDVqFGDTz/9lJ49e+Lg4MCxY8d0x/766y9Wr15dcmlLSan3AN2/AwcXa9dtuZ/bgybbVYhc+84l8uKqI9zLyCbQw44Vw1vg51qM34nkOFj7PFw9BKig8/+0m+HKuCAhRDlXoj1A8fHxNGjQAAB7e3uSkrS7Wvfq1Ytff/21qJkrp5QE2DodPmoAO9/VFj+ugfDEZzDuGLR6UYofQZua7qwbHYKPkzUXbqbS7/N9HL9yt+gndPSBYb9C0+HataN2vKMdKJ1xz5ixhRCiTDO4AKpatSpxcXGQ2xv0xx/adWgOHTqElZWV8RNWRHevwOZJ8HED2PcxZN4Dz3rw5BIYc0jb6yN7dYkHBHs78uNLbajj40hiSibPfPUXW0/fKPoJza2g98fQ+1Mws4SoTbC4i3ZFcSGEqAQMLoD69evH9u3bARg7dixvvfUWQUFBDBkyhPDw8JLIWHEknoWNY+DTRnDwK+12FVWawaDvtdtVNHgKzAxemklUEt5O1qx7MYT2tTy4n5XD/606zMqIS8U7adOhMPw3cPCFxGhY3BmiNhsrshBClFnFngYfERFBREQEQUFB9O7d23jJTKjExgBtmqCdiUPudhXtXoPq7WXshTBIVo6Gtzae5PtDVwB4oX0gk7sFF32aPLmXY9cOhdj92vtdpkHbCfK7KYQoV2QdoGIqsQLoziXYMhXavgp+zY13XlHpKIrCwp3n+PAP7TT5ng18mDegIdYWxdjwNicLfp+q7Z0EaDUGHn8H1AZ3FAshhEmUaAG0cuXKRz4+ZMgQQ05XJpl8HSAhCunHo1d5/Ye/ycpRaObvwuIhzXCxK+b4sYjP4fcp2q8bPqsdlC+XZoUQ5UCJFkAuLi5697OyskhLS8PS0hJbW1tZCFGIUrb/fCL/t+oI99KzCXS3Y9nw5vi7FXOV52PfwU9jtJup1u4BTy0DC2tjRRZCiBJRotPg79y5o3dLSUkhOjqatm3b8t133xUntxCiCFrXcGf96NZUcbbhQmIq/T/fz9HYO8U7aaNB2pWizawgejN88ySkJxsrshBCmJxRLu4HBQXx3nvvMX78eGOcTghhoFpeDvz4Umvq+TpyKzWTQYv/4o9T8cU7aXAPeH4DWDrA5b2wohek3DRWZCGEMCmjjW40Nzfn+vXrxjqdEMJAno7WrP2/EDrW9iA9S8P/fXOE5fsuFu+kAW1h2Cbtlixxx2FZN+06VkIIUc4ZPAbo559/1ruvKApxcXEsWLAAPz8/fvvtN2NnLHUyBkiUZ9k5Gt766RTfHYwFYGTb6kztUad40+QTz8GqvpB0RbtNy/M/gkdt44UWQggjKNFB0Op/TYlVqVR4eHjQuXNn5s2bh4+PT9FSlyFSAInyTlEUvth9nrlbogHoXt+bjwY2Kt40+aSrsKofJMaAjSs89wNUaWq80EIIUUyyDlAxSQEkKoqfjl1j0rq/yczR0KSaM18PbY5rcabJp96Cb5+C65FgaQ/PrIbADsaMLIQQRVais8CEEOVHn0ZVWDmiBY7W5kTG3qX/5/u4lJha9BPaucHQn7UrmGemaIuhM78YM7IQQpQKg3uAJkyYUOi28+fPL0omk5MeIFHRnEu4x9Clh7h29z6udpYsHtKMpv4uhXjmQ2Slw/oR2k1UVWrtpqpNnjdmZCGEMFiJXgLr1KkTR48eJSsri9q1tYMgY2JiMDMzo0mTJv+cWKVix44dRX0PJiUFkKiIEu6lM2L5YU5cS8LKXM0nzzSiW/1ijNnLyYZNr8DRVdr7XWdBm3FGyyuEEIYq0UtgvXv3pn379ly9epXIyEgiIyO5cuUKnTp1olevXuzcuZOdO3eW2+JHiIrK08Ga719oRZdgTzKyNYz+NpIle4sxTd7MXLtNRuvcomfrW7BtBsiwQiFEOWBwD1CVKlX4448/qFevnt7xkydP8vjjj1eItYCkB0hUZNk5Gmb8copv/tJOkx/eJoA3e9bFrDjT5Pd+pC1+AJoMhV4fgboYM86EEKIISrQHKDk5mZs3868Ge/PmTe7du2fo6YQQpczcTM2sPvWZ3D0YgGX7LvHSt0e4n5lT9JO2fRV6f6IdDxS5An4Ih+wM44UWQggjM7gA6tevH8OHD2fDhg1cvXqVq1evsn79ekaMGEH//v1LJqUQwqhUKhUvdqjBZ4MaY2mm5vdTNxi0+C9upRSjaGk6TLtpqtoCTm+E1QMhI8WYsYUQwmgMvgSWlpbGxIkTWbp0KVlZWZC7DcaIESP44IMPsLMr5i7UZYBcAhOVycGLtxm18jBJ97Pwd7Nl2bDmBHrYF/2E53fA989BVipUaQaD14GtqzEjCyFEgUplIcTU1FTOnz8PQI0aNSpE4ZNHCiBR2ZxLSGH48oNcuX0fF1sLZvdrwOP1vIs+LujqYe0aQffvgEewdusMR19jxxZCCD2luhL05cuXSU1NJTg4ON82GeWVFECiMrp5L4ORKw5x/GoSAH6uNgxrXZ0BzariYG1h+AkTzmi3zrgXB87V4PmN4FbD+MGFECJXiQyCXrp0ab6FDV944QUCAwNp0KAB9evX58oV2SVaiPLKw8GK715oxZhONXC2teDK7fvM2nSakDk7mPnLKS7fMnAFac86EP47uAbC3VhYGgZxf5dUfCGEMEihC6CvvvoKF5d/Vo7dsmULy5YtY+XKlRw6dAhnZ2dmzpxZUjmFEKXA1tKcSWHBREzuwux+DajpaU9KRjbL9l2i44e7GLXyMH9duEWhO45d/LVFkHcDSL0Jy3vC5f0l/TaEEOI/FfoSmJubG7t27aJBgwYAjB49mps3b/LDDz8AsGvXLoYPH87Fi8VYWK2MkEtgQmgpisKfZxNZuvciu2P+Wf6iro8j4W2r07uhD1bmhVjvJz0JVj8DsfvB3BoGrIRaYSUbXghR6ZTIJbD79+/rnWz//v20b99edz8wMJD4+PiiZhZClEEqlYoOtTxYEd6CbRPaM7hlNawt1JyOS2biuuO0eW8HH2+L4ea9/5g+b+0Ez62HoDDITofvBsHfa0vrbQghRD6FLoD8/f05cuQIAImJiZw6dYo2bdroHo+Pj8fJyalkUgohTK6mpwPv9mvAX1O68Ea3YLwdrUlMyeTjbWdp894OJq47zunryQ8/gaUtPPMtPDYQlBzYMAoOfFmab0EIIXTMC9tw6NChjBkzhlOnTrFjxw6Cg4Np2rSp7vH9+/dTv379ksophCgjnG0tGd2xBiPbVWfLyXiW7L3IsSt3+eHIVX44cpWQQDfC21anc7Bn/mn0ZhbQdxFYO8PBL+G317VT5Tu8AapibMUhhBAGKnQB9Prrr5OWlsaGDRvw9vZm3bp1eo/v27ePQYMGlURGIUQZZGGmpndDX3o39CUy9g5L917kt5PxRFy4RcSFW/i72TKsdQBPN/PD3uqBf2rUauj+Pti6wa7ZsGsOpN2Gbu9pHxNCiFJQ7HWAKiIZBC1E0Vy/e5+VEZf57mAsSfe1K8U7WJkzoLkfw1oH4Odqq/+EA1/Bb5O0XzcYAH0/1/YSCSFEEZTqQogVkRRAQhRPWmY2GyKvsXTfRS7c1K4fpFZB17pehLepTovqrqjyLnn9vQ42vgiabO0g6aeXa8cLCSGEgaQAKiYpgIQwDo1GYffZmyzde5E9ZxN1x+v5OhLepjq98qbRx/wOa4doZ4hVC4FB34ONs0mzCyHKHymAikkKICGML+bGPZbtu8SGyKtkZGsgd/Xp51v582zLarjfiszdQT4JvBrA8xvA3tPUsYUQ5YgUQMUkBZAQJedOaiarD8ayMuISN5K16wdZmqvp28iXF4PvE/jb85CaoN1C4/mN2tWkhRCiEKQAKiYpgIQoeVk5GjafiGPp3ou6DVgB+vmnMztlGjapV8HBR7uTvGcdk2YVQpQPJVoA5eTksHz5crZv305CQgIajUbv8R07dhQtdRkiBZAQpUdRFCJj7+ZOo49Do4And1hj8z7VlVgUaxdUz/0AVZuZOqoQoowz5PO70OsA5Rk/fjzLly+nZ8+e1K9f/5+ZHEIIUQQqlYqm/i409Xfh6p00VuVOo+97/02WWc6lSfo5Mpf24u4Ty/Bs1N3UcYUQFYTBPUDu7u6sXLmSHj16lFwqE5MeICFMKzUjmw2RV/lu7xkmJ79Le7MTZCpmLPd5k8bdhtHM30X+8yWEyKdENkPNY2lpSc2aNYuTTwghHsnOypznQwLY9Fo3NM98xwGb9liqchgZ9zbrF7/LEwv28ePRq2RmawpxNiGEyM/gHqB58+Zx4cIFFixYUGH/ByY9QEKUMZoc7q4bi/OZbwF4P+sZvsjpjaeDNUNC/BnUohpu9lamTimEMLES7QHau3cv3377LTVq1KB37970799f72aohQsXEhAQgLW1NS1btuTgwYOPbL9u3TqCg4OxtramQYMGbN68We/xDRs28Pjjj+Pm5oZKpeLYsWMGZxJClDFqM5wHLIS2EwB4w+J7ZtmuIeFeOh/+EUPr93Ywef3fRMffM3VSIUQ5YXAB5OzsTL9+/ejQoQPu7u44OTnp3QyxZs0aJkyYwPTp04mMjKRhw4aEhYWRkJBQYPv9+/czaNAgRowYwdGjR+nbty99+/bl5MmTujapqam0bduW999/39C3JoQoy1QqCJ0Oj78DwPOan9ldaz0Nfe3JyNbw/aErhH38J899fYAdUTfQaGSFDyHEw5l0HaCWLVvSvHlzFixYAIBGo8HPz4+xY8cyefLkfO0HDhxIamoqmzZt0h1r1aoVjRo1YtGiRXptL126RPXq1Tl69CiNGjV6ZI6MjAwyMjJ095OTk/Hz85NLYEKUVUe/gZ/HgqJBCe5FZLMPWXLgOltOxpNX91R3t2N4mwCebFIVOyuDJ7wKIcqhEr0EZiyZmZkcOXKE0NDQf8Ko1YSGhhIREVHgcyIiIvTaA4SFhT20fWHNmTNHrxfLz8+vWOcTQpSwxs/BgJVgZokqahNN973A50/VYvekToxqVx0Ha3MuJqYy7adTtJqzndmbz3Dt7n1TpxZClCFFKoB++OEHBgwYQKtWrWjSpInerbASExPJycnBy8tL77iXlxfx8fEFPic+Pt6g9oU1ZcoUkpKSdLcrV64U63xCiFJQpzcM/gEs7eHin7CiN35W9/lfz7pETOnCzCfqEeBmy730bL768wLt5+5kzLeRHLl8G1kAXwhhcAH06aefMnz4cLy8vDh69CgtWrTAzc2NCxcu0L17+VykzMrKCkdHR72bEKIcCOwAQ38GG1e4fhSWdYOkq9hbmTO0dQA7XuvIkqHNaFPTjRyNwq8n4njyiwj6LtzHT8eukZUj0+iFqKwMLoA+//xzvvrqKz777DMsLS15/fXX2bp1K+PGjSMpKakQZ9Byd3fHzMyMGzdu6B2/ceMG3t7eBT7H29vboPZCiEqgSlMI/x0cq0BiDCwJg8SzAKjVKrrU8eLbka3Y8ko7BjSriqW5muNXkxj//THavb+Tz3ed425apqnfhRCilBlcAMXGxtK6dWsAbGxsuHdPO+30+eef57vvviv0eSwtLWnatCnbt2/XHdNoNGzfvp2QkJACnxMSEqLXHmDr1q0PbS+EqCQ8ammLILeakHwVlnaD6/pLYAR7OzL3qYbsn9yZCV1r4W5vRXxyOnO3RBMyZwdvbjzB+ZspJnsLQojSZXAB5O3tze3btwGoVq0af/31FwAXL140+Lr6hAkTWLx4MStWrODMmTOMHj2a1NRUhg8fDsCQIUOYMmWKrv348ePZsmUL8+bNIyoqihkzZnD48GFefvllXZvbt29z7NgxTp8+DUB0dDTHjh0r9jghIUQZ5+ynLYJ8GkJaIizvBRf35Gvmbm/FuC5B7JvciQ+fbkgdH0fuZ+XwzV+xdJm3m+HLDrL3bKKMExKigjO4AOrcuTM///wzAMOHD+fVV1+la9euDBw4kH79+hl0roEDB/Lhhx8ybdo0GjVqxLFjx9iyZYtuoHNsbCxxcXG69q1bt2b16tV89dVXNGzYkB9++IGNGzdSv359XZuff/6Zxo0b07NnTwCeeeYZGjdunG+avBCiArJzh6GbIKAdZN6Db56EqM0FNrUyN+OpplXZPK4tq0e1JLSOJyoV7Iy+yXNLDtD9kz2sPXSF9KycUn8bQoiSZ/A6QBqNBo1Gg7m5dl2N77//nv379xMUFMT//d//YWlpWVJZS41shSFEOZeVDj8Mh+jNoDKDPguh0aD/fNrFxFSW7bvIusNXuZ9b+LjZWTK4lT/Pt/LHw0G22xCiLDPk89ukCyGWVVIACVEB5GRrF0s8vlp7P2w2hIwp1FOT0rL4/lAsK/Zf4npSOgCWZmqeaOTLiLbVqeMj/y4IURaVeAG0Z88evvzyS86fP88PP/xAlSpVWLVqFdWrV6dt27bFyV4mSAEkRAWh0cAfb8JfC7X3202Ezm9qt9UohKwcDVtOxrNk70WOXbmrO966hhsj2lanU21P1OqKuSm0EOVRia4EvX79esLCwrCxseHo0aO6LSSSkpKYPXt20VMLIYSxqdUQ9i50fkt7f8+H8OtroCncuB4LMzW9G/qycUwbNrzUmp6P+WCmVrH//C1GrDhMl/m7WRlxibTM7JJ9H0IIozO4B6hx48a8+uqrDBkyBAcHB44fP05gYCBHjx6le/fuFWK2lfQACVEBHVqiLX5QILAjBPeCaq3Asy6ozQp9mmt377Ni/yW+OxjLvXRt4eNobc6gltUYGhKAr7NNCb4JIcSjlOglMFtbW06fPk1AQIBeAXThwgXq1q1Lenp6cfObnBRAQlRQJ9fDhv8DTdY/x6wcoWoz8GsFfi20X1s5/OepUjOyWXf4Csv2X+LyrTQAzNQqejTwIbxNAI2ruZTkOxFCFMCQz2+Dt0j29vbm3LlzBAQE6B3fu3cvgYGBhqcVQojSUv9JcK8NZ36BKwfg6iHISIbzO7Q3AJUavOpre4f8Wmr/dKqa71R2VuYMa1Od50MC2BGVwJK9F/jrwm1+OX6dX45fp0k1Z0a0DSSsnhfmZibbd1oI8RAG9wDNmTOHb775hqVLl9K1a1c2b97M5cuXefXVV3nrrbcYO3ZsyaUtJdIDJEQlocmBG6e0xVDsX9o/kwrYDNmxyj/FkF9LbYFklv//jyevJbF030V+OX6drBztP61VnG0Y1jqAgS38cLS2KI13JUSlVaKXwBRFYfbs2cyZM4e0NG23r5WVFRMnTmTWrFnFS15GSAEkRCWWdE1bCOUVRfEnQPnXoGlLe+0eZNXyLps1B2sn3cMJyel889dlvjkQy+1U7T5jdpZmPN3Mj+FtAvB3syvtdyVEpVAq6wBlZmZy7tw5UlJSqFu3Lvb29kXNW+ZIASSE0MlIgWtH/imKrhyCjH9v/KwCr3r6vUTO1UjP1rDx6DWW7rtIzA3tPmMqFYTW8WJE2+q0rO6KqpBT8oUQ/00WQiwmKYCEEA+lyYGbUf9cMov9C+5ezt/OwUfbO+TXCsWvJXtTfFgScZVd0Td1Ter5OjKibXV6PeaLpbmMExKiuEqkAAoPDy/Uiy9durRwKcswKYCEEAa5F59bDB2AK39B3HHQ/GttIAtbqNKU226N+fm2H5+fcyEhyxYADwcrhrTyZ3Arf1ztyv92QkKYSokUQGq1Gn9/fxo3bvzIXZJ//PFHwxOXMVIACSGKJTMNrkfm9hId1BZH6XfzNbtlW4Pd6YHsTa/BYaU2N8y86d+kKuFtqhPk9d9T8YUQ+kqkABozZgzfffcd/v7+DB8+nOeeew5XV1djZS5TpAASQhiVRgOJMdreobxeotsX8jW7qThxWFOLI5pa5FRtQaeOXWkX7CvjhIQopBIbA5SRkcGGDRtYunQp+/fvp2fPnowYMYLHH3+8Qv0FlQJICFHiUhJye4e0RZESdwxVTqZek3TFghjzIMz8Qwhq2gXL6iFgWzH/4ymEMZTKIOjLly+zfPlyVq5cSXZ2NqdOnaowM8GkABJClLqsdLh+FK4cIO38foj9C9ucf882g2zXIMz9W+XONmsFbjUKvbmrEBVdia4EnUetVqNSqVAUhZycwm0sKIQQ4iEsrME/BPxDsG37CigKKdfPELl3C8kxe6mTdZoa6jjMb5+F22fh6Crt82zdtdPu/VpoiyKfRtpzCSEeqciXwPbu3UuvXr0YPnw43bp1Q62uOFM4pQdICFGWZOdo2Hr6Bmv/PIb62iGaqWNoqo6hkfoClmTpNzazhCrNoMMkqNHZVJGFMIkSuQT20ksv8f333+Pn50d4eDiDBw/G3d3dWJnLFCmAhBBl1fErd1my9yKbT8Sh1mRSX3WRUPtL9HSJxS/lBOq0f9YZol4/CJsDjj6mjCxEqSmxafDVqlWjcePGjxzwvGHDBsMTlzFSAAkhyrq4pPusjLjM6gOxJN3X9gI5WJnxYgMVQ8z/wOH4UlA0YOkAnaZCixcK3L9MiIqkRAqgYcOGFWqm17JlywqftIySAkgIUV6kZWazPvIay/Ze5EJiKgBmahXj66Xx4r2FWMZHaht61Yee86FaS9MGFqIEyVYYxSQFkBCivNFoFHbH3GTJ3ovsPZcIgI0FfFTzBI/HfYE6byHGJkMgdKZMpxcVkhRAxSQFkBCiPDt86TZzfoviyOU7AFS3SWOR9y/UjvtJ28DGFbrOhEbPQQWawCKEFEDFJAWQEKK8UxSFbWcSeH9LFOcStDvRd3e8yHtWy3G6d1bbqGoL6DUfvBuYNqwQRiIFUDFJASSEqCiyczRsiLzG/K0xxCenY042rzvvJjz7e8yzU0FlBi1fhE5TwEr2HxPlmxRAxSQFkBCioknPymH5/kt8vvMcyenZeHOLj53X0Cp9r7aBgw+EzdZOnZeVpUU5JQVQMUkBJISoqO6mZfLFrvMs23+JzGwNHdTH+cBuFZ5Z17UNanSGHh9qt9gQopyRAqiYpAASQlR01+7e56OtMayPvIqlkslo818YY/ELFkqmdjXptq9qbxY2po4qRKFJAVRMUgAJISqL6Ph7fPB7FNvOJBCgimOW5UraqY5rH3SpDj0+gKCupo4pRKFIAVRMUgAJISqbAxdu8d6WKI7G3qG7+iAzLFfhxW3tg3WegG5zwKmqqWMK8UhSABWTFEBCiMpIURR+P3WDub9HceNmIuPNNxBu/hvmaFAs7FB1nAytRoOZhamjClEgKYCKSQogIURllp2jYd2Rq3y0NQbXlLO8Y7GUZuoYABTPuqh6zgf/EFPHFCIfKYCKSQogIYSA+5k5LN13kS93nSUsewdTzFfjqtIuqkijwdD1bbBzN3VMIXQM+fyWNdCFEEIUyMbSjDGdarL79S44hQynW/ZHrM7upH3w2LfkfNoUDi8FjcbUUYUwmPQAFUB6gIQQIr+rd9KYvzWGi8d28Y75UuqpLwOQ5d0Eiyc+At9Gpo4oKjm5BFZMUgAJIcTDnYlL5sPfTuF3fjWvma/DQXUfDWqym47AsutbYO1k6oiikpICqJikABJCiP8Wcf4WX/26j343P+cJswgA0izdsewxG/OGA2RLDVHqZAyQEEKIEhdSw42lY3tjMWAZr1nP5LzGB9vMRMw3vsDNhWFoEqJNHVGIh5IeoAJID5AQQhgmK0fDugPnubttHuE5P2CtyiILc+LrjcKvzzSwtDV1RFEJyCWwYpICSAghiiYtM5sftu4h4OBM2quOAnDTzIu0LnPwb/2kqeOJCk4KoGKSAkgIIYrn1r10tm1cRrtzH+CrugXA3/ZtcH3yI6pWr23qeKKCkjFAQgghTMrNwZqBz49G89JBtrsOIksx47GUfbgtb8u2r94gMemeqSOKSk4KICGEECWmqpc7XcYt4tJTW4iyaoCNKpPQ64tImt+SH374jtSMbFNHFJWUXAIrgFwCE0KIEqAoxGz9Gs+IWTgrSQBsVrUntcMM+rZrjIWZ/J9cFI+MASomKYCEEKLkaFLvELt+CtUufI8ahWTFliVWzxHUcxw9GlRFrZb1g0TRyBggIYQQZZbazoWAIYvICd/GLce6OKrSeDXzK6qt783ET5ax/1yiqSOKSkAKICGEECZhUa0Zbq/sJePxuWSY2fOY+iIf3p3A+eX/x4tfb+fU9SRTRxQVmFwCK4BcAhNCiFKWkkD65ilYn/4BgJuKI7OzBqM0GMBrYcH4ucpCiuK/yRigYpICSAghTOTiHjJ/fhXLO2cBOKAJZkZOOK1atWVs5yBc7SxNnVCUYTIGSAghRPlUvR2WY/ZD6Aw0Zta0VEfxs/kUPA/MIWzubyzYcZa0TJk6L4pPeoAKID1AQghRBtyNhd8mQ/SvAFxT3JiZNYSjtm1oGeiGuVqFmVqNmRrM1Orc+9qbuVqFOvfP/PfVmKnAzCz3Oaq8NmBBNhZKFhZKJhZkYZ77tbmShQVZmOVkYq5kYqbJ+zMLM00GZkom6hztfbUmE3VOhvbP3K9VOZmoctJRZWdCTiZkp0N2hvaWkwHZucdyMrXHzCzBqar25uwHTrm3vK8dfMDM3NQ/oTJHLoEVkxRAQghRhkT/hvLb66juxgKwPacxf2iaYUUmlmRjSRZWqiysyM53TPc1WViqcv/Ue1x7zIpsrFRZpn6nhaaozFA5+uYWSXmFUVVwqvbP15Z2po5Z6qQAKiYpgIQQoozJTIM981D2fYJKUzqFShbmZGFBpsqS3P4gMrAkE3MysSADCzIVc9KxIEMxJ0OxIF2xIF0xJxPz3Mdz22GhPaY88HXuuTIeaJN3ThtVJr6qRHxVt6iqSqSKKhFftH/6qG5hqcr5z/waa1dUzlVROVfL7UH6V2+SnTuoKtaaS1IAFZMUQEIIUUbdjIG9H8H9O2Bu9c/N7N9fW4K5tfZSkrl17vG8ry312+vaPfB8M0tQF32YrEajkK1RyNEo5CgKOTkK2RqN9muNQnaOgkb5p43+fQ330rNJuJdBQnI6N5IzuJGczo3c+zfv3cdVc5eqqpv4qm5RJa9Ayv2ziioRR9X9/8yYY2ZNtkMV1M5+mLtUQ+Xsp9+b5FgFzCyK/D0wBSmAikkKICGEEGVVjkbhdmomN5LTSbj3QIGUnFsw3UsnJek2NmnX8eGfoujBm5fq7n++joKKdGtPshyqoHKuhqVbNSxd/XN7lHJ7k6wcSuU9F5Yhn99lYgTVwoUL+eCDD4iPj6dhw4Z89tlntGjR4qHt161bx1tvvcWlS5cICgri/fffp0ePHrrHFUVh+vTpLF68mLt379KmTRu++OILgoKCSukdCSGEECXDTK3Cw8EKDwcrwOmh7bJzNCSmZOYWR9oepJjcr28l3UNz9xoWqddxzIijCom5vUk3c3uTbmGlysYm/QY26TfgZiSczf8a980cSLPxIdO+Kiqnqli6+WPnVR0rN//cy2wexepJK0kmL4DWrFnDhAkTWLRoES1btuTjjz8mLCyM6OhoPD0987Xfv38/gwYNYs6cOfTq1YvVq1fTt29fIiMjqV+/PgBz587l008/ZcWKFVSvXp233nqLsLAwTp8+jbW1tQnepRBCCFG6zM3UeDtZ4+306M+9jOwcbt7L0PUg7U5O50byfdJux6HcvYJ5yjVs71/HLfvmP+ORVIk4q1KxybmHTco9SImB+PznzsSCJAtPUm18ybT3BSc/LN38sfUIwKlaXaxcqpbcN+A/mPwSWMuWLWnevDkLFiwAQKPR4Ofnx9ixY5k8eXK+9gMHDiQ1NZVNmzbpjrVq1YpGjRqxaNEiFEXB19eX1157jYkTJwKQlJSEl5cXy5cv55lnnsl3zoyMDDIyMnT3k5OT8fPzk0tgQgghRK60zGwSHhiPdOd2IhmJl9HcvYr5vavY3r+OY8YNvLmJryoRb+6gVj28xPjL42lajfnaqBnLzSWwzMxMjhw5wpQpU3TH1Go1oaGhREREFPiciIgIJkyYoHcsLCyMjRs3AnDx4kXi4+MJDQ3VPe7k5ETLli2JiIgosACaM2cOM2fONOI7E0IIISoWW0tzAtzNCXDPm17vCzym10ZRFFIysrmRnMGBu/dITrhMemIsmruxmCdfxeZ+HI4Zcbjn3CTTqbpJ3kcekxZAiYmJ5OTk4OXlpXfcy8uLqKioAp8THx9fYPv4+Hjd43nHHtbm36ZMmaJXVOX1AAkhhBCi8FQqFQ7WFjhYW1DT0x5q+QCt8rVTFAU/jWnnYJl8DFBZYGVlhZWVlaljCCGEEJWCSqXCwsy0axCZdGi2u7s7ZmZm3LhxQ+/4jRs38Pb2LvA53t7ej2yf96ch5xRCCCFE5WLSAsjS0pKmTZuyfft23TGNRsP27dsJCQkp8DkhISF67QG2bt2qa1+9enW8vb312iQnJ3PgwIGHnlMIIYQQlYvJL4FNmDCBoUOH0qxZM1q0aMHHH39Mamoqw4cPB2DIkCFUqVKFOXPmADB+/Hg6dOjAvHnz6NmzJ99//z2HDx/mq6++gtxutVdeeYV33nmHoKAg3TR4X19f+vbta9L3KoQQQoiyweQF0MCBA7l58ybTpk0jPj6eRo0asWXLFt0g5tjYWNQPLKLUunVrVq9ezZtvvsnUqVMJCgpi48aNujWAAF5//XVSU1N54YUXuHv3Lm3btmXLli2yBpAQQgghoCysA1QWyVYYQgghRPljyOd32VyfWgghhBCiBEkBJIQQQohKRwogIYQQQlQ6UgAJIYQQotKRAkgIIYQQlY4UQEIIIYSodKQAEkIIIUSlIwWQEEIIISodk68EXRblrQ2ZnJxs6ihCCCGEKKS8z+3CrPEsBVAB7t27B4Cfn5+powghhBDCQPfu3cPJyemRbWQrjAJoNBquX7+Og4MDKpXKqOdOTk7Gz8+PK1euVMhtNuT9lX8V/T3K+yv/Kvp7lPdXdIqicO/ePXx9ffX2ES2I9AAVQK1WU7Vq1RJ9DUdHxwr5i51H3l/5V9Hfo7y/8q+iv0d5f0XzXz0/eWQQtBBCCCEqHSmAhBBCCFHpSAFUyqysrJg+fTpWVlamjlIi5P2VfxX9Pcr7K/8q+nuU91c6ZBC0EEIIISod6QESQgghRKUjBZAQQgghKh0pgIQQQghR6UgBJIQQQohKRwqgUjBnzhyaN2+Og4MDnp6e9O3bl+joaFPHMqovvviCxx57TLewVUhICL/99pupY5WY9957D5VKxSuvvGLqKEYxY8YMVCqV3i04ONjUsYzu2rVrPPfcc7i5uWFjY0ODBg04fPiwqWMZRUBAQL6foUqlYsyYMaaOZhQ5OTm89dZbVK9eHRsbG2rUqMGsWbMKtedTeXLv3j1eeeUV/P39sbGxoXXr1hw6dMjUsYrkzz//pHfv3vj6+qJSqdi4caPe44qiMG3aNHx8fLCxsSE0NJSzZ8+WWj4pgErB7t27GTNmDH/99Rdbt24lKyuLxx9/nNTUVFNHM5qqVavy3nvvceTIEQ4fPkznzp3p06cPp06dMnU0ozt06BBffvkljz32mKmjGFW9evWIi4vT3fbu3WvqSEZ1584d2rRpg4WFBb/99hunT59m3rx5uLi4mDqaURw6dEjv57d161YAnn76aVNHM4r333+fL774ggULFnDmzBnef/995s6dy2effWbqaEY1cuRItm7dyqpVqzhx4gSPP/44oaGhXLt2zdTRDJaamkrDhg1ZuHBhgY/PnTuXTz/9lEWLFnHgwAHs7OwICwsjPT29dAIqotQlJCQogLJ7925TRylRLi4uytdff23qGEZ17949JSgoSNm6davSoUMHZfz48aaOZBTTp09XGjZsaOoYJeqNN95Q2rZta+oYpWb8+PFKjRo1FI1GY+ooRtGzZ08lPDxc71j//v2VwYMHmyyTsaWlpSlmZmbKpk2b9I43adJE+d///meyXMYAKD/++KPuvkajUby9vZUPPvhAd+zu3buKlZWV8t1335VKJukBMoGkpCQAXF1dTR2lROTk5PD999+TmppKSEiIqeMY1ZgxY+jZsyehoaGmjmJ0Z8+exdfXl8DAQAYPHkxsbKypIxnVzz//TLNmzXj66afx9PSkcePGLF682NSxSkRmZibffPMN4eHhRt/Q2VRat27N9u3biYmJAeD48ePs3buX7t27mzqa0WRnZ5OTk4O1tbXecRsbmwrXI3vx4kXi4+P1/i11cnKiZcuWRERElEoG2Qy1lGk0Gl555RXatGlD/fr1TR3HqE6cOEFISAjp6enY29vz448/UrduXVPHMprvv/+eyMjIcns9/lFatmzJ8uXLqV27NnFxccycOZN27dpx8uRJHBwcTB3PKC5cuMAXX3zBhAkTmDp1KocOHWLcuHFYWloydOhQU8czqo0bN3L37l2GDRtm6ihGM3nyZJKTkwkODsbMzIycnBzeffddBg8ebOpoRuPg4EBISAizZs2iTp06eHl58d133xEREUHNmjVNHc+o4uPjAfDy8tI77uXlpXuspEkBVMrGjBnDyZMnK1w1D1C7dm2OHTtGUlISP/zwA0OHDmX37t0Vogi6cuUK48ePZ+vWrfn+d1YRPPi/6Mcee4yWLVvi7+/P2rVrGTFihEmzGYtGo6FZs2bMnj0bgMaNG3Py5EkWLVpU4QqgJUuW0L17d3x9fU0dxWjWrl3Lt99+y+rVq6lXrx7Hjh3jlVdewdfXt0L9/FatWkV4eDhVqlTBzMyMJk2aMGjQII4cOWLqaBWOXAIrRS+//DKbNm1i586dVK1a1dRxjM7S0pKaNWvStGlT5syZQ8OGDfnkk09MHcsojhw5QkJCAk2aNMHc3Bxzc3N2797Np59+irm5OTk5OaaOaFTOzs7UqlWLc+fOmTqK0fj4+OQrxuvUqVPhLvVdvnyZbdu2MXLkSFNHMapJkyYxefJknnnmGRo0aMDzzz/Pq6++ypw5c0wdzahq1KjB7t27SUlJ4cqVKxw8eJCsrCwCAwNNHc2ovL29Abhx44be8Rs3bugeK2lSAJUCRVF4+eWX+fHHH9mxYwfVq1c3daRSodFoyMjIMHUMo+jSpQsnTpzg2LFjuluzZs0YPHgwx44dw8zMzNQRjSolJYXz58/j4+Nj6ihG06ZNm3zLT8TExODv72+yTCVh2bJleHp60rNnT1NHMaq0tDTUav2PLDMzMzQajckylSQ7Ozt8fHy4c+cOv//+O3369DF1JKOqXr063t7ebN++XXcsOTmZAwcOlNrYUbkEVgrGjBnD6tWr+emnn3BwcNBd33RycsLGxsbU8YxiypQpdO/enWrVqnHv3j1Wr17Nrl27+P33300dzSgcHBzyjdmys7PDzc2tQozlmjhxIr1798bf35/r168zffp0zMzMGDRokKmjGc2rr75K69atmT17NgMGDODgwYN89dVXfPXVV6aOZjQajYZly5YxdOhQzM0r1j/vvXv35t1336VatWrUq1ePo0ePMn/+fMLDw00dzah+//13FEWhdu3anDt3jkmTJhEcHMzw4cNNHc1gKSkper3IFy9e5NixY7i6ulKtWjVeeeUV3nnnHYKCgqhevTpvvfUWvr6+9O3bt3QClspcs0oOKPC2bNkyU0czmvDwcMXf31+xtLRUPDw8lC5duih//PGHqWOVqIo0DX7gwIGKj4+PYmlpqVSpUkUZOHCgcu7cOVPHMrpffvlFqV+/vmJlZaUEBwcrX331lakjGdXvv/+uAEp0dLSpoxhdcnKyMn78eKVatWqKtbW1EhgYqPzvf/9TMjIyTB3NqNasWaMEBgYqlpaWire3tzJmzBjl7t27po5VJDt37izws2/o0KGKkjsV/q233lK8vLwUKysrpUuXLqX6u6tSKtoymkIIIYQQ/0HGAAkhhBCi0pECSAghhBCVjhRAQgghhKh0pAASQgghRKUjBZAQQgghKh0pgIQQQghR6UgBJIQQQohKRwogIYQQQlQ6UgAJ8QjLly/H2dnZ1DEKZcaMGTRq1Mig56hUKjZu3Fhimcq6W7du4enpyaVLl0r1dTt27Mgrr7xS6Pbl6ffwYYry+/koly5dQqVScezYsWKdZ9GiRfTu3dtouUT5IQWQqNCGDRuGSqVCpVLpdqt/++23yc7ONnU0o5s4caLexoLG8OD3z8LCAi8vL7p27crSpUsrxCaU7777Ln369CEgIAAe+FA1MzPj2rVrem3j4uIwNzdHpVKVesFUWDt37qRHjx64ublha2tL3bp1ee211/K9l4rAz8+PuLg43V58u3btQqVScffuXYPOEx4eTmRkJHv27CmhpKKskgJIVHjdunUjLi6Os2fP8tprrzFjxgw++OADU8cyOnt7e9zc3Ix+3rzv36VLl/jtt9/o1KkT48ePp1evXiVeSGZmZpbYudPS0liyZAkjRozI91iVKlVYuXKl3rEVK1ZQpUqVEstTXF9++SWhoaF4e3uzfv16Tp8+zaJFi0hKSmLevHmmjmd0ZmZmeHt7F3vTV0tLS5599lk+/fRTo2UT5YMUQKLCs7KywtvbG39/f0aPHk1oaCg///wzAHfu3GHIkCG4uLhga2tL9+7dOXv2bIHnuXTpEmq1msOHD+sd//jjj/H390ej0ej+F7p9+3aaNWuGra0trVu3Jjo6Wu85X3zxBTVq1MDS0pLatWuzatUqvcdVKhVffvklvXr1wtbWljp16hAREcG5c+fo2LEjdnZ2tG7dmvPnz+ue8+9LDIcOHaJr1664u7vj5OREhw4diIyMLPL3r0qVKjRp0oSpU6fy008/8dtvv7F8+XJdu7t37zJy5Eg8PDxwdHSkc+fOHD9+XO9c77zzDp6enjg4ODBy5EgmT56sl3nYsGH07duXd999F19fX2rXrg3AlStXGDBgAM7Ozri6utKnT598vTBff/01derUwdramuDgYD7//PNHvq/NmzdjZWVFq1at8j02dOhQli1bpncsb5f1f9u9ezctWrTAysoKHx8fJk+erFcYpqamMmTIEOzt7fHx8SmwGMnIyGDixIlUqVIFOzs7WrZsya5dux6Z/0FXr15l3LhxjBs3jqVLl9KxY0cCAgJo3749X3/9NdOmTdO1Xb9+PfXq1cPKyoqAgIB8eQICAnjnnXd0mf39/fn555+5efMmffr0wd7enscee0zv70HeJbqNGzcSFBSEtbU1YWFhXLly5ZG5H/UzCw8P57HHHiMjIwNyi+HGjRszZMgQ+NclsEuXLtGpUycAXFxcUKlUDBs2jJUrV+Lm5qY7R56+ffvy/PPP6+737t2bn3/+mfv37xf6ey4qgFLbdlUIExg6dKjSp08fvWNPPPGE0qRJE93XderUUf7880/l2LFjSlhYmFKzZk0lMzNTURRFWbZsmeLk5KR7bteuXZWXXnpJ73yPPfaYMm3aNEV5YPfjli1bKrt27VJOnTqltGvXTmndurWu/YYNGxQLCwtl4cKFSnR0tDJv3jzFzMxM2bFjh64NoFSpUkVZs2aNEh0drfTt21cJCAhQOnfurGzZskU5ffq00qpVK6Vbt26650yfPl1p2LCh7v727duVVatWKWfOnFFOnz6tjBgxQvHy8lKSk5P1XufHH3806PuXp2HDhkr37t1190NDQ5XevXsrhw4dUmJiYpTXXntNcXNzU27duqUoiqJ88803irW1tbJ06VIlOjpamTlzpuLo6KiXeejQoYq9vb3y/PPPKydPnlROnjypZGZmKnXq1FHCw8OVv//+Wzl9+rTy7LPPKrVr19btBP7NN98oPj4+yvr165ULFy4o69evV1xdXZXly5c/9L2NGzdO7/unKIpy8eJFBVAOHjyouLu7K3v27FEURVH27NmjeHh4KAcPHlQA5eLFi4qiKMrVq1cVW1tb5aWXXlLOnDmj/Pjjj4q7u7syffp03TlHjx6tVKtWTdm2bZvy999/K7169VIcHByU8ePH69qMHDlSad26tfLnn38q586dUz744APFyspKiYmJUZQCfg//bf78+QqgXL9+/aFtFEVRDh8+rKjVauXtt99WoqOjlWXLlik2NjbKsmXLdG38/f0VV1dXZdGiRUpMTIwyevRoxdHRUenWrZuydu1a3e9jnTp1FI1Go8tnYWGhNGvWTNm/f79y+PBhpUWLFnq/9//+/fyvn9m9e/eUwMBA5ZVXXlEURVEmTpyoBAQEKElJSXo/q6NHjyrZ2dnK+vXrFUCJjo5W4uLilLt37yppaWmKk5OTsnbtWt3r3rhxQzE3N9f7+5aamqqo1Wpl586dj/z+iYpFCiBRoT34Aa7RaJStW7cqVlZWysSJE5WYmBgFUPbt26drn5iYqNjY2Oj+wfz3B8+aNWsUFxcXJT09XVEURTly5IiiUql0H4h5BdC2bdt0z/n1118VQLl//76iKIrSunVrZdSoUXo5n376aaVHjx66+4Dy5ptv6u5HREQogLJkyRLdse+++06xtrbW3f/3B8y/5eTkKA4ODsovv/yi9zpFLYAGDhyo1KlTR1FyCwRHR0fd9yVPjRo1lC+//FJRFEVp2bKlMmbMGL3H27Rpk68A8vLy0hU2iqIoq1atUmrXrq37sFUURcnIyFBsbGyU33//Xfc6q1ev1jv3rFmzlJCQkIe+tz59+ijh4eF6xx78UH3llVeU4cOHK4qiKMOHD1deffVV5ejRo3oF0NSpU/NlW7hwoWJvb6/k5OQo9+7dUywtLfU+gG/duqXY2NjoCqDLly8rZmZmyrVr1/SydOnSRZkyZYqiFKIAyitS/suzzz6rdO3aVe/YpEmTlLp16+ru+/v7K88995zuflxcnAIob731lu5Y3u9jXFycLh+g/PXXX7o2Z86cUQDlwIEDilLA72dhfmb79+9XLCwslLfeeksxNzfXFaTKv35WygN/9+7cuZPve/NgoT5v3jwlMDBQ72emKIri4uLyyIJZVDxyCUxUeJs2bcLe3h5ra2u6d+/OwIEDmTFjBmfOnMHc3JyWLVvq2rq5uVG7dm3OnDlT4Ln69u2LmZkZP/74I+R2/Xfq1Ek3iDbPY489pvvax8cHgISEBADOnDlDmzZt9Nq3adMm32s+eA4vLy8AGjRooHcsPT2d5OTkArPeuHGDUaNGERQUhJOTE46OjqSkpBAbG/uf37PCUBQFlUoFwPHjx0lJScHNzQ17e3vd7eLFi7rLdNHR0bRo0ULvHP++n/ceLS0tdfePHz/OuXPncHBw0J3X1dWV9PR0zp8/T2pqKufPn2fEiBF6r/3OO+/oXSL8t/v372Ntbf3Qx8PDw1m3bh3x8fGsW7eO8PDwfG3OnDlDSEiI7vtA7s8yJSWFq1evcv78eTIzM/V+x1xdXXWX9gBOnDhBTk4OtWrV0su/e/fuR+Z/0IM/i0d52O/e2bNnycnJ0R0rzO8eD/xOA5ibm9O8eXPd/eDgYJydnQv8u1TYn1lISAgTJ05k1qxZvPbaa7Rt27ZQ348HjRo1ij/++EM3EHz58uW6wf0PsrGxIS0tzeDzi/KreKPHhCgHOnXqxBdffIGlpSW+vr7FGjRpaWnJkCFDWLZsGf3792f16tV88skn+dpZWFjovs77h9bQWVMFncOQ8w4dOpRbt27xySef4O/vj5WVFSEhIUYbWHzmzBmqV68OQEpKCj4+PgWOWzF0+radnZ3e/ZSUFJo2bcq3336br62HhwcpKSkALF68WK/QIHeg7MO4u7tz586dhz7eoEEDgoODGTRoEHXq1KF+/frFnnJdkJSUFMzMzDhy5Ei+vPb29oU6R61atUhKSiIuLk5XcBdHcX/3/kthf2YajYZ9+/ZhZmbGuXPnivRajRs3pmHDhqxcuZLHH3+cU6dO8euvv+Zrd/v2bTw8PIr0GqJ8kh4gUeHZ2dlRs2ZNqlWrplf81KlTh+zsbA4cOKA7duvWLaKjo6lbt+5Dzzdy5Ei2bdvG559/TnZ2Nv379zcoT506ddi3b5/esX379j3yNYti3759jBs3jh49eugGvSYmJhrl3Dt27ODEiRM8+eSTADRp0oT4+HjMzc2pWbOm3s3d3R2A2rVrc+jQIb3z/Pt+QZo0acLZs2fx9PTMd24nJye8vLzw9fXlwoUL+R7PK9AK0rhxY06fPv3I1w4PD2fXrl0F9v6Q+7OMiIhAezVRa9++fTg4OFC1alVq1KiBhYWF3u/YnTt3iImJ0cuRk5NDQkJCvvze3t7/+f0BeOqpp7C0tGTu3LkFPp43Nfxhv3u1atV6ZLFYGNnZ2XoDo6Ojo7l79y516tTJ17awP7MPPviAqKgodu/ezZYtW/INTH9QXq/hgz1ZeUaOHMny5ctZtmwZoaGh+Pn56T1+/vx50tPTady4cZHfvyh/pAASlVZQUBB9+vRh1KhR7N27l+PHj/Pcc89RpUoV+vTp89Dn1alTh1atWvHGG28waNAgbGxsDHrdSZMmsXz5cr744gvOnj3L/Pnz2bBhAxMnTjTCu/pHUFAQq1at4syZMxw4cIDBgwcbnJXcGUrx8fFcu3aNyMhIZs+eTZ8+fejVq5duRk5oaCghISH07duXP/74g0uXLrF//37+97//6T4Ux44dy5IlS1ixYgVnz57lnXfe4e+///7PSzeDBw/G3d2dPn36sGfPHi5evMiuXbsYN24cV69eBWDmzJnMmTOHTz/9lJiYGE6cOMGyZcuYP3/+Q88bFhbGqVOnHtkLNGrUKG7evMnIkSMLfPyll17iypUrjB07lqioKH766SemT5/OhAkTUKvV2NvbM2LECCZNmsSOHTs4efIkw4YNQ63+55/eWrVqMXjwYIYMGcKGDRu4ePEiBw8eZM6cOQX2VBTEz8+Pjz76iE8++YQRI0awe/duLl++zL59+/i///s/Zs2aBcBrr73G9u3bmTVrFjExMaxYsYIFCxYY5XfPwsKCsWPHcuDAAY4cOcKwYcNo1apVgZc5KcTP7OjRo0ybNo2vv/6aNm3aMH/+fMaPH8+FCxcKPJ+/vz8qlYpNmzZx8+ZNXS8TwLPPPsvVq1dZvHhxgcXsnj17CAwMpEaNGsX+PohyxNSDkIQoSY8axKsoinL79m3l+eefV5ycnBQbGxslLCxMN/NGecTg0yVLluhmCz2ooIGY/x44qyiK8vnnnyuBgYGKhYWFUqtWLWXlypV65/n34OR/D/gs6LX+Pcg0MjJSadasmWJtba0EBQUp69atU/z9/ZWPPvrooa9T0PcPUADF3Nxc8fDwUEJDQ5WlS5cqOTk5em2Tk5OVsWPHKr6+voqFhYXi5+enDB48WImNjdW1efvttxV3d3fF3t5eCQ8PV8aNG6e0atVK7/UK+nnFxcUpQ4YMUdzd3RUrKyslMDBQGTVqlG5GkKIoyrfffqs0atRIsbS0VFxcXJT27dsrGzZseOh7UxRFadGihbJo0aJHfp8fVNDPcteuXUrz5s0VS0tLxdvbW3njjTeUrKws3eP37t1TnnvuOcXW1lbx8vJS5s6dq3To0EFvFlhmZqYybdo0JSAgQLGwsFB8fHyUfv36KX///beiFGIQdJ6tW7cqYWFhiouLi2Jtba0EBwcrEydO1Jsd9sMPPyh169ZVLCwslGrVqikffPCB3jn+/TuiFOL3MS/f+vXrlcDAQMXKykoJDQ1VLl++rHtOQYP0H/Yzu3//vlK3bl3lhRde0Gv/xBNPKK1bt1ays7ML/Fm9/fbbire3t6JSqZShQ4fqPff5559XXF1d8w3UVxRFefzxx5U5c+b85/dXVCwq5cG+WyFEocyaNYt169bx999/mzpKuda1a1e8vb3zrYNUWn799VcmTZrEyZMn9XplhGGWL1/OK6+8YvAqzKWpS5cu1KtXL9+Ch6dOnaJz587ExMTg5ORksnyi9MkgaCEMkJKSwqVLl1iwYAHvvPOOqeOUK2lpaSxatIiwsDDMzMz47rvv2LZtG1u3bjVZpp49e3L27FmuXbuWb1yIqBju3LnDrl272LVrV4GLY8bFxbFy5UopfiohKYCEMMDLL7/Md999R9++fR86MFYUTKVSsXnzZt59913S09OpXbs269evJzQ01KS5DNmUVJQ/jRs35s6dO7z//vt6yw/kMfXvnzAduQQmhBBCiEpHLnoLIYQQotKRAkgIIYQQlY4UQEIIIYSodKQAEkIIIUSlIwWQEEIIISodKYCEEEIIUelIASSEEEKISkcKICGEEEJUOv8PChtoz02/VUQAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"fig, ax = plt.subplots()\n",
"\n",
"ax.plot(df_results[\"degree\"], df_results[\"mse_train\"], label=\"Training MSE\")\n",
"ax.plot(df_results[\"degree\"], df_results[\"mse_test\"], label=\"Testing MSE\")\n",
"ax.set_xlabel(\"Polynomial Degree (Model Complexity)\")\n",
"ax.set_ylabel(\"Mean Squared Error w.r.t. True Values\")\n",
"ax.legend()"
]
},
{
"cell_type": "markdown",
"id": "5e5b5954",
"metadata": {},
"source": [
"**f)** Interpret the graph. Why do the lines move as they do? What does it tell us about model performance and generalizability?"
]
},
{
"cell_type": "markdown",
"id": "ad2acfb9",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-success\">\n",
"According to the graph the models ability to generalize increases with the degree of the polynomial. This is indicated by the testing MSE decreasing with the increase in polynomial degree. The difference to the graph from Bickel et al stems from the fact, that we use an exponential function which in turn can be written as an infinite polynomial function. Thus we approach the source of the data with an increasing degree of the polynomial inputs. Thus in this case the models generalization ability increases with the degree. We don't see the typical features of overfitting normally arriving with an increase in parameters of the model since our ratios between dataset size, noise in the training data and parameter number is still in favor of good generalizbility. With a further increase in noise or parameter number or a decrease in dataset size we could observe a decrease in model performance with resepect to the testing set. Our model would be overfit.\n",
"\n",
"</div>"
]
},
{
"cell_type": "markdown",
"id": "5994f0c5",
"metadata": {},
"source": [
"## Exercise 5 - Comparing your code with sklearn"
]
},
{
"cell_type": "markdown",
"id": "8f595b7a",
"metadata": {},
"source": [
"When implementing different algorithms for the first time, it can be helpful to double check your results with established implementations before you go on to add more complexity."
]
},
{
"cell_type": "markdown",
"id": "8ab310c1",
"metadata": {},
"source": [
"**a)** Make sure your `polynomial_features` function creates the same feature matrix as sklearns PolynomialFeatures.\n",
"\n",
"(https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "85b964d1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sklearn.preprocessing import PolynomialFeatures\n",
"\n",
"poly_features_sklearn = PolynomialFeatures(degree=5, include_bias=True).fit_transform(x.reshape(-1, 1))\n",
"poly_features_own = polynomial_features(x, 5)\n",
"\n",
"np.allclose(poly_features_sklearn, poly_features_own) # If this is true, our output is identical (within numerical precision)"
]
},
{
"cell_type": "markdown",
"id": "73c32c52",
"metadata": {},
"source": [
"**b)** Make sure your `OLS_parameters` function computes the same parameters as sklearns LinearRegression with fit_intercept set to False, since the intercept is included in the feature matrix. Use `your_model_object.coef_` to extract the computed parameters.\n",
"\n",
"(https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "35b04126",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from sklearn.linear_model import LinearRegression\n",
"\n",
"model = LinearRegression(fit_intercept=False)\n",
"model.fit(poly_features_sklearn, y)\n",
"\n",
"beta_sklearn = model.coef_\n",
"beta_own = OLS_parameters(poly_features_own, y)\n",
"\n",
"np.allclose(beta_sklearn, beta_own) # If this is true, our coefficients are identical (within numerical precision)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "lecture-materials",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Reference in New Issue View Git Blame Copy Permalink