diff --git a/doc/pub/week36/ipynb/week36.ipynb b/doc/pub/week36/ipynb/week36.ipynb index bdaee795a..f76a80988 100644 --- a/doc/pub/week36/ipynb/week36.ipynb +++ b/doc/pub/week36/ipynb/week36.ipynb @@ -3,9 +3,7 @@ { "cell_type": "markdown", "id": "7411173e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "\n", @@ -15,9 +13,7 @@ { "cell_type": "markdown", "id": "f50db913", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "# Week 36: Statistical interpretation of Linear Regression and Resampling techniques\n", "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n", @@ -28,9 +24,7 @@ { "cell_type": "markdown", "id": "4a90e380", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Plans for week 36\n", "\n", @@ -54,9 +48,7 @@ { "cell_type": "markdown", "id": "351d36fd", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Material for the active learning sessions Tuesday and Wednesday\n", "\n", @@ -66,9 +58,7 @@ { "cell_type": "markdown", "id": "0f9f4d2f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Linear Regression and the SVD\n", "\n", @@ -78,9 +68,7 @@ { "cell_type": "markdown", "id": "9ec36bdd", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}^T\\boldsymbol{X}=\\boldsymbol{V}\\boldsymbol{\\Sigma}^T\\boldsymbol{U}^T\\boldsymbol{U}\\boldsymbol{\\Sigma}\\boldsymbol{V}^T=\\boldsymbol{V}\\boldsymbol{\\Sigma}^T\\boldsymbol{\\Sigma}\\boldsymbol{V}^T.\n", @@ -90,9 +78,7 @@ { "cell_type": "markdown", "id": "614493c6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Since the matrices here have dimension $p\\times p$, with $p$ corresponding to the singular values, we defined last week the matrix" ] @@ -100,9 +86,7 @@ { "cell_type": "markdown", "id": "dc4be423", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{\\Sigma}^T\\boldsymbol{\\Sigma} = \\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} & \\boldsymbol{0}\\\\ \\end{bmatrix}\\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} \\\\ \\boldsymbol{0}\\end{bmatrix},\n", @@ -112,9 +96,7 @@ { "cell_type": "markdown", "id": "b262332c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "where the tilde-matrix $\\tilde{\\boldsymbol{\\Sigma}}$ is a matrix of dimension $p\\times p$ containing only the singular values $\\sigma_i$, that is" ] @@ -122,9 +104,7 @@ { "cell_type": "markdown", "id": "ab1829a4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\tilde{\\boldsymbol{\\Sigma}}=\\begin{bmatrix} \\sigma_0 & 0 & 0 & \\dots & 0 & 0 \\\\\n", @@ -139,9 +119,7 @@ { "cell_type": "markdown", "id": "7c82dd9e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "meaning we can write" ] @@ -149,9 +127,7 @@ { "cell_type": "markdown", "id": "5bab8721", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}^T\\boldsymbol{X}=\\boldsymbol{V}\\tilde{\\boldsymbol{\\Sigma}}^2\\boldsymbol{V}^T.\n", @@ -161,9 +137,7 @@ { "cell_type": "markdown", "id": "72273a4f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Multiplying from the right with $\\boldsymbol{V}$ (using the orthogonality of $\\boldsymbol{V}$) we get" ] @@ -171,9 +145,7 @@ { "cell_type": "markdown", "id": "d4af03ea", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)\\boldsymbol{V}=\\boldsymbol{V}\\tilde{\\boldsymbol{\\Sigma}}^2.\n", @@ -183,9 +155,7 @@ { "cell_type": "markdown", "id": "c0d35853", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## What does it mean?\n", "\n", @@ -197,9 +167,7 @@ { "cell_type": "markdown", "id": "db5ea827", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)\\boldsymbol{v}_i=\\boldsymbol{v}_i\\sigma_i^2.\n", @@ -209,9 +177,7 @@ { "cell_type": "markdown", "id": "a6b341ce", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "In other words, each non-zero singular value of $\\boldsymbol{X}$ is a positive\n", "square root of an eigenvalue of $\\boldsymbol{X}^T\\boldsymbol{X}$. It means also that\n", @@ -231,9 +197,7 @@ { "cell_type": "markdown", "id": "db723709", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{C}[\\boldsymbol{X}]=\\frac{1}{n}\\boldsymbol{X}^T\\boldsymbol{X},\n", @@ -243,9 +207,7 @@ { "cell_type": "markdown", "id": "1e600268", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "meaning that every squared non-singular value of $\\boldsymbol{X}$ divided by $n$ (\n", "the number of samples) are the eigenvalues of the covariance\n", @@ -258,9 +220,7 @@ { "cell_type": "markdown", "id": "b9452938", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## And finally $\\boldsymbol{X}\\boldsymbol{X}^T$\n", "\n", @@ -270,9 +230,7 @@ { "cell_type": "markdown", "id": "53fea1e6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}\\boldsymbol{X}^T=\\boldsymbol{U}\\boldsymbol{\\Sigma}\\boldsymbol{V}^T\\boldsymbol{V}\\boldsymbol{\\Sigma}^T\\boldsymbol{U}^T=\\boldsymbol{U}\\boldsymbol{\\Sigma}^T\\boldsymbol{\\Sigma}\\boldsymbol{U}^T.\n", @@ -282,9 +240,7 @@ { "cell_type": "markdown", "id": "ed10baad", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Since the matrices here have dimension $n\\times n$, we have" ] @@ -292,9 +248,7 @@ { "cell_type": "markdown", "id": "a4d0db96", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{\\Sigma}\\boldsymbol{\\Sigma}^T = \\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} \\\\ \\boldsymbol{0}\\\\ \\end{bmatrix}\\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} \\boldsymbol{0}\\\\ \\end{bmatrix}=\\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} & \\boldsymbol{0} \\\\ \\boldsymbol{0} & \\boldsymbol{0}\\\\ \\end{bmatrix},\n", @@ -304,9 +258,7 @@ { "cell_type": "markdown", "id": "fae3cadc", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "leading to" ] @@ -314,9 +266,7 @@ { "cell_type": "markdown", "id": "e973c775", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}\\boldsymbol{X}^T=\\boldsymbol{U}\\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} & \\boldsymbol{0} \\\\ \\boldsymbol{0} & \\boldsymbol{0}\\\\ \\end{bmatrix}\\boldsymbol{U}^T.\n", @@ -326,9 +276,7 @@ { "cell_type": "markdown", "id": "f88a33f4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Multiplying with $\\boldsymbol{U}$ from the right gives us the eigenvalue problem" ] @@ -336,9 +284,7 @@ { "cell_type": "markdown", "id": "5bc7214d", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "(\\boldsymbol{X}\\boldsymbol{X}^T)\\boldsymbol{U}=\\boldsymbol{U}\\begin{bmatrix} \\tilde{\\boldsymbol{\\Sigma}} & \\boldsymbol{0} \\\\ \\boldsymbol{0} & \\boldsymbol{0}\\\\ \\end{bmatrix}.\n", @@ -348,9 +294,7 @@ { "cell_type": "markdown", "id": "561d4c91", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "It means that the eigenvalues of $\\boldsymbol{X}\\boldsymbol{X}^T$ are again given by\n", "the non-zero singular values plus now a series of zeros. The column\n", @@ -365,9 +309,7 @@ { "cell_type": "markdown", "id": "0eda2d1e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Code for SVD and Inversion of Matrices\n", "\n", @@ -379,10 +321,7 @@ "cell_type": "code", "execution_count": 1, "id": "c95dbc32", - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -393,9 +332,7 @@ { "cell_type": "markdown", "id": "77079d25", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Let us first look at a matrix which does not causes problems and write our own function where we just use the SVD." ] @@ -404,10 +341,7 @@ "cell_type": "code", "execution_count": 2, "id": "cd0ecc44", - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -444,9 +378,7 @@ { "cell_type": "markdown", "id": "41d896dc", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Inverse of Rectangular Matrix\n", "\n", @@ -464,9 +396,7 @@ { "cell_type": "markdown", "id": "f594852c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{A}_{\\mathrm{PI}}= \\boldsymbol{V}\\boldsymbol{D}_{\\mathrm{PI}}\\boldsymbol{U}^T,\n", @@ -476,9 +406,7 @@ { "cell_type": "markdown", "id": "b78264f2", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "where $\\boldsymbol{D}_{\\mathrm{PI}}$ can be calculated by creating a diagonal matrix from $\\boldsymbol{\\Sigma}$ where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD." ] @@ -487,10 +415,7 @@ "cell_type": "code", "execution_count": 3, "id": "d05eff95", - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -521,9 +446,7 @@ { "cell_type": "markdown", "id": "2fae32c6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "As you can see from this example, our own decomposition based on the SVD agrees with the pseudoinverse algorithm provided by **Numpy**." ] @@ -531,9 +454,7 @@ { "cell_type": "markdown", "id": "8631011a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Ridge and LASSO Regression\n", "\n", @@ -544,9 +465,7 @@ { "cell_type": "markdown", "id": "8e7701de", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in {\\mathbb{R}}^{p}}}\\frac{1}{n}\\left\\{\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right)^T\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right)\\right\\}.\n", @@ -556,9 +475,7 @@ { "cell_type": "markdown", "id": "aea7230b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "or we can state it as" ] @@ -566,9 +483,7 @@ { "cell_type": "markdown", "id": "c5029213", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n", @@ -579,9 +494,7 @@ { "cell_type": "markdown", "id": "535ad70e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "where we have used the definition of a norm-2 vector, that is" ] @@ -589,9 +502,7 @@ { "cell_type": "markdown", "id": "159ae526", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\vert\\vert \\boldsymbol{x}\\vert\\vert_2 = \\sqrt{\\sum_i x_i^2}.\n", @@ -601,9 +512,7 @@ { "cell_type": "markdown", "id": "2f25dfa4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## From OLS to Ridge and Lasso\n", "\n", @@ -616,9 +525,7 @@ { "cell_type": "markdown", "id": "9d187806", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n", @@ -629,9 +536,7 @@ { "cell_type": "markdown", "id": "6c692db5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which leads to the Ridge regression minimization problem where we\n", "require that $\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_2^2\\le t$, where $t$ is\n", @@ -643,9 +548,7 @@ { "cell_type": "markdown", "id": "3404e0b4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{X},\\boldsymbol{\\beta})=\\frac{1}{n}\\vert\\vert \\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\vert\\vert_2^2+\\lambda\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_1,\n", @@ -655,9 +558,7 @@ { "cell_type": "markdown", "id": "4bc49666", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "we have a new optimization equation" ] @@ -665,9 +566,7 @@ { "cell_type": "markdown", "id": "3f08dfa1", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n", @@ -678,9 +577,7 @@ { "cell_type": "markdown", "id": "30bcc0cd", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which leads to Lasso regression. Lasso stands for least absolute shrinkage and selection operator. \n", "\n", @@ -690,9 +587,7 @@ { "cell_type": "markdown", "id": "d311a7d8", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\vert\\vert \\boldsymbol{x}\\vert\\vert_1 = \\sum_i \\vert x_i\\vert.\n", @@ -702,9 +597,7 @@ { "cell_type": "markdown", "id": "32c46b73", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Deriving the Ridge Regression Equations\n", "\n", @@ -714,9 +607,7 @@ { "cell_type": "markdown", "id": "3cceda87", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{X},\\boldsymbol{\\beta})=\\left\\{(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})^T(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\right\\}+\\lambda\\boldsymbol{\\beta}^T\\boldsymbol{\\beta},\n", @@ -726,9 +617,7 @@ { "cell_type": "markdown", "id": "10a16675", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and \n", "taking the derivatives with respect to $\\boldsymbol{\\beta}$ we obtain then\n", @@ -740,9 +629,7 @@ { "cell_type": "markdown", "id": "6dba20bf", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}_{\\mathrm{Ridge}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}+\\lambda\\boldsymbol{I}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y},\n", @@ -752,9 +639,7 @@ { "cell_type": "markdown", "id": "9325f924", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "with $\\boldsymbol{I}$ being a $p\\times p$ identity matrix with the constraint that" ] @@ -762,9 +647,7 @@ { "cell_type": "markdown", "id": "745d61e9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\sum_{i=0}^{p-1} \\beta_i^2 \\leq t,\n", @@ -774,9 +657,7 @@ { "cell_type": "markdown", "id": "3866af64", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "with $t$ a finite positive number." ] @@ -784,9 +665,7 @@ { "cell_type": "markdown", "id": "710de288", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Note on Scikit-Learn\n", "\n", @@ -796,9 +675,7 @@ { "cell_type": "markdown", "id": "5893fba0", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}_{\\mathrm{Ridge}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}+n\\lambda\\boldsymbol{I}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n", @@ -808,9 +685,7 @@ { "cell_type": "markdown", "id": "5ef8c59f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "In our codes where we compare our own codes with **Scikit-Learn**, we do thus not include the $1/n$ factor in the cost function." ] @@ -818,9 +693,7 @@ { "cell_type": "markdown", "id": "79b257ae", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Comparison with OLS\n", "When we compare this with the ordinary least squares result we have" @@ -829,9 +702,7 @@ { "cell_type": "markdown", "id": "46270996", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}_{\\mathrm{OLS}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y},\n", @@ -841,9 +712,7 @@ { "cell_type": "markdown", "id": "785aebe9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which can lead to singular matrices. However, with the SVD, we can always compute the inverse of the matrix $\\boldsymbol{X}^T\\boldsymbol{X}$.\n", "\n", @@ -857,9 +726,7 @@ { "cell_type": "markdown", "id": "443f3a3f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## SVD analysis\n", "\n", @@ -870,9 +737,7 @@ { "cell_type": "markdown", "id": "cfb8d4c0", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\tilde{\\boldsymbol{y}}_{\\mathrm{OLS}}=\\boldsymbol{X}\\boldsymbol{\\beta} =\\boldsymbol{U}\\boldsymbol{U}^T\\boldsymbol{y}.\n", @@ -882,9 +747,7 @@ { "cell_type": "markdown", "id": "2be1ca1f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "For Ridge regression this becomes" ] @@ -892,9 +755,7 @@ { "cell_type": "markdown", "id": "9ed0aec2", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\tilde{\\boldsymbol{y}}_{\\mathrm{Ridge}}=\\boldsymbol{X}\\boldsymbol{\\beta}_{\\mathrm{Ridge}} = \\boldsymbol{U\\Sigma V^T}\\left(\\boldsymbol{V}\\boldsymbol{\\Sigma}^2\\boldsymbol{V}^T+\\lambda\\boldsymbol{I} \\right)^{-1}(\\boldsymbol{U\\Sigma V^T})^T\\boldsymbol{y}=\\sum_{j=0}^{p-1}\\boldsymbol{u}_j\\boldsymbol{u}_j^T\\frac{\\sigma_j^2}{\\sigma_j^2+\\lambda}\\boldsymbol{y},\n", @@ -904,9 +765,7 @@ { "cell_type": "markdown", "id": "943b1cd5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "with the vectors $\\boldsymbol{u}_j$ being the columns of $\\boldsymbol{U}$ from the SVD of the matrix $\\boldsymbol{X}$." ] @@ -914,9 +773,7 @@ { "cell_type": "markdown", "id": "4a8a6bff", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Interpreting the Ridge results\n", "\n", @@ -926,9 +783,7 @@ { "cell_type": "markdown", "id": "b0371965", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\frac{\\sigma_j^2}{\\sigma_j^2+\\lambda} \\leq 1.\n", @@ -938,9 +793,7 @@ { "cell_type": "markdown", "id": "38e7b497", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Ridge regression finds the coordinates of $\\boldsymbol{y}$ with respect to the\n", "orthonormal basis $\\boldsymbol{U}$, it then shrinks the coordinates by\n", @@ -954,9 +807,7 @@ { "cell_type": "markdown", "id": "682f90db", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## More interpretations\n", "\n", @@ -966,9 +817,7 @@ { "cell_type": "markdown", "id": "b5fb7771", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}^T\\boldsymbol{X}=(\\boldsymbol{X}^T\\boldsymbol{X})^{-1} =\\boldsymbol{I}.\n", @@ -978,9 +827,7 @@ { "cell_type": "markdown", "id": "0c228525", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "In this case the standard OLS results in" ] @@ -988,9 +835,7 @@ { "cell_type": "markdown", "id": "08a400af", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{\\beta}^{\\mathrm{OLS}} = \\boldsymbol{X}^T\\boldsymbol{y}=\\sum_{i=0}^{n-1}\\boldsymbol{u}_i\\boldsymbol{u}_i^T\\boldsymbol{y},\n", @@ -1000,9 +845,7 @@ { "cell_type": "markdown", "id": "ec0808b1", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and" ] @@ -1010,9 +853,7 @@ { "cell_type": "markdown", "id": "92d3d92f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{\\beta}^{\\mathrm{Ridge}} = \\left(\\boldsymbol{I}+\\lambda\\boldsymbol{I}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}=\\left(1+\\lambda\\right)^{-1}\\boldsymbol{\\beta}^{\\mathrm{OLS}},\n", @@ -1022,9 +863,7 @@ { "cell_type": "markdown", "id": "07486e29", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "that is the Ridge estimator scales the OLS estimator by the inverse of a factor $1+\\lambda$, and\n", "the Ridge estimator converges to zero when the hyperparameter goes to\n", @@ -1039,9 +878,7 @@ { "cell_type": "markdown", "id": "9868d4f4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Deriving the Lasso Regression Equations\n", "\n", @@ -1051,9 +888,7 @@ { "cell_type": "markdown", "id": "1793a225", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{X},\\boldsymbol{\\beta})=\\frac{1}{n}\\left\\{(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})^T(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\right\\}+\\lambda\\vert\\vert\\boldsymbol{\\beta}\\vert\\vert_1,\n", @@ -1063,9 +898,7 @@ { "cell_type": "markdown", "id": "c7d3fa52", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Taking the derivative with respect to $\\boldsymbol{\\beta}$ and recalling that the derivative of the absolute value is (we drop the boldfaced vector symbol for simplicity)" ] @@ -1073,9 +906,7 @@ { "cell_type": "markdown", "id": "47894977", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\frac{d \\vert \\beta\\vert}{d \\beta}=\\mathrm{sgn}(\\beta)=\\left\\{\\begin{array}{cc} 1 & \\beta > 0 \\\\-1 & \\beta < 0, \\end{array}\\right.\n", @@ -1085,9 +916,7 @@ { "cell_type": "markdown", "id": "d1a11dc8", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "we have that the derivative of the cost function is" ] @@ -1095,9 +924,7 @@ { "cell_type": "markdown", "id": "a6cf82de", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\frac{\\partial C(\\boldsymbol{X},\\boldsymbol{\\beta})}{\\partial \\boldsymbol{\\beta}}=-\\frac{2}{n}\\boldsymbol{X}^T(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})+\\lambda sgn(\\boldsymbol{\\beta})=0,\n", @@ -1107,9 +934,7 @@ { "cell_type": "markdown", "id": "bbf7ccc7", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and reordering we have" ] @@ -1117,9 +942,7 @@ { "cell_type": "markdown", "id": "bc6dbb83", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}^T\\boldsymbol{X}\\boldsymbol{\\beta}+\\lambda sgn(\\boldsymbol{\\beta})=\\boldsymbol{X}^T\\boldsymbol{y}.\n", @@ -1129,9 +952,7 @@ { "cell_type": "markdown", "id": "ef5a8292", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "This equation does not lead to a nice analytical equation as in Ridge regression or ordinary least squares. We have absorbed the factor $2/n$ in a redefinition of the parameter $\\lambda$. We will solve this type of problems using libraries like **scikit-learn**." ] @@ -1139,9 +960,7 @@ { "cell_type": "markdown", "id": "d77bf73a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Simple example to illustrate Ordinary Least Squares, Ridge and Lasso Regression\n", "\n", @@ -1154,9 +973,7 @@ { "cell_type": "markdown", "id": "265ce93f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta})=\\sum_{i=0}^{p-1}(y_i-\\beta_i)^2,\n", @@ -1166,9 +983,7 @@ { "cell_type": "markdown", "id": "c19e720c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and minimizing we have that" ] @@ -1176,9 +991,7 @@ { "cell_type": "markdown", "id": "8cad2962", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\beta}_i^{\\mathrm{OLS}} = y_i.\n", @@ -1188,9 +1001,7 @@ { "cell_type": "markdown", "id": "bf6316a2", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Ridge Regression\n", "\n", @@ -1200,9 +1011,7 @@ { "cell_type": "markdown", "id": "78413ab6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta})=\\sum_{i=0}^{p-1}(y_i-\\beta_i)^2+\\lambda\\sum_{i=0}^{p-1}\\beta_i^2,\n", @@ -1212,9 +1021,7 @@ { "cell_type": "markdown", "id": "a77ca41b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and minimizing we have that" ] @@ -1222,9 +1029,7 @@ { "cell_type": "markdown", "id": "67b3f67e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\beta}_i^{\\mathrm{Ridge}} = \\frac{y_i}{1+\\lambda}.\n", @@ -1234,9 +1039,7 @@ { "cell_type": "markdown", "id": "d3571c8b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Lasso Regression\n", "\n", @@ -1246,9 +1049,7 @@ { "cell_type": "markdown", "id": "943f822e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta})=\\sum_{i=0}^{p-1}(y_i-\\beta_i)^2+\\lambda\\sum_{i=0}^{p-1}\\vert\\beta_i\\vert=\\sum_{i=0}^{p-1}(y_i-\\beta_i)^2+\\lambda\\sum_{i=0}^{p-1}\\sqrt{\\beta_i^2},\n", @@ -1258,9 +1059,7 @@ { "cell_type": "markdown", "id": "2ec0ba01", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and minimizing we have that" ] @@ -1268,9 +1067,7 @@ { "cell_type": "markdown", "id": "85bf4c9b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "-2\\sum_{i=0}^{p-1}(y_i-\\beta_i)+\\lambda \\sum_{i=0}^{p-1}\\frac{(\\beta_i)}{\\vert\\beta_i\\vert}=0,\n", @@ -1280,9 +1077,7 @@ { "cell_type": "markdown", "id": "c81d05f4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which leads to" ] @@ -1290,9 +1085,7 @@ { "cell_type": "markdown", "id": "ce36390d", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}_i^{\\mathrm{Lasso}} = \\left\\{\\begin{array}{ccc}y_i-\\frac{\\lambda}{2} &\\mathrm{if} & y_i> \\frac{\\lambda}{2}\\\\\n", @@ -1304,9 +1097,7 @@ { "cell_type": "markdown", "id": "a47dd697", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Plotting these results shows clearly that Lasso regression suppresses (sets to zero) values of $\\beta_i$ for specific values of $\\lambda$. Ridge regression reduces on the other hand the values of $\\beta_i$ as function of $\\lambda$." ] @@ -1314,9 +1105,7 @@ { "cell_type": "markdown", "id": "15fe5a0c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Yet another Example\n", "\n", @@ -1326,9 +1115,7 @@ { "cell_type": "markdown", "id": "62d7f0e9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{y}=\\begin{bmatrix}4 \\\\ 2 \\\\3\\end{bmatrix},\n", @@ -1338,9 +1125,7 @@ { "cell_type": "markdown", "id": "f2d4576f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and our inputs as a $3\\times 2$ design matrix" ] @@ -1348,9 +1133,7 @@ { "cell_type": "markdown", "id": "fb2bc0ff", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}=\\begin{bmatrix}2 & 0\\\\ 0 & 1 \\\\ 0 & 0\\end{bmatrix},\n", @@ -1360,9 +1143,7 @@ { "cell_type": "markdown", "id": "807ac009", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "meaning that we have two features and two unknown parameters $\\beta_0$ and $\\beta_1$ to be determined either by ordinary least squares, Ridge or Lasso regression." ] @@ -1370,9 +1151,7 @@ { "cell_type": "markdown", "id": "fe72889b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## The OLS case\n", "\n", @@ -1382,9 +1161,7 @@ { "cell_type": "markdown", "id": "48844e10", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}^{\\mathrm{OLS}}=\\left( \\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n", @@ -1394,9 +1171,7 @@ { "cell_type": "markdown", "id": "72a4bc4c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Inserting the above values we obtain that" ] @@ -1404,9 +1179,7 @@ { "cell_type": "markdown", "id": "cfe1db6f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}^{\\mathrm{OLS}}=\\begin{bmatrix}2 \\\\ 2\\end{bmatrix},\n", @@ -1416,9 +1189,7 @@ { "cell_type": "markdown", "id": "05b6e6f6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "The code which implements this simpler case is presented after the discussion of Ridge and Lasso." ] @@ -1426,9 +1197,7 @@ { "cell_type": "markdown", "id": "0e58f37a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## The Ridge case\n", "\n", @@ -1438,9 +1207,7 @@ { "cell_type": "markdown", "id": "78593de5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}}=\\left( \\boldsymbol{X}^T\\boldsymbol{X}+\\lambda\\boldsymbol{I}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n", @@ -1450,9 +1217,7 @@ { "cell_type": "markdown", "id": "26b08a56", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Inserting the above values we obtain that" ] @@ -1460,9 +1225,7 @@ { "cell_type": "markdown", "id": "3df661a4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}}=\\begin{bmatrix}\\frac{8}{4+\\lambda} \\\\ \\frac{2}{1+\\lambda}\\end{bmatrix},\n", @@ -1472,9 +1235,7 @@ { "cell_type": "markdown", "id": "68460e68", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "There is normally a constraint on the value of $\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_2$ via the parameter $\\lambda$.\n", "Let us for simplicity assume that $\\beta_0^2+\\beta_1^2=1$ as constraint. This will allow us to find an expression for the optimal values of $\\beta$ and $\\lambda$.\n", @@ -1485,9 +1246,7 @@ { "cell_type": "markdown", "id": "6c660112", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Writing the Cost Function\n", "\n", @@ -1497,9 +1256,7 @@ { "cell_type": "markdown", "id": "6ff54cf1", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}\\boldsymbol{\\beta}=\\begin{bmatrix} 2\\beta_0 \\\\ \\beta_1 \\\\0 \\end{bmatrix},\n", @@ -1509,9 +1266,7 @@ { "cell_type": "markdown", "id": "f3f88f94", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta})=(4-2\\beta_0)^2+(2-\\beta_1)^2+\\lambda(\\beta_0^2+\\beta_1^2),\n", @@ -1521,9 +1276,7 @@ { "cell_type": "markdown", "id": "6918fc7c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and taking the derivative with respect to $\\beta_0$ we get" ] @@ -1531,9 +1284,7 @@ { "cell_type": "markdown", "id": "1b2321a3", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\beta_0=\\frac{8}{4+\\lambda},\n", @@ -1543,9 +1294,7 @@ { "cell_type": "markdown", "id": "d5e21457", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and for $\\beta_1$ we obtain" ] @@ -1553,9 +1302,7 @@ { "cell_type": "markdown", "id": "e03c3617", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\beta_1=\\frac{2}{1+\\lambda},\n", @@ -1565,9 +1312,7 @@ { "cell_type": "markdown", "id": "6c4d9a4a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Using the constraint for $\\beta_0^2+\\beta_1^2=1$ we can constrain $\\lambda$ by solving" ] @@ -1575,9 +1320,7 @@ { "cell_type": "markdown", "id": "584bef2e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\left(\\frac{8}{4+\\lambda}\\right)^2+\\left(\\frac{2}{1+\\lambda}\\right)^2=1,\n", @@ -1587,9 +1330,7 @@ { "cell_type": "markdown", "id": "87273242", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which gives $\\lambda=4.571$ and $\\beta_0=0.933$ and $\\beta_1=0.359$." ] @@ -1597,9 +1338,7 @@ { "cell_type": "markdown", "id": "f7e8277f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Lasso case\n", "\n", @@ -1610,9 +1349,7 @@ { "cell_type": "markdown", "id": "e7167c28", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta})=(4-2\\beta_0)^2+(2-\\beta_1)^2+\\lambda(\\vert\\beta_0\\vert+\\vert\\beta_1\\vert),\n", @@ -1622,9 +1359,7 @@ { "cell_type": "markdown", "id": "6a3fa7b2", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\frac{\\partial C(\\boldsymbol{\\beta})}{\\partial \\beta_0}=-4(4-2\\beta_0)+\\lambda\\mathrm{sgn}(\\beta_0)=0,\n", @@ -1634,9 +1369,7 @@ { "cell_type": "markdown", "id": "92c33067", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and" ] @@ -1644,9 +1377,7 @@ { "cell_type": "markdown", "id": "1b7be5cc", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\frac{\\partial C(\\boldsymbol{\\beta})}{\\partial \\beta_1}=-2(2-\\beta_1)+\\lambda\\mathrm{sgn}(\\beta_1)=0.\n", @@ -1656,9 +1387,7 @@ { "cell_type": "markdown", "id": "9e2db51b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We have now four cases to solve besides the trivial cases $\\beta_0$ and/or $\\beta_1$ are zero, namely\n", "1. $\\beta_0 > 0$ and $\\beta_1 > 0$,\n", @@ -1673,9 +1402,7 @@ { "cell_type": "markdown", "id": "7ea5c023", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## The first Case\n", "\n", @@ -1685,9 +1412,7 @@ { "cell_type": "markdown", "id": "efff1845", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "-4(4-2\\beta_0)+\\lambda=0,\n", @@ -1697,9 +1422,7 @@ { "cell_type": "markdown", "id": "80496d8c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and" ] @@ -1707,9 +1430,7 @@ { "cell_type": "markdown", "id": "bb33463a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "-2(2-\\beta_1)+\\lambda=0.\n", @@ -1719,9 +1440,7 @@ { "cell_type": "markdown", "id": "f16b1457", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which yields" ] @@ -1729,9 +1448,7 @@ { "cell_type": "markdown", "id": "4bee5a43", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\beta_0=\\frac{16+\\lambda}{8},\n", @@ -1741,9 +1458,7 @@ { "cell_type": "markdown", "id": "17176990", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and" ] @@ -1751,9 +1466,7 @@ { "cell_type": "markdown", "id": "ec2c8f34", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\beta_1=\\frac{4+\\lambda}{2}.\n", @@ -1763,9 +1476,7 @@ { "cell_type": "markdown", "id": "5934144d", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Using the constraint on $\\beta_0$ and $\\beta_1$ we can then find the optimal value of $\\lambda$ for the different cases. We leave this as an exercise to you." ] @@ -1773,9 +1484,7 @@ { "cell_type": "markdown", "id": "fe097c51", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Simple code for solving the above problem\n", "\n", @@ -1788,10 +1497,7 @@ "cell_type": "code", "execution_count": 4, "id": "92cd1361", - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", @@ -1850,9 +1556,7 @@ { "cell_type": "markdown", "id": "9e4426f4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We see here that we reach a plateau. What is actually happening?" ] @@ -1860,9 +1564,7 @@ { "cell_type": "markdown", "id": "5f2b6935", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## With Lasso Regression" ] @@ -1871,10 +1573,7 @@ "cell_type": "code", "execution_count": 5, "id": "409bbeac", - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import os\n", @@ -1938,22 +1637,39 @@ { "cell_type": "markdown", "id": "67a24e77", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Another Example, now with a polynomial fit" ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 3, "id": "8da16a74", - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[2.0000000e+00 7.1942452e-14 5.0000000e+00]\n", + "Training MSE for OLS\n", + "9.647251180883923e-28\n", + "Test MSE OLS\n", + "1.2167390602128887e-27\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAG1CAYAAAAMU3WaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABzlElEQVR4nO3dd3gUVdvH8e+m91BCCCWh996bUlSkiYioqIgCyiMiIKIIqIDyqAiPAjbECqivgAVQERBEepMWQJr0DqGmENJ25/1jk4WFJBBIMtnk97muvbIzc8/uvcuSvXPOmXMshmEYiIiIiLgoN7MTEBEREbkdKmZERETEpamYEREREZemYkZERERcmooZERERcWkqZkRERMSlqZgRERERl6ZiRkRERFyaihkRERFxaSpmRERExKWZWsysWLGCzp07U7JkSSwWC3Pnzs0w9tlnn8VisTBp0qRcy09ERETyPlOLmUuXLlGnTh0+/vjjTOPmzp3L+vXrKVmyZC5lJiIiIq7Cw8wn79ChAx06dMg05vjx4wwYMIA//viDTp06Zfk5bDYbJ06cIDAwEIvFcqupioiISC4yDIPY2FhKliyJm1vmbS+mFjM3YrPZ6NmzJ0OHDqVGjRo3dU5iYiKJiYmO7ePHj1O9evWcSlFERERy0NGjRyldunSmMXm6mBk3bhweHh4MGjTops8ZO3Ysb7755nX7jx49SlBQUHamJyIiIjkkJiaG8PBwAgMDbxibZ4uZTZs28cEHH7B58+YsdQ+NGDGCIUOGOLbT3oygoCAVMyIiIi7mZmqAPHtp9sqVK4mKiiIiIgIPDw88PDw4fPgwL730EmXLls3wPG9vb0fhogJGREQk/8uzLTM9e/bknnvucdrXrl07evbsSe/evU3KSkRERPIaU4uZuLg49u3b59g+ePAgkZGRFClShIiICIoWLeoU7+npSVhYGFWqVMntVEVERCSPMrWY2bhxI23atHFsp411eeqpp5g2bVqu5mK1WklOTs7V55SCzdPTE3d3d7PTEBFxeaYWM61bt8YwjJuOP3ToULbnYBgGp06d4uLFi9n+2CI3UqhQIcLCwjQHkojIbcizY2ZyS1ohExoaip+fn75UJFcYhkF8fDxRUVEAlChRwuSMRERcV4EuZqxWq6OQuXZ8jkhO8/X1BSAqKorQ0FB1OYmI3KI8e2l2bkgbI+Pn52dyJlJQpX32NF5LROTWFehiJo26lsQs+uyJiNw+FTMiIiLi0lTMSI46dOgQFouFyMjIDGOWLVuGxWJxuSvKWrduzeDBg81OQ0SkwFMx44J69eqFxWKhX79+1x3r378/FouFXr16OfZFRUXx7LPPEhERgbe3N2FhYbRr1461a9c6YsqWLYvFYrnu9u6772aYR+vWrR1xXl5eVKhQgREjRjitWh4eHs7JkyepWbNm9rz429CrVy8eeOCBbHu82bNn89///jfbHk9ERG5Ngb6ayZWFh4czc+ZMJk6c6LgqJiEhgRkzZhAREeEU261bN5KTk5k+fTrly5fn9OnTLFmyhPPnzzvFjRkzhr59+zrtu9FqpX379mXMmDEkJSWxYcMGx1ITY8eOBcDd3Z2wsLDbeq25LTk5GU9PzxvGFSlSJBeyERHJuwzDYM/pWEICvAkJ8DYtD7XMuKj69esTERHB7NmzHftmz55NeHg49erVc+y7ePEiq1atYty4cbRp04YyZcrQuHFjRowYQadOnZweMzAwkLCwMKebv79/pnn4+fkRFhZGREQE3bp1o23btixatMhxPL1upvnz51O5cmV8fX1p06ZNupMhfvHFF4SHh+Pn50fXrl2ZMGEChQoVcor57bffaNCgAT4+PpQvX54333yTlJSUdPN84403mD59Or/88oujNWnZsmWO/H744Qdat26Nj48P3333HefOneOxxx6jdOnS+Pn5UatWLWbMmOH0mNd2M5UtW5Z33nmHPn36EBgYSEREBJ9//nmm75+IiKtJttpYve8sb/y6gzvHL6X9pJXM3XLc1JzUMnMNwzC4nGzN9ef19XTP8pUtvXv3ZurUqfTo0QOAr7/+mj59+rBs2TJHTEBAAAEBAcydO5emTZvi7Z1zlfPWrVtZvXp1pquaHz16lAcffJB+/frx3HPPsXHjRl566SWnmNWrV9OvXz/GjRvH/fffz59//snIkSOdYv744w+eeOIJPvzwQ+68807279/Pf/7zHwBGjx593fO+/PLL7Nq1i5iYGKZOnQrYW1ZOnDgBwLBhw3j//feZOnUq3t7eJCQk0KBBA4YNG0ZQUBC///47PXv2pHz58jRp0iTD1/f+++/z3//+l1dffZWffvqJ5557jpYtW1K1atWbeg9FRPKiFKuNlfvOMnfLcf7aHUVswpU/HL093LgQn2RidipmrnM52Ur1UX/k+vPuHNMOP6+s/XP07NmTESNGOFoXVq9ezcyZM52KGQ8PD6ZNm0bfvn2ZMmUK9evXp1WrVjz66KPUrl3b6fGGDRvG66+/7rRv3rx5tG7dOsMcJk+ezJdffklycjJJSUm4ubnxySefZBj/6aefUr58eSZOnIjFYqFKlSps376dcePGOWI++ugjOnTowMsvvwxA5cqVWbNmDfPmzXPEvP322wwfPpynnnoKgPLly/Pf//6XV155Jd1iJiAgAF9fXxITE9Pt9ho8eDAPPvig07605wcYOHAgCxcu5Mcff8y0mOnYsSP9+/cH7O/nxIkTWbZsmYoZEXE5hmGw40QMc7Yc55fIE5yNuzIeMiTAi7uqhnJPteLcUSkky99f2U3FjAsLCQmhU6dOTJ8+HcMw6NSpEyEhIdfFdevWjU6dOrFy5UrWrl3LwoULGT9+PF9++aXTQOGhQ4c6bQOUKlUq0xx69OjBa6+9RkxMDOPGjSMoKIhu3bplGL9r1y6aNm3q1ArVrFkzp5g9e/bQtWtXp32NGzd2KmY2bdrEhg0bePvttx37rFYrCQkJxMfHZ3kixIYNGzptW61W3n33XWbNmsXx48dJTEwkMTHxht1uVxeIFouFsLAwx5IFIiKuICnFxvztJ/l69UG2HYt27C/i70Xn2iW4v25J6oYXxt0t78yTpWLmGr6e7uwc086U570Vffr0YcCAAQCZtoj4+PjQtm1b2rZty6hRo3jmmWcYPXq0U/ESEhJCxYoVs/T8wcHBjnO+++47atSowVdffcXTTz+dbvzNLCxqGMZ1XW7Xnmez2XjzzTeva00B+2vNqmuLlPfff5+JEycyadIkatWqhb+/P4MHDyYpKfOm1GsHDlssFmw2W5bzERHJbWfjEvl+/RG+W3eYqFh7K4yXhxttqxWna71StKpSDE/3vDnUVsXMNSwWi+nNZVnRvn17xxdsu3Y3X4RVr16duXPnZmsunp6evPrqq4wYMYLHHnss3daR9J533bp1TttVq1bl77//dtq3ceNGp+369euzZ8+eLBVfXl5eWK03Nx5q5cqVdOnShSeeeAKwF0979+6lWrVqN/18IiKu4PylJKYs38/0NYdITLH/8VU8yJueTcvwWOMIipp4ldLNcp1vbUmXu7s7u3btcty/1rlz53j44Yfp06cPtWvXJjAwkI0bNzJ+/Hi6dOniFBsbG8upU6ec9vn5+REUFHTT+Tz++OO8+uqrTJ482WnMSZp+/frx/vvvM2TIEJ599lk2bdrEtGnTnGIGDhxIy5YtmTBhAp07d+avv/5iwYIFTq01o0aN4r777iM8PJyHH34YNzc3tm3bxvbt23nrrbfSza1s2bL88ccf7Nmzh6JFixIcHJzh66hYsSI///wza9asoXDhwkyYMIFTp06pmBGRfCP6cjJfrTzAV6sOcinJ/odendLB9LmjHB1qlsDLI2+2wqTHdTKVDAUFBWVYcAQEBNCkSRMmTpxIy5YtqVmzJiNHjqRv3758/PHHTrGjRo2iRIkSTrdXXnklS7l4eXkxYMAAxo8fT1xc3HXHIyIi+Pnnn/ntt9+oU6cOU6ZM4Z133nGKadGiBVOmTGHChAnUqVOHhQsX8uKLLzp1H7Vr14558+axePFiGjVqRNOmTZkwYQJlypTJMLe+fftSpUoVGjZsSLFixVi9enWGsSNHjqR+/fq0a9eO1q1bExYWlq0T7omImCXFauPLlQdoOX4pH/61j0tJVmqVCmZa70bMfb4FXeqWcqlCBsBi3MwgBhcWExNDcHAw0dHR133hJyQkcPDgQcqVK3dL4ywk9/Tt25fdu3ezcuVKs1PJVvoMikhu2nzkAq/N+YddJ2MAqBQawEv3VqZdjbA8t/BtZt/f11I3k+RJ7733Hm3btsXf358FCxYwffp0Jk+ebHZaIiIuKTo+mXF/7GbG30cwDCjk58nw9lV5uGF4nroq6VapmJE86e+//2b8+PHExsZSvnx5PvzwQ5555hmz0xIRcTl/7T7NKz9t42yc/WKRhxuUZkTHahTx9zI5s+yjYkbypB9++MHsFEREXFpiipV3F+xm6upDgL1L6a0HatKkfFFzE8sBKmZERETymQNn4hg4Yws7TtjHxvRpUY5hHarg7XFrc5rldSpmRERE8pFfIo8zYvZ24pOsFPbz5L2H63B3teJmp5WjVMyIiIjkA4Zh8MGSvUz6cy8ATcsXYVL3eoQF5/8rJVXMiIiIuLikFBvDZ29j9ubjAPRrVYGh7arkiyuVboaKGRERERcWHZ9Mv+82sfbAOdzdLLz1QE0eaxxhdlq5SsWMiIiIizpx8TJPfv03+6LiCPD24JMe9WlVuZjZaeU615qvWFzOoUOHsFgsREZGZhizbNkyLBYLFy9ezLW8RERc3anoBB77Yh37ouIIC/Lhh2ebFchCBlTMuKRevXphsVjo16/fdcf69++PxWKhV69ejn1RUVE8++yzRERE4O3tTVhYGO3atWPt2rWOmLJly2KxWK67vfvuuxnm0bp1a0ecl5cXFSpUYMSIESQmJjpiwsPDOXnyJDVr1syeF3+LMnp9abfWrVvf8mO3bt2awYMHZ1uuIiI3EhWbwONfrOPwuXjCi/gyu39zqpe8+UWB8xt1M7mo8PBwZs6cycSJE/H19QXs6/zMmDGDiAjnvtJu3bqRnJzM9OnTKV++PKdPn2bJkiWcP3/eKW7MmDH07dvXaV9gYGCmefTt25cxY8aQlJTEhg0b6N27NwBjx44F7Ct5h4WF3dZrzQ4bNmzAarWvCrtmzRq6devGnj17HOt9eHnln5kwRSR/OxuXyONfrOfA2UuUKuTLjL5NKVnI1+y0TKWWGRdVv359IiIimD17tmPf7NmzCQ8Pp169eo59Fy9eZNWqVYwbN442bdpQpkwZGjduzIgRI+jUqZPTYwYGBhIWFuZ08/f3zzQPPz8/wsLCiIiIoFu3brRt25ZFixY5jqfXzTR//nwqV66Mr68vbdq04dChQ9c97hdffEF4eDh+fn507dqVCRMmUKhQIaeY3377jQYNGuDj40P58uV58803SUlJSTfPYsWKOV5TkSJFAAgNDXXs2717Ny1btsTX15fw8HAGDRrEpUuXHOdPnjyZSpUq4ePjQ/HixXnooYcAeyvZ8uXL+eCDDxytPOm9HhGR7HD+UhJPfLne0bU0o29TShf2Mzst06mYycilSxnfEhJuPvby5RvH3qLevXszdepUx/bXX39Nnz59nGICAgIICAhg7ty5Tt0/OWHr1q2sXr0aT0/PDGOOHj3Kgw8+SMeOHYmMjOSZZ55h+PDhTjGrV6+mX79+vPDCC0RGRtK2bVvefvttp5g//viDJ554gkGDBrFz504+++wzpk2bdl3czdi+fTvt2rXjwQcfZNu2bcyaNYtVq1YxYMAAADZu3MigQYMYM2YMe/bsYeHChbRs2RKADz74gGbNmtG3b19OnjzJyZMnCQ8Pz3IOIiI3cikxhSe/Xs/uU7GEBnoz4z9NiSiqQgYAI5+Ljo42ACM6Ovq6Y5cvXzZ27txpXL58+foTIeNbx47OsX5+Gce2auUcGxJyfUwWPfXUU0aXLl2MM2fOGN7e3sbBgweNQ4cOGT4+PsaZM2eMLl26GE899ZQj/qeffjIKFy5s+Pj4GM2bNzdGjBhhbN261ekxy5QpY3h5eRn+/v5Ot6VLl2aYR6tWrQxPT0/D39/f8PLyMgDDzc3N+OmnnxwxBw8eNABjy5YthmEYxogRI4xq1aoZNpvNETNs2DADMC5cuGAYhmF0797d6NSpk9Nz9ejRwwgODnZs33nnncY777zjFPPtt98aJUqUuOH7t3TpUqfn69mzp/Gf//zHKWblypWGm5ubcfnyZePnn382goKCjJiYmAzfhxdeeOGGz5ueTD+DIiKprFab0Xf6BqPMsHlG/TGLjL2nY81OKcdl9v19LY2ZcWEhISF06tSJ6dOnYxgGnTp1IiQk5Lq4bt260alTJ1auXMnatWtZuHAh48eP58svv3QaKDx06FCnbYBSpUplmkOPHj147bXXiImJYdy4cQQFBdGtW7cM43ft2kXTpk2xWK5M5NSsWTOnmD179tC1a1enfY0bN2bevHmO7U2bNrFhwwanlhir1UpCQgLx8fH4+d38XyubNm1i3759/N///Z9jn2EY2Gw2Dh48SNu2bSlTpgzly5enffv2tG/fnq5du2bpOUREbsf7i/ewaOdpvNzd+PzJhlQMDTA7pTxFxUxG4uIyPuZ+zUJdUVEZx7pd05OXzeMp+vTp4+gO+eSTTzKM8/HxoW3btrRt25ZRo0bxzDPPMHr0aKfiJSQkhIoVK2bp+YODgx3nfPfdd9SoUYOvvvqKp59+Ot14wzBu+JiGYTgVO+mdZ7PZePPNN3nwwQevO9/HJ2tTd9tsNp599lkGDRp03bGIiAi8vLzYvHkzy5YtY9GiRYwaNYo33niDDRs2XDeOR0Qku83ZcoxPlu4H4N1utWhQprDJGeU9KmYycoOBr7kSexPat29PUlISAO3atbvp86pXr87cuXOzNRdPT09effVVRowYwWOPPZZuy0V6z7tu3Tqn7apVq/L333877du4caPTdv369dmzZ0+Wi6/01K9fnx07dmT6WB4eHtxzzz3cc889jB49mkKFCvHXX3/x4IMP4uXl5bhSSkQkO206fIFhP28H4LnWFXiwfmmTM8qbNADYxbm7u7Nr1y527dqF+7UtRsC5c+e46667+O6779i2bRsHDx7kxx9/ZPz48XTp0sUpNjY2llOnTjndYmJispTP448/jsViYfLkyeke79evH/v372fIkCHs2bOH77//nmnTpjnFDBw4kPnz5zNhwgT27t3LZ599xoIFC5xaa0aNGsU333zDG2+8wY4dO9i1axezZs3i9ddfz1K+AMOGDWPt2rU8//zzREZGsnfvXn799VcGDhwIwLx58/jwww+JjIzk8OHDfPPNN9hsNqpUqQLY57BZv349hw4d4uzZs9hstiznICJyreMXL/PstxtJSrFxb/XiDL23itkpXXH+PMyaBb17w9KlZmejYiY/CAoKcsyXcq2AgACaNGnCxIkTadmyJTVr1mTkyJH07duXjz/+2Cl21KhRlChRwun2yiuvZCkXLy8vBgwYwPjx44lLp6suIiKCn3/+md9++406deowZcoU3nnnHaeYFi1aMGXKFCZMmECdOnVYuHAhL774olP3Ubt27Zg3bx6LFy+mUaNGNG3alAkTJlCmTJks5QtQu3Ztli9fzt69e7nzzjupV68eI0eOpESJEgAUKlSI2bNnc9ddd1GtWjWmTJnCjBkzqFGjBgAvv/wy7u7uVK9enWLFinHkyJEs5yAicrUUq41BM7ZwNi6JaiWCmNi9Lm55YdHIr76C5s2hWDF49FGYNg3mzDE7KyzGzQxicGExMTEEBwcTHR193Rd+QkICBw8epFy5clkeZyG5q2/fvuzevZuVK1eanUq20mdQRNLzwZ97mfjnvwR6ezD/hTsJL2LSBQcHD0KJEpD2+2nUKPjvf+33a9SA9u2hWze45kKO7JDZ9/e1NGZG8qT33nuPtm3b4u/vz4IFC5g+fXqGXVciIvnJpsMX+PCvvQD894GauV/IJCbC7NkwZQqsWGG/n3aF6ZNPQqlS0KEDROSdlblN7WZasWIFnTt3pmTJklgsFqeBocnJyQwbNoxatWrh7+9PyZIlefLJJzlx4oR5CUuu+fvvv2nbti21atViypQpfPjhhzzzzDNmpyUikqNiE5IZPGsLVpvBA3VL8kC9zKfHyFb79sErr0Dp0vD44/ZCxs0Nduy4ElOxIjz7bJ4qZMDklplLly5Rp04devfufd3cJPHx8WzevJmRI0dSp04dLly4wODBg7n//vuvu7JF8p8ffvjB7BRERHLd6F92cPT8ZUoX9mXMA7m0QG9sLDzzDPz4o30qV7AXNH37wtNP21ti8jhTi5kOHTrQoUOHdI8FBwezePFip30fffQRjRs35siRI9ctpigiIuLKfok8zuwtx3GzwKTudQnyyXhpmGwVEGBvlTEM+xiY556Djh3Bw3VGorhOpkB0dDQWi0UTlYmISL4SFZPA63P/AWDAXZVoWLZIzj3Znj3w7rvw8cf2uc8sFpg8GXx9oXbtnHveHOQyxUxCQgLDhw/n8ccfz3RUc2JiotOCilmdJ0VERCS3/ff3XcQmpFCndDCD7rr9yUDTdfEijBkDH30EKSlQpgy88Yb9WJMmOfOcucQl5plJTk7m0UcfxWaz3fCKlrFjxxIcHOy4aQVjERHJy1buPcNvW0/gZoG3u9bCwz2bv5ptNvjsM6hUCSZOtBcynTrBY49l7/OYKM8XM8nJyTzyyCMcPHiQxYsX3/Ba8xEjRhAdHe24HT16NJcyFRERyZqEZCsjU7uXnmxWlpqlgrP3CQ4ehDZtoF8/OHsWqlWDhQth3jyokodmFL5NebqbKa2Q2bt3L0uXLqVo0aI3PMfb2xtvb+9cyE5EROT2TFm+n0Pn4gkN9Oaleytn/xMMHWq/xNrfH955xz641zOXBhbnIlOLmbi4OPbt2+fYPnjwIJGRkRQpUoSSJUvy0EMPsXnzZubNm4fVauXUqVMAFClSBC8vL7PSFuGNN95g7ty5REZGmp2KiLioQ2cvMXmZfTXskfdVJzAnrl764ANISoJJk6B8+ex//DzC1G6mjRs3Uq9ePerVqwfAkCFDqFevHqNGjeLYsWP8+uuvHDt2jLp16zqtF7RmzRoz0zZdr169sFgs9OvX77pj/fv3x2Kx0KtXL8e+qKgonn32WSIiIvD29iYsLIx27dqxdu1aR0zZsmWxWCzX3d59990M82jdujWDBw/OzpeWY6ZNm5atV8G9/PLLLFmyJNseT0QKFsMwGPnLPySl2LizUgj31S6RPQ+8fDm8+eaV7VKl4Ndf83UhAya3zLRu3ZrMlobK58tG3Zbw8HBmzpzJxIkT8fX1BexXfM2YMeO6OXi6detGcnIy06dPp3z58pw+fZolS5Zw/vx5p7gxY8bQt29fp32BgYE5+0LymKSkpJtq9QsICCAgICAXMhKR/Gj+9lOs3HsWLw83xnSpicVym4tIGob98urBg+0DfGvXvrIEQQGQ5wcAS/rq169PREQEs2fPduybPXs24eHhjpYugIsXL7Jq1SrGjRtHmzZtKFOmDI0bN2bEiBF06tTJ6TEDAwMJCwtzuvn7+99yjsOGDaNy5cr4+flRvnx5Ro4cSXJysuP41q1badOmDYGBgQQFBdGgQQPH7M6HDx+mc+fOFC5cGH9/f2rUqMH8+fMd5y5fvpzGjRvj7e1NiRIlGD58OCkpKenmsWzZMnr37u2Yp8hisfBG6uWIZcuW5a233qJXr14EBwc7irkb5f7GG29Qt25dx3avXr144IEHeO+99yhRogRFixbl+eefdzpHRAQgKcXG2AW7AHiuVQXKhdz671nAvpbSM8/AgAH2Quaxx6Bdu2zI1HXk6QHApjAMSI7P/ef19LNPXJQFvXv3ZurUqfTo0QOAr7/+mj59+rBs2TJHTFoLwty5c2natGmuDo4ODAxk2rRplCxZku3bt9O3b18CAwN55ZVXAOjRowf16tXj008/xd3dncjISDxTB6Y9//zzJCUlsWLFCvz9/dm5c6ejJeT48eN07NiRXr168c0337B792769u2Lj4+Po0i5WvPmzZk0aRKjRo1iz549AE6tKv/73/8YOXIkr7/++k3nnp6lS5dSokQJli5dyr59++jevTt169a9rrVLRAq2GX8f4diFy4QGetOvVYXbe7DTp+GBB2DdOvs6SuPGwUsvZfn7xNWpmLlWcjy8UzL3n/fVE+CVteq8Z8+ejBgxgkOHDmGxWFi9ejUzZ850KmY8PDyYNm0affv2ZcqUKdSvX59WrVrx6KOPUvuamR6HDRvm9IUOMG/ePFq3bn1LL+nqxypbtiwvvfQSs2bNchQER44cYejQoVStWhWASpUqOeKPHDlCt27dqFWrFgDlr+rvnTx5MuHh4Xz88cdYLBaqVq3KiRMnGDZsGKNGjcLNzbnB0cvLi+DgYCwWC2FhYdfledddd/Hyyy9nKff0FC5cmI8//hh3d3eqVq1Kp06dWLJkiYoZEXG4lJjCR6krYg+6uxK+Xu63/mBHj8Ldd8PevVC4MMycCffem02ZuhYVMy4sJCSETp06MX36dAzDoFOnToSEhFwX161bNzp16sTKlStZu3YtCxcuZPz48Xz55ZdOA4WHDh3qtA1Q6jYWGPvpp5+YNGkS+/btIy4ujpSUFKd5goYMGcIzzzzDt99+yz333MPDDz9MhQr2v1IGDRrEc889x6JFi7jnnnvo1q2bo/jatWsXzZo1c+pjbtGiBXFxcRw7dizL63Y1bNgwy7mnp0aNGri7X/nFVKJECbZv356lXEQkf/t61UHOxiVRpqgf3Rvd5qSuq1fbC5mICPjzT/ukeAWUiplrefrZW0nMeN5b0KdPHwYMGADAJ598kmGcj48Pbdu2pW3btowaNYpnnnmG0aNHOxUvISEhVKyYPdNor1u3jkcffZQ333yTdu3aERwczMyZM3n//fcdMW+88QaPP/44v//+OwsWLGD06NHMnDmTrl278swzz9CuXTt+//13Fi1axNixY3n//fcZOHAghmFcN1gubbD4rQyiu3Zc0M3knh7Pa+ZusFgs2Gy2LOcjIvnT+UtJfL7iAAAv3VsFz9ud6ffRR+2XXbdubS9oCjAVM9eyWLLc3WOm9u3bk5SUBEC7LAz4ql69OnPnzs2hrGD16tWUKVOG1157zbHv8OHD18VVrlyZypUr8+KLL/LYY48xdepUuqaOwA8PD6dfv37069ePESNG8MUXXzBw4ECqV6/Ozz//7FTUrFmzhsDAwAxbkry8vLBardmau4hIVny6bB+xiSlULxHEfbVu8VLsyEgoUQKKF7dvP/lktuXnylTMuDh3d3d27drluH+tc+fO8fDDD9OnTx9q165NYGAgGzduZPz48XTp0sUpNjY21jExYRo/P79Mu1fOnDlz3cRxYWFhVKxYkSNHjjBz5kwaNWrE77//zpw5cxwxly9fZujQoTz00EOUK1eOY8eOsWHDBrp16wbA4MGD6dChA5UrV+bChQv89ddfVKtWDbDPpTNp0iQGDhzIgAED2LNnD6NHj2bIkCHXjZdJU7ZsWeLi4liyZAl16tTBz88PP7/0W8NulLuISFaduHiZ6WvtfxS90r4Kbm63MEB3xw770gTh4bBsGRTJwZW1XYwuzc4HgoKCMiw4AgICaNKkCRMnTqRly5bUrFmTkSNH0rdvXz7++GOn2FGjRjlNTliiRIlMB7wCfP/9946JD9NuU6ZMoUuXLrz44osMGDCAunXrsmbNGkaOHOk4z93dnXPnzvHkk09SuXJlHnnkETp06MCbqZM9Wa1Wnn/+eapVq0b79u2pUqWKY5HRUqVKMX/+fP7++2/q1KlDv379ePrpp68bvHy15s2b069fP7p3706xYsUYP358hrE3yl1EJKs++HMvSSk2mpQrQqvKxbL+AMeOQfv29pWv/f3BQ20RV7MY+XxmupiYGIKDg4mOjr7uCz8hIYGDBw9Srlw5fHx8TMpQCjJ9BkXyv4NnL3H3+8uwGTC7f3PqRxTO2gNcuAB33mlvmalaFVatgptYq9DVZfb9fS21zIiIiOSgz1fsx2bA3VVDs17IJCTY55HZsQNKlrSveF0ACpmsUjEjIiKSQ6JiEvh503EAnmudxQnyrFZ44gn7qtdBQbBgAZQpkwNZuj4VMyIiIjnk69WHSLLaaFCmMA3LZnHA7uHDsHYteHnB3Ln29ZYkXRpBJCIikgNiEpL5v3X2K5huadmC8uVhyxbYutV+FZNkSC0zIiIiOWDG+iPEJqZQKTSAu6uG3vyJV1+XExoKbdtmf3L5jIoZERGRbJaYYuWrVQcB+E/L8jc/r0xCArRsCd99l4PZ5T8qZkRERLLZ3C3HiYpNJCzIhy51s7DG3Qsv2C+9fvFF+5wyclNUzIiIiGQjm83gs9Q1mJ6+oxxeHjf5VfvNN/D55/Zldb7/HgoVyrkk8xkVMyIiItlo0c7THDhziSAfDx5rcpMLQO7bB/362e+/8YbGyWSRihkX1KtXLywWC/3SPvhX6d+/PxaLxWk17KioKJ599lkiIiLw9vYmLCyMdu3asXbtWkdM2bJlsVgs193efffdDPNo3bo1gwcPzs6XliPS3q/MbrfqjTfeoG7dutmXrIi4vC9X2ltlejYrQ4D3TVw0bLPB00/D5ctw112QydIskj5dmu2iwsPDmTlzJhMnTsTX1xewT40/Y8YMIq5ZCr5bt24kJyczffp0ypcvz+nTp1myZAnnz593ihszZgx9+/Z12hcYGJizLyQXfPDBB05FWYkSJZg6dSrt27c3MSsRyY92nYxh4+ELeLhZeKpZ2Zs7afJk+8R4/v7w1VeQwYK5kjG9Yy6qfv36REREMHv2bMe+2bNnEx4eTr169Rz7Ll68yKpVqxg3bhxt2rShTJkyNG7cmBEjRtCpUyenxwwMDCQsLMzp5u/vf8s5Dhs2jMqVK+Pn50f58uUZOXIkycnJjuNbt26lTZs2BAYGEhQURIMGDdi4cSMAhw8fpnPnzhQuXBh/f39q1KjB/PnzHecuX76cxo0b4+3tTYkSJRg+fDgpKSnp5hEcHOz0mgAKFSrk2LZarXTv3p3ChQtTtGhRunTpwqFDhxznL1u2jMaNG+Pv70+hQoVo0aIFhw8fZtq0abz55pts3brV0cIzbdq0W36/RMT1/d96+7wy99YoTmjQTa63duqUfZzMuHFQtmzOJZePqZhxYb1792bq1KmO7a+//po+ffo4xQQEBBAQEMDcuXNJTEzM1fwCAwOZNm0aO3fu5IMPPuCLL75g4sSJjuM9evSgdOnSbNiwgU2bNjF8+HA8PT0BeP7550lMTGTFihVs376dcePGERAQAMDx48fp2LEjjRo1YuvWrXz66ad89dVXvPXWW1nOMT4+njZt2hAQEMCKFStYtWoVAQEBtG/fnqSkJFJSUnjggQdo1aoV27ZtY+3atfznP//BYrHQvXt3XnrpJWrUqMHJkyc5efIk3bt3z543T0RcTlxiCnM225cueKJJFpYdeOst2LwZnnsuhzIrAIx8Ljo62gCM6Ojo645dvnzZ2Llzp3H58uXrT4yLy/h2bXxmsfHxN47Noqeeesro0qWLcebMGcPb29s4ePCgcejQIcPHx8c4c+aM0aVLF+Opp55yxP/0009G4cKFDR8fH6N58+bGiBEjjK1btzo9ZpkyZQwvLy/D39/f6bZ06dIM82jVqpXxwgsv3HTe48ePNxo0aODYDgwMNKZNm5ZubK1atYw33ngj3WOvvvqqUaVKFcNmszn2ffLJJ0ZAQIBhtVpvmAdgzJkzxzAMw/jqq6+ue6zExETD19fX+OOPP4xz584ZgLFs2bJ0H2v06NFGnTp1bvicGcn0MygiLuW7dYeMMsPmGW3+t9Tpd4rcmsy+v6+llpmMBARkfOvWzTk2NDTj2A4dnGPLlr0+5haFhITQqVMnpk+fztSpU+nUqRMhISHXxXXr1o0TJ07w66+/0q5dO5YtW0b9+vWv6xIZOnQokZGRTrcmTZrccn4//fQTd9xxB2FhYQQEBDBy5EiOHDniOD5kyBCeeeYZ7rnnHt59913279/vODZo0CDeeustWrRowejRo9m2bZvj2K5du2jWrJnTwN0WLVoQFxfHsWPHspTjpk2b2LdvH4GBgY5WrCJFipCQkMD+/fspUqQIvXr1ol27dnTu3JkPPviAkydP3vJ7IiL5k2EYfLfO/vvt8SYRN76w4MgR6NgR/v03F7LL/1TMuLg+ffowbdo0pk+ffl0X09V8fHxo27Yto0aNYs2aNfTq1YvRo0c7xYSEhFCxYkWnW9rg4qxat24djz76KB06dGDevHls2bKF1157jaSkJEfMG2+8wY4dO+jUqRN//fUX1atXZ86cOQA888wzHDhwgJ49e7J9+3YaNmzIRx99BNh/aVz7i8JInf47q1cm2Ww2GjRocF0R9++///L4448DMHXqVNauXUvz5s2ZNWsWlStXZt26dbf0vohI/rT5yEV2nYzB28ONhxqUvvEJAwfaV8FW11K20NVMGYmLy/iYu7vzdlRUxrHXjkq/amBpdkgb2wHQrl27mz6vevXqzJ07N1tzudrq1aspU6YMr732mmPf4cOHr4urXLkylStX5sUXX+Sxxx5j6tSpdO3aFbBfsdWvXz/69evHiBEj+OKLLxg4cCDVq1fn559/dipq1qxZQ2BgIKVKZWGmTewDqWfNmkVoaChBQUEZxtWrV4969eoxYsQImjVrxvfff0/Tpk3x8vLCarVm6TlFJP9JW1Cyc52SFPLzyjz4zz/h11/t3yWpf6TJ7VExk5GsXMWTU7E3wd3dnV27djnuX+vcuXM8/PDD9OnTh9q1axMYGMjGjRsZP348Xbp0cYqNjY3l1KlTTvv8/Pwy/ZI/c+YMkZGRTvvCwsKoWLEiR44cYebMmTRq1Ijff//d0eoCcPnyZYYOHcpDDz1EuXLlOHbsGBs2bKBbahfe4MGD6dChA5UrV+bChQv89ddfVKtWDbDPpTNp0iQGDhzIgAED2LNnD6NHj2bIkCG4ZfGSxh49evC///2PLl26MGbMGEqXLs2RI0eYPXs2Q4cOJTk5mc8//5z777+fkiVLsmfPHv7991+efPJJwD4/z8GDB4mMjKR06dIEBgbi7e2dpRxExLVduJTEvO327ucnmt5g4G9Kin2pAoD+/aF69RzOroDI4fE7prvlAcB5WNoA4IxcPQA4ISHBGD58uFG/fn0jODjY8PPzM6pUqWK8/vrrRvxVg5PLlCljANfdnn322Qyfp1WrVumeM3r0aMMwDGPo0KFG0aJFjYCAAKN79+7GxIkTjeDgYMMw7INsH330USM8PNzw8vIySpYsaQwYMMDxbzFgwACjQoUKhre3t1GsWDGjZ8+extmzZx3PvWzZMqNRo0aGl5eXERYWZgwbNsxITk6+qfePqwYAG4ZhnDx50njyySeNkJAQw9vb2yhfvrzRt29fIzo62jh16pTxwAMPGCVKlDC8vLyMMmXKGKNGjXIMNE5ISDC6detmFCpUyACMqVOn3lQOaVz1MygiV3y+fL9RZtg8o+MHK2488PfTTw0DDKNwYcO46neaXC8rA4AthnH1WuP5T0xMDMHBwURHR1/XwpCQkMDBgwcpV64cPj43OR+ASDbSZ1DEtdlsBne9v4xD5+IZ+2AtHmucyfIFFy9CpUpw9ix8+KF93IxkKLPv72tpALCIiMgtWnvgHIfOxRPo7cH9dUpmHvzxx/ZCpmrVK+swSbbQmBkREZFb9PNm+3QQneuWxP9G6zC98gr4+UHNmpA6QahkDxUzIiIityA+KYWF/9gvmuhW/yaupPTygiFDcjirgkndTCIiIrfgjx2niE+yUqaoH/UjCmccePo0aAqHHKViRkRE5BbMTl2HqWu9UplP2Pn441CjBmiyzRyjbiauzB4rktv02RNxTVExCazedxawFzMZWrkS/vrLPkam5A0GCMstK9AtM2krNMfHx5uciRRUaZ89Tw0GFHEpv0SewGZAgzKFKVM0k8lQ33zT/rNPH4jI5LJtuS0FumXG3d2dQoUKEZW6HIGfn1+W1/YRuRWGYRAfH09UVBSFChVKd/ZmEcm7Zm+50sWUoZUrYckSe6vMiBG5lFnBVKCLGbBPvQ84ChqR3FSoUCHHZ1BEXMOukzHsOhmDl7sb99UukXHg1a0yZW6wzIHclgJfzFgsFkqUKEFoaCjJyclmpyMFiKenp1pkRFzQnNRWmbuqhma8qOSqVWqVyUWmFjMrVqzgf//7H5s2beLkyZPMmTOHBx54wHHcMAzefPNNPv/8cy5cuECTJk345JNPqFGjRrbn4u7uri8WERHJlNVm8EtkahdTZnPLpC2s27u3WmVygakDgC9dukSdOnX4+OOP0z0+fvx4JkyYwMcff8yGDRsICwujbdu2xMbG5nKmIiIisGb/WU7HJFLIz5M2VUIzDnzvPfjzT3jttdxLrgAztWWmQ4cOdOjQId1jhmEwadIkXnvtNR588EEApk+fTvHixfn+++959tlnczNVERER5qTOLXNf7RJ4eWTSHmCxwN1351JWkmcvzT548CCnTp3i3nvvdezz9vamVatWrFmzxsTMRESkIEpItvLHDvvyBV3rlU4/6PRpiI7OxawE8nAxc+qU/QNTvHhxp/3Fixd3HEtPYmIiMTExTjcREZHbteLfM1xKslIy2If6EYXSD3r9dQgPh+nTczW3gi7PFjNprp33xTCMTOeCGTt2LMHBwY5beHh4TqcoIiIFwILURSU71CqR/vfQmTPw7bcQGwsVKuRydgVbni1m0ubeuLYVJioq6rrWmquNGDGC6Ohox+3o0aM5mqeIiOR/iSlW/tx5GoCOtTKYG+qzzyAxERo2hBYtcjE7ybPFTLly5QgLC2Px4sWOfUlJSSxfvpzmzZtneJ63tzdBQUFONxERkduxet9ZYhNTKB7kTb3wdFbITkyETz6x3x882D4AWHKNqVczxcXFsW/fPsf2wYMHiYyMpEiRIkRERDB48GDeeecdKlWqRKVKlXjnnXfw8/Pj8ccfNzFrEREpaH7fltrFVLMEbm7pFCo//ACnTtkXk3z44VzOTkwtZjZu3EibNm0c20OGDAHgqaeeYtq0abzyyitcvnyZ/v37OybNW7RoEYGBgWalLCIiBUxSio3FO9OKmXS6mAwDJk6033/+efDKYFZgyTEWwzAMs5PISTExMQQHBxMdHa0uJxERybJle6LoNXUDIQHerH/1btyvbZnZuRNq1bIXMceOQdGi5iSaz2Tl+7vAr80kIiKSmQXb7a0y7WsWv76QAaheHQ4cgI0bVciYRMWMiIhIBpKtNv5I7WLqWDOTFbLLlNEaTCbKs1cziYiImG39gfNcjE+mqL8XjcsVuT4gPj73k5LrqJgRERHJwPx/TgJwb40wPNyv+cpMToZKlaBTJzh50oTsJI26mURERNJhtRn8kTrrb7oT5f3+O5w4ASkpGitjMrXMiIiIpOPvg+c5dymJQn6eNC2fTrHy5Zf2n7166XJsk6mYERERSUfaCtltqxXH89oupmPHYMEC+/2nn87lzORaKmZERESuYRgGS3bb12K6p3o66wFOnQo2G7RsCZUr53J2ci0VMyIiItfYFxXH0fOX8XJ3446KIc4HbTb46iv7/b59cz85uY6KGRERkWv8uSsKgGYViuLvfc21Mn/+CYcPQ6FC0K1b7icn19HVTCIiItdYsiu1i6la6PUHmzaFzz6DS5fA1zeXM5P0qJgRERG5yvlLSWw+cgGAu6qlM14mKAj+859czkoyo24mERGRqyzbE4XNgGolgihVSC0vrkDFjIiIyFWWpI6Xua6LyTCge3eYPNnexSR5hooZERGRVEkpNpb/ewaAu6peU8xERsIPP8CQIfZZfyXPUDEjIiKS6u+D54lLTCEkwJs6pQs5H/z+e/vPzp0hODjXc5OMqZgRERFJ9WfqVUx3VS2Gm5vlygGrFWbMsN9//HETMpPMqJgRERHBedbfu6+9imnlSjh+3D63TMeOuZ+cZErFjIiICFfN+uvhxp2Vrpn19//+z/7zoYfA2zv3k5NMqZgRERHhyqy/zSsUxc/rqmnYEhPhp5/s99XFlCepmBEREeHKrL/XdTGdOwd33gnlytkXlpQ8RzMAi4hIgXcx/qpZf6+9JLtkSfj1V0hKAnd3E7KTG1HLjIiIFHir9p3FZkDl4gEZz/rr5ZW7SclNUzEjIiIF3orUifJaVirmfGDrVjhwwISMJCtUzIiISIFmGIZj1t9WVa4pZoYOhQoV4MsvTchMbpaKGRERKdD+PR3H6ZhEfDzdaFS2yJUDp07BkiX2+23amJOc3BQVMyIiUqCldTE1KVcUH8+rBvj++CPYbNC0qb11RvIsFTMiIlKgObqYKl/TxZQ2t8wjj+RyRpJVKmZERKTAik9K4e+D5wFoeXUxc/q0fQkDgAcfNCEzyQoVMyIiUmCtP3CeJKuNUoV8qVDM/8qBX34Bw4CGDaFMGfMSlJuiYkZERAqstC6mlpWLYbFctUr2r7/af3brZkJWklWaAVhERAqsFY7xMtcsLDlzJixcCI0bm5CVZJWKGRERKZCOno/nwNlLuLtZaF7xmmImIMC+Qra4BHUziYhIgbRir71Vpn5EIYJ8PE3ORm6HihkRESmQlu9JZwmD2Fh719KYMfaFJcUlqJgREZECJ9lqY83+c8A1SxjMnw8bNsB334GnWmtchYoZEREpcDYfvkBcYgpF/L2oWTL4yoHZs+0/u3WDq69ukjxNxYyIiBQ4K/eeBeCOiiG4uaUWLZcvw++/2+9rojyXomJGREQKnFX7UouZSlddxbRoEVy6BBER9snyxGWomBERkQIl+nIy245dBOwtMw4//2z/+eCD6mJyMXm6mElJSeH111+nXLly+Pr6Ur58ecaMGYPNZjM7NRERcVHrD5zDZkD5EH9KFvK170xOht9+s9/XrL8uJ09Pmjdu3DimTJnC9OnTqVGjBhs3bqR3794EBwfzwgsvmJ2eiIi4oLSrmJpXLHpl58WL0KEDbNoEzZqZk5jcsjxdzKxdu5YuXbrQqVMnAMqWLcuMGTPYuHGjyZmJiIircoyXubqLqVgx+P57++KS6mJyOXm6m+mOO+5gyZIl/PvvvwBs3bqVVatW0bFjxwzPSUxMJCYmxukmIiICcDomgX1RcVgs0LR80esDVMi4pDzdMjNs2DCio6OpWrUq7u7uWK1W3n77bR577LEMzxk7dixvvvlmLmYpIiKuYnVqq0zNksEU8vOy7zx9Gs6cgRo1VMy4qDzdMjNr1iy+++47vv/+ezZv3sz06dN57733mD59eobnjBgxgujoaMft6NGjuZixiIjkZav32cfLtLi6i+mbb6BWLXjqKZOyktuVp1tmhg4dyvDhw3n00UcBqFWrFocPH2bs2LE8lcGHztvbG29v79xMU0REXIBhGKzZb2+ZaXH14N+0q5gaNzYhK8kOebplJj4+Hjc35xTd3d11abaIiGTZgbOXOBmdgJeHG43KFrHvPHcOVq+23+/c2bzk5Lbk6ZaZzp078/bbbxMREUGNGjXYsmULEyZMoE+fPmanJiIiLmZN6niZBhGF8fF0t+9csABsNns3U5kyJmYntyNPFzMfffQRI0eOpH///kRFRVGyZEmeffZZRo0aZXZqIiLiYq6Ml0mni0mtMi4tTxczgYGBTJo0iUmTJpmdioiIuDCr7erxMqmDf5OTYeFC+30VMy4tT4+ZERERyQ47TkQTk5BCoLcHtUoF23euXAkxMRAaqsG/Li5Pt8yIiIhkh7RZf5tWKIqHe+rf8S1a2Ftmzp4FN/1t78pUzIiISL63Jm28TIWrxst4e0O7diZlJNlJpaiIiORrCclWNhw6D1wzWZ7kGypmREQkX9ty5CKJKTaKBXpTMTTAvnPaNHj5Zdi61dTcJHuom0lERPK1tQfsXUzNyhfFkrb20vTpsGwZlC0LdeqYlptkD7XMiIhIvrZuf2oxkzZeJiYGVq2y3+/QwaSsJDupmBERkXzrcpKVLUcvAPaWGQCWLIGUFKhUCSpUMDE7yS4qZkREJN/adPgCyVaDEsE+lCnqZ9+5YIH9Z/v25iUm2UrFjIiI5FtrD9jnl3GMlzGMK7P+qosp31AxIyIi+dba1PEyTdPGy+zcCUePgo8PtG5tXmKSrXQ1k4iI5EuXElPYdiwauGq8zNGjEBYGdeuCr695yUm2UjEjIiL50oZD50mxGZQu7Et4kdTxMu3bw/HjcPGiqblJ9lI3k4iI5EtXzy/jxM0NihQxISPJKSpmREQkX7pufpnYWLDZTMxIckqWipnx48dz+fJlx/aKFStITEx0bMfGxtK/f//sy05EROQWxCQks/146niZtGLmlVfs42WmTzcxM8kJWSpmRowYQWxsrGP7vvvu4/jx447t+Ph4Pvvss+zLTkRE5BZsOHgemwFli/pRItj3yiXZZ85A0aI3fgBxKVkqZgzDyHRbREQkL1h7bRfTnj1w6BB4eUGbNuYlJjlCY2ZERCTfSRv82zRt8G/arL+tWoG/v0lZSU5RMSMiIvnKxfgkdp6MAa66kklLGORrWZ5n5ssvvyQgIACAlJQUpk2bRkhICIDTeBoREREzrD94HsOACsX8CQ3ygfh4WLHCflBLGORLWSpmIiIi+OKLLxzbYWFhfPvtt9fFiIiImGXdgWvGyyxfDomJEBEBVauamJnklCwVM4cOHcqhNERERLLHugPngavGy5QvD8OHQ1AQWCwmZiY5RcsZiIhIvnExPondp+zjZRzFTJUqMHasiVlJTsvSAOD169ezIG0QVapvvvmGcuXKERoayn/+8x+nSfRERERyU9p4mUqhAYQEeJudjuSSLBUzb7zxBtu2bXNsb9++naeffpp77rmH4cOH89tvvzFW1a+IiJhk3bWXZK9bB/Pnw6VLJmYlOS1LxUxkZCR33323Y3vmzJk0adKEL774giFDhvDhhx/yww8/ZHuSIiIiNyNtvEyT8qkLSX7wAXTqBOPGmZiV5LQsFTMXLlygePHiju3ly5fT/qpr9hs1asTRo0ezLzsREZGbdPV4mSblioLVCosX2w+2a2diZpLTslTMFC9enIMHDwKQlJTE5s2badasmeN4bGwsnp6e2ZuhiIjITUgbL1MxNIBigd6wZQucO2e/iqlxY7PTkxyUpWKmffv2DB8+nJUrVzJixAj8/Py48847Hce3bdtGhQoVsj1JERGRG7kyXia1i2nRIvvPu+4C/aGdr2Xp0uy33nqLBx98kFatWhEQEMC0adPw8vJyHP/666+59957sz1JERGRG1l/7fwyacWMvpfyvSwVM8WKFWPlypVER0cTEBCAu7u70/Eff/yRwMDAbE1QRETkRi7GJ7Hr6vEysbGwZo39oIqZfC9LxUyfPn1uKu7rr7++pWRERERuxd/XjpdZuBSSk+2z/2r4Q76XpWJm2rRplClThnr16mEYRk7lJCIikiVXljBIHS/Trh3s2gUnTpiYleSWLBUz/fr1Y+bMmRw4cIA+ffrwxBNPUKRIkZzKTURE5KakDf5tUi51vIzFYl9UUgtLFghZuppp8uTJnDx5kmHDhvHbb78RHh7OI488wh9//KGWGhERMYXTeJny+gO7IMpSMQPg7e3NY489xuLFi9m5cyc1atSgf//+lClThri4uJzIUUREJENp42UqFPMnNNAHZsyARx6BefPMTk1ySZaLmatZLBYsFguGYWCz2bIrJxERkZu27tpLsmfPhh9/hE2bTMxKclOWi5nExERmzJhB27ZtqVKlCtu3b+fjjz/myJEjBAQE5ESOIiIiGVp/8KrFJa1WWLLEfkCXZBcYWSpm+vfvT4kSJRg3bhz33Xcfx44d48cff6Rjx464ud1WI0+Gjh8/zhNPPEHRokXx8/Ojbt26bFK1LSIiQHR8MjtPXjVeZtMmuHABgoOhUSOTs5PckqWrmaZMmUJERATlypVj+fLlLF++PN242bNnZ0tyFy5coEWLFrRp04YFCxYQGhrK/v37KVSoULY8voiIuLb1B89hGFA+bbxM2sKSbdqAR5a+4sSFZelf+sknn8RiseRULtcZN24c4eHhTJ061bGvbNmyufb8IiKSt6WNl2mWNl7mzz/tP9u2NSkjMUOWJ83LTb/++ivt2rXj4YcfZvny5ZQqVYr+/fvTt2/fDM9JTEwkMTHRsR0TE5MbqYqIiAmuLC5ZFC5dgtWr7QdUzBQoOTPQJZscOHCATz/9lEqVKvHHH3/Qr18/Bg0axDfffJPhOWPHjiU4ONhxCw8Pz8WMRUQkt1w3v8zJk1C/vn0Jg4oVTc5OcpPFyMOz3Xl5edGwYUPWpC0WBgwaNIgNGzawdu3adM9Jr2UmPDyc6OhogoKCcjxnERHJHYt2nOI/326iYmgAfw5pdeVAYiJ4e5uXmGSLmJgYgoODb+r7O0+3zJQoUYLq1as77atWrRpHjhzJ8Bxvb2+CgoKcbiIikv+sdXQxXTPrrwqZAidPFzMtWrRgz549Tvv+/fdfypQpY1JGIiKSVzhNlhcfDxojWWDl6WLmxRdfZN26dbzzzjvs27eP77//ns8//5znn3/e7NRERMREF+OT2J02XqZcUfj5ZyhSBDK5QETyrzxdzDRq1Ig5c+YwY8YMatasyX//+18mTZpEjx49zE5NRERMtD51PaaKoQEUC/S2zy9jtUJIiNmpiQny/IxC9913H/fdd5/ZaYiISB6ydr99vEyz8kXBMK5MlqdLsgukPN0yIyIikh6n+WX++QdOnQJfX2jRwuTMxAwqZkRExKVcuJTE7lOxQOr8MmmtMi1b6kqmAkrFjIiIuJT1B+1XMVUKDSAkwPtKMaNVsgssFTMiIuJS0rqYmlUoap8gL23RY42XKbDy/ABgERGRqzmNl7Fa4f33Yf16qFnT5MzELGqZERERl3H+qvEyjcsVAT8/eO45mDYNLBZzkxPTqJgRERGX8fdBe6tM5eKp42VEUDEjIiIuJG1+mabli8L58/Dpp7Bvn8lZidlUzIiIiMtYe/V4mcWLoX9/ePBBk7MSs6mYERERl3AmNpF/T8cBqcXMokX2A7qKqcBTMSMiIi4hrVWmeokgivh5an4ZcVAxIyIiLmHt/rMANK9QFPbsgaNHwcsL7rzT5MzEbCpmRETEJaxJHfzbvGLRK60yd95pvzxbCjQVMyIikucduxDP4XPxuLtZaFS2iMbLiBMVMyIikuelXZJdu3QwgZ5usHq1/YDGywhazkBERFxAWjHTvEJRcHeHgwftazLVqWNyZpIXqJgREZE8zTCMK+NlKoTYdwYHw/33m5iV5CXqZhIRkTzt4NlLnIpJwMvdjQZlCpudjuRBKmZERCRPW53aKlO/TCF8YqOheXMYOdK+YrYIKmZERCSPuzK/TAj89ResXQuzZ9vHzoigYkZERPIwm81wHvybNr+MLsmWq6iYERGRPGv3qVguxCfj5+VO7VLBml9G0qViRkRE8qw1qV1MjcsVwevAPjh0yL6EQevWpuYleYuKGRERybOcupgWLrTvbNkS/P1NzEryGhUzIiKSJ6VYbaw/eB5IHfybVsy0b29iVpIXqZgREZE8afvxaOISUwj29aRaiSAIC4PChVXMyHU0A7CIiORJq/bax8s0LV8EdzcLTJ1qn1vGTX+HizMVMyIikietTC1m7qxU7MpOzS0j6VB5KyIieU5cYgqbj1wAoGWlYrBvHxiGyVlJXqViRkRE8px1+8+RYjMoU9SPiEtnoVIlqFgRkpLMTk3yIBUzIiKS56zcewaAOyuFwB9/2HeGhdnnmBG5hooZERHJc5zGy+iSbLkBFTMiIpKnHLsQz4Gzl3B3s9AsIgj+/NN+oF07cxOTPEvFjIiI5Clpl2TXDS9EUOQmiImBokWhQQOTM5O8SsWMiIjkKVe6mK4aL3PvvbosWzKkYkZERPIMq81g1b50xsuoi0kyoUnzREQkz9h+PJroy8kE+nhQp3QwvPUWLFigYkYypWJGRETyjJX/2i/Jbl6hKB7ubvYrmHQVk9yAuplERCTPWLkvnSUMRG7ApYqZsWPHYrFYGDx4sNmpiIhINotLTGHz4dQlDMoWghEjYNkysNlMzUvyPpcpZjZs2MDnn39O7dq1zU5FRERygNMSBru3wLvvwsMPa00muSGXKGbi4uLo0aMHX3zxBYULFzY7HRERyQFOSxj8/rt9Z4cOuiRbbsglipnnn3+eTp06cc8995idioiI5JC0+WXuqFgM5s2z77zvPhMzEleR569mmjlzJps3b2bDhg03FZ+YmEhiYqJjOyYmJqdSExGRbHLw7CUOnL2Eh5uFFpYLsHs3eHjYJ8sTuYE83TJz9OhRXnjhBb777jt8fHxu6pyxY8cSHBzsuIWHh+dwliIicrv+2h0FQONyRQj8c5F95x13QKFC5iUlLiNPFzObNm0iKiqKBg0a4OHhgYeHB8uXL+fDDz/Ew8MDq9V63TkjRowgOjracTt69KgJmYuISFYsTS1m7qoaemW8jLqY5Cbl6W6mu+++m+3btzvt6927N1WrVmXYsGG4pzMozNvbG29v79xKUUREblNcYgrrD54D4K6KReDff+0HVMzITcrTxUxgYCA1a9Z02ufv70/RokWv2y8iIq5p1d4zJFsNyhb1o3yJQnDgAGzbBpUrm52auIg8XcyIiEj+lzZepk3VUPsONzeoW9e8hMTluFwxs2zZMrNTEBGRbGKzGSzdY59f5u4qofbZft3y9HBOyYP0iREREdP8cyKaM7GJ+Hu50+TcAShRAgYNMjstcTEqZkRExDRpXUx3VArBc+F8iIqCEydMzkpcjYoZERExTVoxc3fV4pr1V26ZihkRETFFVGwC245FA3CXXzxs3gwWC3TsaHJm4mpUzIiIiCmWpQ78rVUqmJDFC+w777wTQkNNzEpckYoZERExxV+7rpr1d/Zs+84HHzQxI3FVKmZERCTXJaXYWLXPvkr2vUWBVavsB1TMyC1wuXlmRETE9W04dJ64xBRCArypFhYIw4bBwYOgxYHlFqiYERGRXPfHjlMAtKlSDLeSJWDsWJMzElembiYREclVNpvBwn/sxUyHWmEmZyP5gYoZERHJVZuPXCAqNpFAbw/uiPrXPr9MQoLZaYkLUzEjIiK5akFqq8zd1ULxmjQROndWN5PcFhUzIiKSawzjShfTfRWCYeFC+4GuXU3MSlydihkREck1245Fc/ziZfy83Gl5cDNcvgzlykGdOmanJi5MxYyIiOSatC6mNlVD8fp1rn3ngw/alzEQuUUqZkREJFcYhsGCf04C0LFKEfjtN/sBTZQnt0nFjIiI5IpdJ2M5fC4ebw837j7+D8TEQIkS0LSp2amJi1MxIyIiuSKtVaZV5WL4rF9r39m1K7jpq0huj2YAFhGRXJE2XqZjrRLw5Nvw+OPg62tyVpIfqJgREZEct/d0LPui4vB0t3BXtVD7zho1zE1K8g217YmISI5La5W5s1IxgtwMk7OR/EbFjIiI5Lj52+3jZe4v7QWhofDYY1rCQLKNihkREclRe07FsvtULJ7uFtr+swKio2H3bvDxMTs1ySdUzIiISI6avfkYAG2qhOL/0yz7zh49TMxI8hsVMyIikmNSrDbmbDkOQI8wA1avts/2+9hjJmcm+YmKGRERyTGr958jKjaRwn6etPh7kX1nmzZQqpS5iUm+omJGRERyzM+b7F1M99cugceM7+071cUk2UzFjIiI5IiYhGT+2GG/JPtxn4uwcyd4e0O3buYmJvmOJs0TEZEcsWD7SRJTbFQMDaByvcrw7rtw/jwEB5udmuQzKmZERCRH/LzJPvC3W/3SWIoXh2HDTM5I8it1M4mISLY7ci6evw+dx2KBB+qVNDsdyedUzIiISLabvcU+8PeOiiGUGPdf+P57SEw0OSvJr9TNJCIi2cowDGZvvmpumb7vgmFAw4ZQubLJ2Ul+pJYZERHJVhsPX+DI+Xj8vdy5e/Vv9kLmrrtUyEiOUTEjIiLZasbfRwDoVC0Ez6lf23f262diRpLfqZgREZFsczYukXlb7StkPxe7C06dsq+S3aWLyZlJfqZiRkREss2M9UdIstqoE16Icj99a9/59NPg5WVuYpKvqZgREZFskWy18d36wwA8Hw4sXmxfVLJvX3MTk3xPVzOJiEi2WPjPKU7HJBIS4E2bEj7QqhX4+kK5cmanJvlcnm6ZGTt2LI0aNSIwMJDQ0FAeeOAB9uzZY3ZaIiKSjmlrDgHQo0kEno0awLJlMGeOqTlJwZCni5nly5fz/PPPs27dOhYvXkxKSgr33nsvly5dMjs1ERG5yvZj0Ww6fAFPdws9mkRcOeDjY15SUmDk6W6mhQsXOm1PnTqV0NBQNm3aRMuWLU3KSkRErpXWKtOxZhih33wJ3btDsWLmJiUFRp5umblWdHQ0AEWKFDE5ExERSXM2LpHftp4AYKBxGAYOhOrVISHB5MykoMjTLTNXMwyDIUOGcMcdd1CzZs0M4xITE0m8av2PmJiY3EhPRKTAmvn3lcuxK3492r7z8cfVxSS5xmVaZgYMGMC2bduYMWNGpnFjx44lODjYcQsPD8+lDEVECp6kFBvfrrNfjj04+AIsWQIeHjBkiMmZSUHiEsXMwIED+fXXX1m6dCmlS5fONHbEiBFER0c7bkePHs2lLEVECp6fNh3jdEwixQK9aTk7demCxx+HMmXMTUwKlDzdzWQYBgMHDmTOnDksW7aMcjcxV4G3tzfe3t65kJ2ISMGWmGLl47/2AjCiHLiPnGs/MGyYeUlJgZSni5nnn3+e77//nl9++YXAwEBOnToFQHBwML6+viZnJyI35cJxWD4blv0Bp49Cx/JXjv24E+KS0j8v2BserHZle85uuJjBgFJ/T3ikxpXt3/6Fs/Hpx3q5Q49aV7YX7IVTGUz34G6BJ+tc2V58AI5lMg6vVx37jLcASw/BoYsZx/aoZc8FYOVh2Hch49hHa4Cvp/3+2mOw+2zGsQ9Vg8DUP+g2HId/zmQc26UKFEn9XbrlFESeyjj2vkpQzN9+f3sUbDxBSmIKP19Oxs3NQujCeKjgDr4+MG8QbAq0x+46A+uOZ/y4d5eDiGD7/b3nYdWRjGNblYHyhe33D16EZYcyjr0jHCoVtd8/Gg1/Hsw4tkkpqJ565dXJOFi4L+PYBiWgdnH7/TOXYN7ejGPrFIf6Jez3L1yGuZnMk1ajGDQuZb8fmwQ/7cw4tkpRaJ46hCIhBWb8k3FshcLQMrWVLMUG327LOLZMMNx1VaPB1MiMY0sFwr0VrmzX7Ab1nsg4Pofl6WLm008/BaB169ZO+6dOnUqvXr1yPyERuTmGAS/cB/OXwZF4SE7dX8INqhy6ErcyFi4a6T9GiBvUOXlle1UcnLGlHxtsgQanr2yvjoMTGcT6WqDpVV/way/BIWv6sR5Ai3NXttfFw96U9GMB7jx3pZhZHw87M4lteg68UmM3XIbI5IxjG5wF/9RRARsvw4ZMYmudhcKpsZsSYE0GxSJA1TMQmlpQbU6E5YkZx1Y4AzGpsZGJsDQRf8D/usA42LEW4lO/XrYmwdJMrmoKOw3JqYXathvEFj0NRmrszmRYejnj2MBT4Ja6HtTeG8T6nALv1NhDKbA0g0IYwP0k+KcWiyessDSTec+sJyB4l/3+mRvEJpyAoqnFzkUbLI3LODb2OBRPLaLibxB7/jiU2m+/n2zA0tiMY6t5QJmrir6lmRTuFT2gwuEr2yXqZhybC/J0MWMYGfySE5G8be3v8NH8K9s+7lClJNzdAro+cGV/wkKIy+CLo1AgdG17Zdv6J1zI4JdrgB90bX9l27IUzmTQ0uHtBV3vu7LtuQJOZtDS4e4OXa9a7dl3DRzNpPWia9crxUzgejiYSYvEg/eDZ+qv4MIbYW8mLRLdOtpbPABCI6HxgYxjH2oHQaklRsntUC+TloNud0OR1FaRiF1Qc1fGsfe3huKp02KU/5dDRdez40QMvp7utKpSDHeLBRITwdsbOt4BpULtsZX3Q7mtGT/uvc2gTGrrRfXDUHpTxrF3NYYKqeMm6xyH0PUZx7ZsAFVSWySOnILCazKObVYHaqa2Mpw4AwErM45tVAPqVrHfjzoPXssyjq1XFRpWt98/HwOWPzOOrVUJmqa2GMZcgpQ/Mo6tVg7uqGe/fzkR4n/POLZSBLRuaL+fnALRv2YcW7YktG16Zfvs7IxjSxeHDi2ubIdWzTg2F1iMfF4xxMTEEBwcTHR0NEFBQWanI1IwPNMWvvoTKhaBWfOhTkN7YSD5QkKylZbjlxIVm8hbD9TkiaYa7CvZLyvf3y5xNZOIuBDDgHUb7Pd7PAz1m6iQyWe+X3+EqNhEShXy5dH1v8LmzWanJAVcnu5mEhEXdGY3dLPBHcHw7Aizs5FsdjnJyuRl9jEYr5a14fHEi2Czwa5dUKWKydlJQaWWGRHJXjt/sY8badkeSqj7Ib/5dt0hzsYlUrqwLx2+nQhWK9x/vwoZMZWKGRHJXtvn2H9W75J5nLicqJgEPlxiv2z5rcDTuM2fb5/td/x4kzOTgk7FjIhkn01LYNgG+C0RKt1rdjaSzd76fRdxiSnUKxVIqy9SC5h+/aByZXMTkwJPxYyIZJ+vJsJlAy4HgX9Rs7ORbLRm31l+3XoCNwt8ZN2BJTISgoJg9GizUxPRAGARyUYLltt/3t/B3DwkWyWl2Hj9F/sss33qhVJ6yLP2A6+9BiEhJmYmYqdiRkSyx9aVcCh1JtL/DDc3F8lWX6w8wIEzlwgJ8GZgh5pw7jWYNg0GDTI7NRFA3Uwikl2+fM/+s0pRKF8t81hxGUfPx/NR6mKSr3WqSnCgLzz3HKxbBz4+JmcnYqdiRkSyx/xl9p+d22caJq5lzLydJCTbaB4RxAPlAq4cSFu2QSQPUDEjIrdvxzo4kLpu0n+GmZuLZJu5W46zeOdpPNwsfHhwAZZateDPTNYXEjGJxsyIyO3bvxju9AJrMahUy+xsJBscOBPHa3O2A/BGeBIhg/4HKSlw7twNzhTJfSpmROT2WJNh/w9wlw/c91+zs5FskJBsZcD3W7iUZKV5eCA9pgy0FzLdusEjj5idnsh11M0kIrdn+08QfQT8i0GdR83ORrLB2Pm72HkyhiL+Xnx+eCGWbdvsl2BPnqyxMpInqZgRkVtntcJzg2BPMjTtD56+Zmckt2nhP6eYvvYwAF+VvUTAhNSZfidPhtBQEzMTyZiKGRG5dV+OgXVnYXYCVHjQ7GzkNh27EM8rP20FYGhVb+q9+Iy9YH3iCXj4YZOzE8mYxsyIyK2x2eB/k+z3uzSHkmXNzEZuU2xCMs9+u4mYhBTqhhfiP93qwMJ7YN8++Pxzs9MTyZSKGRG5NTM/gP0x4A6M/czsbOQ2JKZYefbbTew4EUNIgBcfPVYPz0A/mDEDoqPBV92Hkrepm0lEbs3Yd+w/29WDCjXMzUVumdVm8OKsSNbsP0eAtwc/VkkkvFDqzL4WCxQqZGp+IjdDxYyIZN28qfDPWftvkHc+MTsbuUWGYTD613+Yv/0UXu5u/Fz0KOUe6Qzdu9vHyoi4CBUzIpJ1Y0baf7asAnWamZuL3LIPluzlu3VHsFhgRtHjVHnpOfuBiAhwdzc3OZEsUDEjIllzcAWUPQ8l3OHtSWZnI7fAMAwmLv6XSX/aF5CcFnCIBq/0s0+M16MHjB9vcoYiWaMBwCJy8xKiYc5zUN0TnugLzbWopKtJttp4bc52fth4DIApHntpNfple7fSk0/C11+rVUZcjooZEbl5M1+AmGNQuCy0e9vsbCSL4pNSeP7/NrN0zxncLPCN337ueOMl+2X2vXrBl1+qkBGXpGJGRG7OlNEwaDrc6wuTPwPvALMzkiw4G5fI09M2sPVYND6ebnz0WH3uOOwGHh7Qs6d9Lhk3jTwQ16RiRkRu7NBOeOVtSAa8akNEU7MzkixYf+AcL86K5ER0AoV9Pfiqd2PqRxSG6h3g77+hVi0VMuLS9OkVkczZbPBIO4i1QpgvTP/D7IzkJiVbbby/aA+PfrGOE9EJtOAia+eNon7siStBdeqokBGXp5YZEcncG0/DhmP2P32++QYCg83OSG7CkXPxDJq5hcijFwEYzX56fTYay8WL8NxzsGKFqfmJZCcVMyKSsf8Nhrem2e8/3QnaPmRmNnITkq02vll7mImL/yUuMYUIaxwzds6i1Pw59oBmzWDmTHOTFMlmKmZEJH0//w+GfQAGcFc1+GSO2RnJDSzbE8V/5+1k/5lLYBi8FPU3/X+djPv5c/aupMGD4Z13wNvb7FRFspWKGRG5XuT3sP1tuNMLPCvCgkjw9DQ7K8nAvqg43pm/i792RwFQ1N+LST6HuHP8f+0BderAF19Ao0YmZimSc1TMiIizzd/BrwMAA17uDx3f09wjedT2Y9FMWb6f+f+cxDAgKCWB7m2qMfDuSgR5ucOiH+Duu+Gll1SMSr6mYkZE7OKi4ZmO8Pd6eNQPWvS1FzIWi9mZyVUMw2DN/nN8umw/q/adBSDiwknG7vyFJns34vHGXvBJLVwWLdK/nxQIKmZEBP6cBb16w/HL9u3oZipk8pjTMQn8vPkYP206xoEzl3CzWbnn4GaG7F9CtcjVWAzDHvj77/ZVr0H/flJgqJgRKciSEmBAV/h6IVgBfzd491UY8F+zMxMgJiGZpbujmLPlOCv+PYPNgKCEOAZtW0ifHYsoFHXVfDEdOsDo0dCkiXkJi5hExYxIQWRNhveGwIdfwYnU1piGpWHWIihfzdzcCrgTFy/z567TLN55mnUHzpFsNfBKScbm4UmjsoXpWSaMzh9Ot7fEFC4MffpAv35QsaLZqYuYRsWMSEGSFA9bvoM1H8I3e+BECnhb4LXn4LWPNBOsCaJiE1h/4DzrDpxj3YFz7D9zCe+UJOof38XgQ5HcfWwbviGFsS3+k3Ih/vaTBg2CunXt3Um+vqbmL5IXqJgRye+SEuHnKfD9NCgbBUXj7PvvLQHNa8KYT6FEGVNTLChiE5LZcSKG7cei2X48mm3HLnLoXDwATY5sp/ORbTQ8votGx3finZx05cTTnuCZcmV70qTcTVwkj3OJYmby5Mn873//4+TJk9SoUYNJkyZx5513mp2WSN51cCfMnwG/zYHVuyDOZt/f1Au6V4YWL0DdHuCpv+qzm81mcCYukaPn49l/Jo59UXHsjbL/PHkujrIXTlDp7BEqnT3CvGaPYHF3p3qJIEatWkaN1Vete1WyJNxzj/3Wti0EaxkJkYzk+WJm1qxZDB48mMmTJ9OiRQs+++wzOnTowM6dO4mIiDA7PRFz2Wxw7ihcOgpRu2DXShg6C84kOcf5ukGzytCrDzz+Irjn+f/6edLlJCvn45M4H5dEVGwCp2MSOR2TQFRsAicuJnDqzEUOxSSTaLXHt96/kbv2b6DFxVOER5+mdPQpvK1XWljuHDmQKnfUJ9jPE4Ifg/Kh0Lgx3HUXVKumq5FEbpLFMNKu58ubmjRpQv369fn0008d+6pVq8YDDzzA2LFjb3h+TEwMwcHBREdHExQUlJOpimQvmw0SY+DyeVizEvbvgiMH4fgxOBkFx8/AiVgoZYEeqWMpDAPej4NLhn2F67qV4NEnoPsA8Ck4rTCGYZBiM0hKsZGUYiMxxUZiipXEFBuXk6wkJFu5nGz/GZ9kJT4+kYS4eBLj4kmIiycpNo7k6FhSYuNIib3E+rDKnEh253KyleaHIml5cDOFL8dSOCGWwvExFL4cQ9H4aAonxHLXM1M4FFKaEsG+vLzme7r+9pVzcn5+UL061KgBr7+ugbsiGcjK93ee/vMsKSmJTZs2MXz4cKf99957L2vWrDEpK7uYi+e4vOcf3E+cyDAmuXwFKGRvGnaLOov7sSMZxlojymILKWKPPXse9yOHMo4tHYEtNAQAy/kLeBw6mElsaWyhofaNi9F4HtifcWzJktjCwuwbcXF4/vtvxrHFi2MrVcq+cSkezz27M44NLYatdLh94/JlPHftcg64qp62hoRgK5M6fiMpCc/t26+LSWMrHIy1fAX7RkoKXpu3pPPs9u4VW3Aw1kqVAANsNjw3bHA8psWa+me0YQPDwBoUhLV6VXusYeC9ajUkJ9vjbTYwrPZzrFZswQGkNKyDxbBiWK34/L4IS0ICFmsyluQUSE6ClCQsyUkYhfxIurcRlpQELCkJ+HyxALeYS1gSErFcToSEJCwJyVgSUiDEHR71wWKk5jYhFmLT/7vDOGshKTCCxEIViA+pg/H6JS6H1iI5OBRsBobVCnMXY6RYsXp4EtOoqf2lGAaBm/7G/cJ5sNmwpVjBZsOw2sBmxermzvG292EYBjbDoNjKpficOp4aY8Ww2SDFCjYrNgO2PdIHq83+uOWXzKPw4X32mNTHw2oFm/09XPjkiyRZ3LHZDOovmUPEv9vs723qzbDZsNhsYNj4oudwLnn4YLUZtF36E3V3rQebYf83MFLjrDYsho0Xu73KaZ8gkq02nl4/mwd3LMXdZsXDZsPNsOJntRJkS8HDZuWBnu9zrJD98z582VT6rf85w8/wvX0+5nKxsgA0ObGLfn/PzjB2RsdwitzfHk93N1jhAbWKQ7ly9luFChARoYHWItksTxczZ8+exWq1Urx4caf9xYsX59SpU+mek5iYSGJiomM7JiYmR3Lb8csEms2cAAsTMg7q4QcVU9/izUnwWyaxD/tC9dRZO/9Jhp8vZxz7gA/U8bLf35MMMzOJ7eQDDVNjD6TAt/EZx7b1huapC9AdS4GvMolt7Q2tUmNPW2HKpYxjW3jBPT72++dt8FFcxrGNPKFjagtCnM3eypCRup7QJTU20YB3YzOOreEBD/nZ79sM+G8msZU94DG/K9tvxdjnYElPOXd40v/K9rgYyOifubQ7/oF/XtlekXGBAm5g2D8PlwxvfIvE4+ZjAy/A3QIWIMWAeIOjlmK0PPMunAH2wu9TB1Ej6v10HzXKvzD3DfjWsf3jd69Q9/jOdGNjvXzpdraUY/ubWRNpeCi9ghGsFjceDWju2J4y5wea/rs2g9cGT1XqSpKH/fU1WbmcZjuWZhjbr0lvon0DAXjw3100+SfjP2Ss8fEkeQQAEBZ7jupRGRf6NYt4UTy8MH5e7pTeU8Sx3+bmRoqPL1ZfP2y+fhh+fnzwSB38G9ansL8nASu9YGFxKFIEQkKgaFH7rVgxKFmS4oULX+keatnSfhORHJWni5k0lmv6jQ3DuG5fmrFjx/Lmm2/mfE7uHlgtFtwyeQdTDHesqV9IHpaUTIcppODmiHXHikcmsVbcSXHE2m4Q63ZVrJFprM3iRnJarMENYi2OWDfA8waxSYZXamwKXpnEGm5uJKbGYtgcs7KnH2vhsmEvqCyGDV9P7F/y6cW6uxFn2AsUwzAI8rmmmLnq82TzcOO8EYyBGwYQGhCb1sBjf3zLlfspPh4cMCKw4YYNC9VCdmBJsoGbveAw0uItFhIK+fCH0ZpEvEjEi4dK/4pXgn1si82w2BuNDDcMG5zwL8ZzlnHEEECyuzdfGS8SfvkEXLZgD7VgWOwPfqFwMcKCfLCk1jjni4Zx1Jpof34sGG5uGBYLBhai/YOpXDwAt9TXezaiAns8DGxu7vYYixtY7OckevvQrHxR3NzAzWLhXO0GbC4cCG5uGBY3DDd3DDcLhrs7uLnTtV4p3CwW3N0g/uI9rCsfYV/Tyc3NfnN3t7/P7u70v7sSFk8v3N3AN+BhNhyqiyU1xuJmwc3DA4u7OxY3N956qCEWPz883NwoWuM59hxui5ubG24eHrh5uuPm7o6bpwfu7u7MaNcBj+BAPN0t+HQPJ/HEs3h4e+Hu5Wn/QLu7g5cXeHkxpWLFK6tHP14bbBPB2xs3Dw+8rvn8OM2807at/SYieUaeHjOTlJSEn58fP/74I127dnXsf+GFF4iMjGT58uXXnZNey0x4eLjGzIiIiLiQrIyZydMdt15eXjRo0IDFixc77V+8eDHNmzdP9xxvb2+CgoKcbiIiIpJ/5flupiFDhtCzZ08aNmxIs2bN+Pzzzzly5Aj9+vUzOzURERHJA/J8MdO9e3fOnTvHmDFjOHnyJDVr1mT+/PmUSbviRURERAq0PD1mJjtonhkRERHXk2/GzIiIiIjciIoZERERcWkqZkRERMSlqZgRERERl6ZiRkRERFyaihkRERFxaSpmRERExKWpmBERERGXpmJGREREXJqKGREREXFpeX5tptuVtlpDTEyMyZmIiIjIzUr73r6ZVZfyfTETGxsLQHh4uMmZiIiISFbFxsYSHBycaUy+X2jSZrNx4sQJAgMDsVgs2frYMTExhIeHc/To0Xy5iKVen+vL769Rr8/15ffXqNd36wzDIDY2lpIlS+LmlvmomHzfMuPm5kbp0qVz9DmCgoLy5Yc0jV6f68vvr1Gvz/Xl99eo13drbtQik0YDgEVERMSlqZgRERERl6Zi5jZ4e3szevRovL29zU4lR+j1ub78/hr1+lxffn+Nen25I98PABYREZH8TS0zIiIi4tJUzIiIiIhLUzEjIiIiLk3FTDZLTEykbt26WCwWIiMjzU4n29x///1ERETg4+NDiRIl6NmzJydOnDA7rWxz6NAhnn76acqVK4evry8VKlRg9OjRJCUlmZ1atnn77bdp3rw5fn5+FCpUyOx0btvkyZMpV64cPj4+NGjQgJUrV5qdUrZasWIFnTt3pmTJklgsFubOnWt2Stlm7NixNGrUiMDAQEJDQ3nggQfYs2eP2Wllq08//ZTatWs75l9p1qwZCxYsMDutHDN27FgsFguDBw825flVzGSzV155hZIlS5qdRrZr06YNP/zwA3v27OHnn39m//79PPTQQ2anlW12796NzWbjs88+Y8eOHUycOJEpU6bw6quvmp1atklKSuLhhx/mueeeMzuV2zZr1iwGDx7Ma6+9xpYtW7jzzjvp0KEDR44cMTu1bHPp0iXq1KnDxx9/bHYq2W758uU8//zzrFu3jsWLF5OSksK9997LpUuXzE4t25QuXZp3332XjRs3snHjRu666y66dOnCjh07zE4t223YsIHPP/+c2rVrm5eEIdlm/vz5RtWqVY0dO3YYgLFlyxazU8oxv/zyi2GxWIykpCSzU8kx48ePN8qVK2d2Gtlu6tSpRnBwsNlp3JbGjRsb/fr1c9pXtWpVY/jw4SZllLMAY86cOWankWOioqIMwFi+fLnZqeSowoULG19++aXZaWSr2NhYo1KlSsbixYuNVq1aGS+88IIpeahlJpucPn2avn378u233+Ln52d2Ojnq/Pnz/N///R/NmzfH09PT7HRyTHR0NEWKFDE7DblGUlISmzZt4t5773Xaf++997JmzRqTspLbER0dDZBv/79ZrVZmzpzJpUuXaNasmdnpZKvnn3+eTp06cc8995iah4qZbGAYBr169aJfv340bNjQ7HRyzLBhw/D396do0aIcOXKEX375xeyUcsz+/fv56KOP6Nevn9mpyDXOnj2L1WqlePHiTvuLFy/OqVOnTMpKbpVhGAwZMoQ77riDmjVrmp1Ottq+fTsBAQF4e3vTr18/5syZQ/Xq1c1OK9vMnDmTzZs3M3bsWLNTUTGTmTfeeAOLxZLpbePGjXz00UfExMQwYsQIs1POkpt9fWmGDh3Kli1bWLRoEe7u7jz55JMYeXzOxay+RoATJ07Qvn17Hn74YZ555hmTMr85t/L68guLxeK0bRjGdfsk7xswYADbtm1jxowZZqeS7apUqUJkZCTr1q3jueee46mnnmLnzp1mp5Utjh49ygsvvMB3332Hj4+P2eloBuDMnD17lrNnz2YaU7ZsWR599FF+++03p1+kVqsVd3d3evTowfTp03M61Vtys68vvQ/qsWPHCA8PZ82aNXm62TSrr/HEiRO0adOGJk2aMG3atBsuO2+2W/k3nDZtGoMHD+bixYs5nF3OSEpKws/Pjx9//JGuXbs69r/wwgtERkayfPlyE7PLGRaLhTlz5vDAAw+YnUq2GjhwIHPnzmXFihWUK1fO7HRy3D333EOFChX47LPPzE7lts2dO5euXbvi7u7u2Ge1WrFYLLi5uZGYmOh0LKd55NozuaCQkBBCQkJuGPfhhx/y1ltvObZPnDhBu3btmDVrFk2aNMnJFG/Lzb6+9KTVwImJidmZUrbLyms8fvw4bdq0oUGDBkydOjXPFzJwe/+GrsrLy4sGDRqwePFip2Jm8eLFdOnSxcTM5GYZhsHAgQOZM2cOy5YtKxCFDNhfd17/nXmz7r77brZv3+60r3fv3lStWpVhw4blaiEDKmayRUREhNN2QEAAABUqVKB06dJmpJSt/v77b/7++2/uuOMOChcuzIEDBxg1ahQVKlTI060yWXHixAlat25NREQE7733HmfOnHEcCwsLMzGz7HPkyBHOnz/PkSNHsFqtjnmQKlas6PjMuoohQ4bQs2dPGjZsSLNmzfj88885cuRIvhrjFBcXx759+xzbBw8eJDIykiJFilz3O8fVPP/883z//ff88ssvBAYGOsY6BQcH4+vra3J22ePVV1+lQ4cOhIeHExsby8yZM1m2bBkLFy40O7VsERgYeN0Yp7QxlaaMfTLlGqp87uDBg/nq0uxt27YZbdq0MYoUKWJ4e3sbZcuWNfr162ccO3bM7NSyzdSpUw0g3Vt+8dRTT6X7+pYuXWp2arfkk08+McqUKWN4eXkZ9evXz3eX9S5dujTdf6+nnnrK7NRuW0b/16ZOnWp2atmmT58+js9nsWLFjLvvvttYtGiR2WnlKDMvzdaYGREREXFpeX9QgIiIiEgmVMyIiIiIS1MxIyIiIi5NxYyIiIi4NBUzIiIi4tJUzIiIiIhLUzEjIiIiLk3FjIiIiLg0FTMiBVTr1q0ZPHiw2Wmk69y5c4SGhnLo0CEAli1bhsViyfHFMW/1eaZNm0ahQoWydE6jRo2YPXt2ls4RkfSpmBGRbHHy5Ekef/xxqlSpgpubW4aF0s8//0z16tXx9vamevXqzJkz57qYsWPH0rlzZ8qWLZuzSZto5MiRDB8+HJvNZnYqIi5PxYyIZIvExESKFSvGa6+9Rp06ddKNWbt2Ld27d6dnz55s3bqVnj178sgjj7B+/XpHzOXLl/nqq6945plncit1U3Tq1Ino6Gj++OMPs1MRcXkqZkSECxcu8OSTT1K4cGH8/Pzo0KEDe/fudYr54osvCA8Px8/Pj65duzJhwgSnrpWyZcvywQcf8OSTTxIcHJzu80yaNIm2bdsyYsQIqlatyogRI7j77ruZNGmSI2bBggV4eHhkuiL7uXPneOyxxyhdujR+fn7UqlWLGTNmOMW0bt2agQMHMnjwYAoXLkzx4sX5/PPPuXTpEr179yYwMJAKFSqwYMGC6x5/9erV1KlTBx8fH5o0acL27dudjk+bNo2IiAjHe3Hu3Dmn4/v376dLly4UL16cgIAAGjVqxJ9//ukU4+7uTseOHa/LW0SyTsWMiNCrVy82btzIr7/+ytq1azEMg44dO5KcnAzYv9z79evHCy+8QGRkJG3btuXtt9/O8vOsXbuWe++912lfu3btWLNmjWN7xYoVNGzYMNPHSUhIoEGDBsybN49//vmH//znP/Ts2dOphQdg+vTphISE8PfffzNw4ECee+45Hn74YZo3b87mzZtp164dPXv2JD4+3um8oUOH8t5777FhwwZCQ0O5//77He/F+vXr6dOnD/379ycyMpI2bdrw1ltvOZ0fFxdHx44d+fPPP9myZQvt2rWjc+fOHDlyxCmucePGrFy58ubePBHJmClrdYuI6Vq1amW88MILxr///msAxurVqx3Hzp49a/j6+ho//PCDYRiG0b17d6NTp05O5/fo0cMIDg7O9LGv5enpafzf//2f077/+7//M7y8vBzbXbp0Mfr06eMUs3TpUgMwLly4kOHr6dixo/HSSy855XDHHXc4tlNSUgx/f3+jZ8+ejn0nT540AGPt2rVOzzNz5kxHzLlz5wxfX19j1qxZhmEYxmOPPWa0b9/e6bm7d++e4XuRpnr16sZHH33ktO+XX34x3NzcDKvVmum5IpI5tcyIFHC7du3Cw8ODJk2aOPYVLVqUKlWqsGvXLgD27NlD48aNnc67dvtmWSwWp23DMJz2Xb58GR8fn0wfw2q18vbbb1O7dm2KFi1KQEAAixYtuq7lo3bt2o777u7uFC1alFq1ajn2FS9eHICoqCin867u4ipSpIjTe7Fr167rusCu3b506RKvvPIK1atXp1ChQgQEBLB79+7r8vP19cVms5GYmJjp6xWRzHmYnYCImMswjAz3pxUZ1xYcmZ2XmbCwME6dOuW0LyoqylFUAISEhHDhwoVMH+f9999n4sSJTJo0iVq1auHv78/gwYNJSkpyivP09HTatlgsTvvSXtPNXFF09XtxI0OHDuWPP/7gvffeo2LFivj6+vLQQw9dl9/58+fx8/PD19f3ho8pIhlTy4xIAVe9enVSUlKcxpucO3eOf//9l2rVqgFQtWpV/v77b6fzNm7cmOXnatasGYsXL3bat2jRIpo3b+7YrlevHjt37sz0cVauXEmXLl144oknqFOnDuXLl79uwPLtWLduneP+hQsX+Pfff6latSpgf7+uPn5tfFp+vXr1omvXrtSqVYuwsDDHnDlX++eff6hfv3625S1SUKmYESngKlWqRJcuXejbty+rVq1i69atPPHEE5QqVYouXboAMHDgQObPn8+ECRPYu3cvn332GQsWLLiutSYyMpLIyEji4uI4c+YMkZGRToXJCy+8wKJFixg3bhy7d+9m3Lhx/Pnnn05z0rRr144dO3Zk2jpTsWJFFi9ezJo1a9i1axfPPvvsdS0+t2PMmDEsWbKEf/75h169ehESEsIDDzwAwKBBg1i4cCHjx4/n33//5eOPP2bhwoXX5Td79mwiIyPZunUrjz/+eLqtPytXrrxuQLSIZJ2KGRFh6tSpNGjQgPvuu49mzZphGAbz5893dMm0aNGCKVOmMGHCBOrUqcPChQt58cUXrxvbUq9ePerVq8emTZv4/vvvqVevHh07dnQcb968OTNnzmTq1KnUrl2badOmMWvWLKfxOrVq1aJhw4b88MMPGeY7cuRI6tevT7t27WjdujVhYWGOYiM7vPvuu7zwwgs0aNCAkydP8uuvv+Ll5QVA06ZN+fLLL/noo4+oW7cuixYt4vXXX3c6f+LEiRQuXJjmzZvTuXNn2rVrd10LzPHjx1mzZg29e/fOtrxFCiqLcSsd3yJS4PXt25fdu3fnyKXF8+fP5+WXX+aff/7BzS1//s01dOhQoqOj+fzzz81ORcTlaQCwiNyU9957j7Zt2+Lv78+CBQuYPn06kydPzpHn6tixI3v37uX48eOEh4fnyHOYLTQ0lJdfftnsNETyBbXMiMhNeeSRR1i2bBmxsbGUL1+egQMH0q9fP7PTEhFRMSMiIiKuLX92RouIiEiBoWJGREREXJqKGREREXFpKmZERETEpamYEREREZemYkZERERcmooZERERcWkqZkRERMSlqZgRERERl/b/Zz0Lo7z18f8AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "import os\n", "import numpy as np\n", @@ -1974,7 +1690,7 @@ "np.random.seed(3155)\n", "\n", "x = np.random.rand(100)\n", - "y = 2.0+5*x*x+0.1*np.random.randn(100)\n", + "y = 2.0+5*x*x#+0.1*np.random.randn(100)\n", "\n", "# number of features p (here degree of polynomial\n", "p = 3\n", @@ -2010,7 +1726,7 @@ " lmb = lambdas[i]\n", " Ridgebeta = np.linalg.inv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train\n", " # include lasso using Scikit-Learn\n", - " RegLasso = linear_model.Lasso(lmb,fit_intercept=False)\n", + " RegLasso = linear_model.Lasso(lmb,fit_intercept=True)\n", " RegLasso.fit(X_train,y_train)\n", " # and then make the prediction\n", " ytildeRidge = X_train @ Ridgebeta\n", @@ -2038,9 +1754,7 @@ { "cell_type": "markdown", "id": "9d2563fa", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Material for lecture Thursday September 7" ] @@ -2048,9 +1762,7 @@ { "cell_type": "markdown", "id": "1232037e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Linking the regression analysis with a statistical interpretation\n", "\n", @@ -2077,9 +1789,7 @@ { "cell_type": "markdown", "id": "e50205be", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\begin{align*} \n", @@ -2093,9 +1803,7 @@ { "cell_type": "markdown", "id": "70eccda9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "The randomness of $\\varepsilon_i$ implies that\n", "$\\mathbf{y}_i$ is also a random variable. In particular,\n", @@ -2112,9 +1820,7 @@ { "cell_type": "markdown", "id": "e6bdcbf2", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Assumptions made\n", "\n", @@ -2126,9 +1832,7 @@ { "cell_type": "markdown", "id": "d74471b8", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{y} = f(\\boldsymbol{x})+\\boldsymbol{\\varepsilon}\n", @@ -2138,9 +1842,7 @@ { "cell_type": "markdown", "id": "f7b178d9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We approximate this function with our model from the solution of the linear regression equations, that is our\n", "function $f$ is approximated by $\\boldsymbol{\\tilde{y}}$ where we want to minimize $(\\boldsymbol{y}-\\boldsymbol{\\tilde{y}})^2$, our MSE, with" @@ -2149,9 +1851,7 @@ { "cell_type": "markdown", "id": "8ea16c62", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{\\tilde{y}} = \\boldsymbol{X}\\boldsymbol{\\beta}.\n", @@ -2161,9 +1861,7 @@ { "cell_type": "markdown", "id": "4a8d7bfd", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Expectation value and variance\n", "\n", @@ -2173,9 +1871,7 @@ { "cell_type": "markdown", "id": "5367d943", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\begin{align*} \n", @@ -2189,9 +1885,7 @@ { "cell_type": "markdown", "id": "c7eefd85", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "while\n", "its variance is" @@ -2200,9 +1894,7 @@ { "cell_type": "markdown", "id": "27ce343a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\begin{align*} \\mbox{Var}(y_i) & = \\mathbb{E} \\{ [y_i\n", @@ -2223,9 +1915,7 @@ { "cell_type": "markdown", "id": "c67ea20b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Hence, $y_i \\sim \\mathcal{N}( \\mathbf{X}_{i, \\ast} \\, \\boldsymbol{\\beta}, \\sigma^2)$, that is $\\boldsymbol{y}$ follows a normal distribution with \n", "mean value $\\boldsymbol{X}\\boldsymbol{\\beta}$ and variance $\\sigma^2$ (not be confused with the singular values of the SVD)." @@ -2234,9 +1924,7 @@ { "cell_type": "markdown", "id": "b4da473a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Expectation value and variance for $\\boldsymbol{\\beta}$\n", "\n", @@ -2246,9 +1934,7 @@ { "cell_type": "markdown", "id": "af812d38", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\mathbb{E}(\\boldsymbol{\\hat{\\beta}}) = \\mathbb{E}[ (\\mathbf{X}^{\\top} \\mathbf{X})^{-1}\\mathbf{X}^{T} \\mathbf{Y}]=(\\mathbf{X}^{T} \\mathbf{X})^{-1}\\mathbf{X}^{T} \\mathbb{E}[ \\mathbf{Y}]=(\\mathbf{X}^{T} \\mathbf{X})^{-1} \\mathbf{X}^{T}\\mathbf{X}\\boldsymbol{\\beta}=\\boldsymbol{\\beta}.\n", @@ -2258,9 +1944,7 @@ { "cell_type": "markdown", "id": "9d3aaff5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "This means that the estimator of the regression parameters is unbiased.\n", "\n", @@ -2272,9 +1956,7 @@ { "cell_type": "markdown", "id": "eebea8b7", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\begin{eqnarray*}\n", @@ -2303,9 +1985,7 @@ { "cell_type": "markdown", "id": "714f281a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "where we have used that $\\mathbb{E} (\\mathbf{Y} \\mathbf{Y}^{T}) =\n", "\\mathbf{X} \\, \\boldsymbol{\\beta} \\, \\boldsymbol{\\beta}^{T} \\, \\mathbf{X}^{T} +\n", @@ -2325,9 +2005,7 @@ { "cell_type": "markdown", "id": "66015ca6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\boldsymbol{\\beta}^{\\mathrm{OLS}}.\n", @@ -2337,9 +2015,7 @@ { "cell_type": "markdown", "id": "267f01c0", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We see clearly that \n", "$\\mathbb{E} \\big[ \\boldsymbol{\\beta}^{\\mathrm{Ridge}} \\big] \\not= \\boldsymbol{\\beta}^{\\mathrm{OLS}}$ for any $\\lambda > 0$. We say then that the ridge estimator is biased.\n", @@ -2350,9 +2026,7 @@ { "cell_type": "markdown", "id": "998de4b6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\mbox{Var}[\\boldsymbol{\\beta}^{\\mathrm{Ridge}}]=\\sigma^2[ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1} \\mathbf{X}^{T} \\mathbf{X} \\{ [ \\mathbf{X}^{\\top} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}\\}^{T},\n", @@ -2362,9 +2036,7 @@ { "cell_type": "markdown", "id": "06a8d452", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and it is easy to see that if the parameter $\\lambda$ goes to infinity then the variance of Ridge parameters $\\boldsymbol{\\beta}$ goes to zero. \n", "\n", @@ -2374,9 +2046,7 @@ { "cell_type": "markdown", "id": "6cd3ed9b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\mbox{Var}[\\boldsymbol{\\beta}^{\\mathrm{OLS}}]-\\mbox{Var}(\\boldsymbol{\\beta}^{\\mathrm{Ridge}})=\\sigma^2 [ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}[ 2\\lambda\\mathbf{I} + \\lambda^2 (\\mathbf{X}^{T} \\mathbf{X})^{-1} ] \\{ [ \\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I} ]^{-1}\\}^{T}.\n", @@ -2386,9 +2056,7 @@ { "cell_type": "markdown", "id": "70f9aa69", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "The difference is non-negative definite since each component of the\n", "matrix product is non-negative definite. \n", @@ -2398,9 +2066,7 @@ { "cell_type": "markdown", "id": "3bb18388", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Deriving OLS from a probability distribution\n", "\n", @@ -2421,9 +2087,7 @@ { "cell_type": "markdown", "id": "84b546af", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "y_i\\sim \\mathcal{N}(\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta}, \\sigma^2)=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}.\n", @@ -2433,9 +2097,7 @@ { "cell_type": "markdown", "id": "0e90935e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Independent and Identically Distrubuted (iid)\n", "\n", @@ -2446,9 +2108,7 @@ { "cell_type": "markdown", "id": "b40b33f5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(y_i, \\boldsymbol{X}\\vert\\boldsymbol{\\beta})=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]},\n", @@ -2458,9 +2118,7 @@ { "cell_type": "markdown", "id": "23e0c7ac", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which reads as finding the likelihood of an event $y_i$ with the input variables $\\boldsymbol{X}$ given the parameters (to be determined) $\\boldsymbol{\\beta}$.\n", "\n", @@ -2470,9 +2128,7 @@ { "cell_type": "markdown", "id": "dcd4ed9c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{y},\\boldsymbol{X}\\vert\\boldsymbol{\\beta})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}=\\prod_{i=0}^{n-1}p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta}).\n", @@ -2482,9 +2138,7 @@ { "cell_type": "markdown", "id": "5c37adb5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We will write this in a more compact form reserving $\\boldsymbol{D}$ for the domain of events, including the ouputs (targets) and the inputs. That is\n", "in case we have a simple one-dimensional input and output case" @@ -2493,9 +2147,7 @@ { "cell_type": "markdown", "id": "27e4eea9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\\dots, (x_{n-1},y_{n-1})].\n", @@ -2505,9 +2157,7 @@ { "cell_type": "markdown", "id": "5dd95110", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "In the more general case the various inputs should be replaced by the possible features represented by the input data set $\\boldsymbol{X}$. \n", "We can now rewrite the above probability as" @@ -2516,9 +2166,7 @@ { "cell_type": "markdown", "id": "bc7cfda7", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{D}\\vert\\boldsymbol{\\beta})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}.\n", @@ -2528,9 +2176,7 @@ { "cell_type": "markdown", "id": "2bf6dbf4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "It is a conditional probability (see below) and reads as the likelihood of a domain of events $\\boldsymbol{D}$ given a set of parameters $\\boldsymbol{\\beta}$." ] @@ -2538,9 +2184,7 @@ { "cell_type": "markdown", "id": "e8ff9386", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Maximum Likelihood Estimation (MLE)\n", "\n", @@ -2569,9 +2213,7 @@ { "cell_type": "markdown", "id": "cd665e76", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## A new Cost Function\n", "\n", @@ -2581,9 +2223,7 @@ { "cell_type": "markdown", "id": "144c5c88", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta}=-\\log{\\prod_{i=0}^{n-1}p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta})}=-\\sum_{i=0}^{n-1}\\log{p(y_i,\\boldsymbol{X}\\vert\\boldsymbol{\\beta})},\n", @@ -2593,9 +2233,7 @@ { "cell_type": "markdown", "id": "76c7c5c6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which becomes" ] @@ -2603,9 +2241,7 @@ { "cell_type": "markdown", "id": "f78405b0", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta}=\\frac{n}{2}\\log{2\\pi\\sigma^2}+\\frac{\\vert\\vert (\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\vert\\vert_2^2}{2\\sigma^2}.\n", @@ -2615,9 +2251,7 @@ { "cell_type": "markdown", "id": "dd7b69b5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Taking the derivative of the *new* cost function with respect to the parameters $\\beta$ we recognize our familiar OLS equation, namely" ] @@ -2625,9 +2259,7 @@ { "cell_type": "markdown", "id": "a1430a20", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{X}^T\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right) =0,\n", @@ -2637,9 +2269,7 @@ { "cell_type": "markdown", "id": "dfcd0c7b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which leads to the well-known OLS equation for the optimal paramters $\\beta$" ] @@ -2647,9 +2277,7 @@ { "cell_type": "markdown", "id": "043d37a6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}}^{\\mathrm{OLS}}=\\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}!\n", @@ -2659,9 +2287,7 @@ { "cell_type": "markdown", "id": "28fe8edc", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Before we make a similar analysis for Ridge and Lasso regression, we need a short reminder on statistics." ] @@ -2669,9 +2295,7 @@ { "cell_type": "markdown", "id": "b600d21c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## More basic Statistics and Bayes' theorem\n", "\n", @@ -2689,9 +2313,7 @@ { "cell_type": "markdown", "id": "d513cf2c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X \\cup Y)= p(X)+p(Y)-p(X \\cap Y).\n", @@ -2701,9 +2323,7 @@ { "cell_type": "markdown", "id": "d18008b5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "**The product rule (aka joint probability) is given by.**" ] @@ -2711,9 +2331,7 @@ { "cell_type": "markdown", "id": "22f55ec4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X \\cup Y)= p(X,Y)= p(X\\vert Y)p(Y)=p(Y\\vert X)p(X),\n", @@ -2723,9 +2341,7 @@ { "cell_type": "markdown", "id": "29e1bff6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "where we read $p(X\\vert Y)$ as the likelihood of obtaining $X$ given $Y$.\n", "\n", @@ -2735,9 +2351,7 @@ { "cell_type": "markdown", "id": "333869a5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Marginal Probability\n", "\n", @@ -2747,9 +2361,7 @@ { "cell_type": "markdown", "id": "bf8de993", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X)=\\sum_{i=0}^{n-1}p(X,Y=y_i)=\\sum_{i=0}^{n-1}p(X\\vert Y=y_i)p(Y=y_i)=\\sum_{i=0}^{n-1}p(X\\vert y_i)p(y_i).\n", @@ -2759,9 +2371,7 @@ { "cell_type": "markdown", "id": "1c5349d1", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Conditional Probability\n", "\n", @@ -2771,9 +2381,7 @@ { "cell_type": "markdown", "id": "dbede038", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X\\vert Y)= \\frac{p(X,Y)}{p(Y)}=\\frac{p(X,Y)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)}.\n", @@ -2783,9 +2391,7 @@ { "cell_type": "markdown", "id": "c1334a9c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Bayes' Theorem\n", "\n", @@ -2795,9 +2401,7 @@ { "cell_type": "markdown", "id": "ac8dcf04", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X\\vert Y)= \\frac{p(X,Y)}{p(Y)},\n", @@ -2807,9 +2411,7 @@ { "cell_type": "markdown", "id": "1992e354", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which we can rewrite as" ] @@ -2817,9 +2419,7 @@ { "cell_type": "markdown", "id": "24e5179f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X\\vert Y)= \\frac{p(X,Y)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)}=\\frac{p(Y\\vert X)p(X)}{\\sum_{i=0}^{n-1}p(Y\\vert X=x_i)p(x_i)},\n", @@ -2829,9 +2429,7 @@ { "cell_type": "markdown", "id": "42979e5a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which is Bayes' theorem. It allows us to evaluate the uncertainty in in $X$ after we have observed $Y$. We can easily interchange $X$ with $Y$." ] @@ -2839,9 +2437,7 @@ { "cell_type": "markdown", "id": "86a8cd07", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Interpretations of Bayes' Theorem\n", "\n", @@ -2858,9 +2454,7 @@ { "cell_type": "markdown", "id": "b9a4f9a9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Example of Usage of Bayes' theorem\n", "\n", @@ -2879,9 +2473,7 @@ { "cell_type": "markdown", "id": "888c64ec", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X=1\\vert Y=1) =0.8.\n", @@ -2891,9 +2483,7 @@ { "cell_type": "markdown", "id": "5086b046", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of $80\\%$ for having cancer.\n", "It is however not correct, as the following Bayesian analysis shows." @@ -2902,9 +2492,7 @@ { "cell_type": "markdown", "id": "609ea2b0", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Doing it correctly\n", "\n", @@ -2915,9 +2503,7 @@ { "cell_type": "markdown", "id": "687bfdf3", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(Y=1) =0.004.\n", @@ -2927,9 +2513,7 @@ { "cell_type": "markdown", "id": "c2af1d2d", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have" ] @@ -2937,9 +2521,7 @@ { "cell_type": "markdown", "id": "b8ec8d93", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(X=1\\vert Y=0) =0.1.\n", @@ -2949,9 +2531,7 @@ { "cell_type": "markdown", "id": "48a6fb5b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute" ] @@ -2959,9 +2539,7 @@ { "cell_type": "markdown", "id": "ba37d4fd", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(Y=1\\vert X=1)=\\frac{p(X=1\\vert Y=1)p(Y=1)}{p(X=1\\vert Y=1)p(Y=1)+p(X=1\\vert Y=0)p(Y=0)}=\\frac{0.8\\times 0.004}{0.8\\times 0.004+0.1\\times 0.996}=0.031.\n", @@ -2971,9 +2549,7 @@ { "cell_type": "markdown", "id": "824f1789", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "That is, in case of a positive test, there is only a $3\\%$ chance of having breast cancer!" ] @@ -2981,9 +2557,7 @@ { "cell_type": "markdown", "id": "07f89fd5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Bayes' Theorem and Ridge and Lasso Regression\n", "\n", @@ -2998,9 +2572,7 @@ { "cell_type": "markdown", "id": "8fcd73d9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Test Function for what happens with OLS, Ridge and Lasso\n", "\n", @@ -3017,10 +2589,7 @@ "cell_type": "code", "execution_count": 7, "id": "4402e113", - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -3091,9 +2660,7 @@ { "cell_type": "markdown", "id": "b700e45e", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "How can we understand this?" ] @@ -3101,9 +2668,7 @@ { "cell_type": "markdown", "id": "69659d9d", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Invoking Bayes' theorem\n", "\n", @@ -3115,9 +2680,7 @@ { "cell_type": "markdown", "id": "80b423eb", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "\\boldsymbol{D}=[(x_0,y_0), (x_1,y_1),\\dots, (x_{n-1},y_{n-1})],\n", @@ -3127,9 +2690,7 @@ { "cell_type": "markdown", "id": "de547949", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "is given by" ] @@ -3137,9 +2698,7 @@ { "cell_type": "markdown", "id": "99d9476a", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{D}\\vert\\boldsymbol{\\beta})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}.\n", @@ -3149,9 +2708,7 @@ { "cell_type": "markdown", "id": "aa5613ea", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "In Bayes' theorem this function plays the role of the so-called likelihood. We could now ask the question what is the posterior probability of a parameter set $\\boldsymbol{\\beta}$ given a domain of events $\\boldsymbol{D}$? That is, how can we define the posterior probability" ] @@ -3159,9 +2716,7 @@ { "cell_type": "markdown", "id": "7dea4e11", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{\\beta}\\vert\\boldsymbol{D}).\n", @@ -3171,9 +2726,7 @@ { "cell_type": "markdown", "id": "99ab1d57", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Bayes' theorem comes to our rescue here since (omitting the normalization constant)" ] @@ -3181,9 +2734,7 @@ { "cell_type": "markdown", "id": "170ee68b", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{\\beta}\\vert\\boldsymbol{D})\\propto p(\\boldsymbol{D}\\vert\\boldsymbol{\\beta})p(\\boldsymbol{\\beta}).\n", @@ -3193,9 +2744,7 @@ { "cell_type": "markdown", "id": "3388d6a4", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We have a model for $p(\\boldsymbol{D}\\vert\\boldsymbol{\\beta})$ but need one for the **prior** $p(\\boldsymbol{\\beta}$!" ] @@ -3203,9 +2752,7 @@ { "cell_type": "markdown", "id": "a5bdaa16", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Ridge and Bayes\n", "\n", @@ -3219,9 +2766,7 @@ { "cell_type": "markdown", "id": "cbf8e7e6", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{\\beta})=\\prod_{j=0}^{p-1}\\exp{\\left(-\\frac{\\beta_j^2}{2\\tau^2}\\right)}.\n", @@ -3231,9 +2776,7 @@ { "cell_type": "markdown", "id": "52c9e48f", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Our posterior probability becomes then (omitting the normalization factor which is just a constant)" ] @@ -3241,9 +2784,7 @@ { "cell_type": "markdown", "id": "1ae37a89", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{\\beta\\vert\\boldsymbol{D})}=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}\\prod_{j=0}^{p-1}\\exp{\\left(-\\frac{\\beta_j^2}{2\\tau^2}\\right)}.\n", @@ -3253,9 +2794,7 @@ { "cell_type": "markdown", "id": "dd8fee27", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "We can now optimize this quantity with respect to $\\boldsymbol{\\beta}$. As we\n", "did for OLS, this is most conveniently done by taking the negative\n", @@ -3266,9 +2805,7 @@ { "cell_type": "markdown", "id": "de603b44", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta})=\\frac{\\vert\\vert (\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\vert\\vert_2^2}{2\\sigma^2}+\\frac{1}{2\\tau^2}\\vert\\vert\\boldsymbol{\\beta}\\vert\\vert_2^2,\n", @@ -3278,9 +2815,7 @@ { "cell_type": "markdown", "id": "3f08e4a0", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and replacing $1/2\\tau^2$ with $\\lambda$ we have" ] @@ -3288,9 +2823,7 @@ { "cell_type": "markdown", "id": "b65ee2f0", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta})=\\frac{\\vert\\vert (\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\vert\\vert_2^2}{2\\sigma^2}+\\lambda\\vert\\vert\\boldsymbol{\\beta}\\vert\\vert_2^2,\n", @@ -3300,9 +2833,7 @@ { "cell_type": "markdown", "id": "6526b426", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which is our Ridge cost function! Nice, isn't it?" ] @@ -3310,9 +2841,7 @@ { "cell_type": "markdown", "id": "29190fe5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "## Lasso and Bayes\n", "\n", @@ -3322,9 +2851,7 @@ { "cell_type": "markdown", "id": "71b1585c", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{\\beta})=\\prod_{j=0}^{p-1}\\exp{\\left(-\\frac{\\vert\\beta_j\\vert}{\\tau}\\right)}.\n", @@ -3334,9 +2861,7 @@ { "cell_type": "markdown", "id": "e18baaa9", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Our posterior probability becomes then (omitting the normalization factor which is just a constant)" ] @@ -3344,9 +2869,7 @@ { "cell_type": "markdown", "id": "1c58b8db", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "p(\\boldsymbol{\\beta}\\vert\\boldsymbol{D})=\\prod_{i=0}^{n-1}\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left[-\\frac{(y_i-\\boldsymbol{X}_{i,*}\\boldsymbol{\\beta})^2}{2\\sigma^2}\\right]}\\prod_{j=0}^{p-1}\\exp{\\left(-\\frac{\\vert\\beta_j\\vert}{\\tau}\\right)}.\n", @@ -3356,9 +2879,7 @@ { "cell_type": "markdown", "id": "238ee8c5", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "Taking the negative\n", "logarithm of the posterior probability and leaving out the\n", @@ -3368,9 +2889,7 @@ { "cell_type": "markdown", "id": "46389b82", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta}=\\frac{\\vert\\vert (\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\vert\\vert_2^2}{2\\sigma^2}+\\frac{1}{\\tau}\\vert\\vert\\boldsymbol{\\beta}\\vert\\vert_1,\n", @@ -3380,9 +2899,7 @@ { "cell_type": "markdown", "id": "0ecbb2d2", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "and replacing $1/\\tau$ with $\\lambda$ we have" ] @@ -3390,9 +2907,7 @@ { "cell_type": "markdown", "id": "27368dfb", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "$$\n", "C(\\boldsymbol{\\beta}=\\frac{\\vert\\vert (\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta})\\vert\\vert_2^2}{2\\sigma^2}+\\lambda\\vert\\vert\\boldsymbol{\\beta}\\vert\\vert_1,\n", @@ -3402,15 +2917,31 @@ { "cell_type": "markdown", "id": "27082f48", - "metadata": { - "editable": true - }, + "metadata": {}, "source": [ "which is our Lasso cost function!" ] } ], - "metadata": {}, + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.9.10" + } + }, "nbformat": 4, "nbformat_minor": 5 }