diff --git a/doc/LectureNotes/exercisesweek35.ipynb b/doc/LectureNotes/exercisesweek35.ipynb index 79e8f035b..7736eafb5 100644 --- a/doc/LectureNotes/exercisesweek35.ipynb +++ b/doc/LectureNotes/exercisesweek35.ipynb @@ -1168,7 +1168,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.7" + "version": "3.13.3" } }, "nbformat": 4, diff --git a/doc/LectureNotes/exercisesweek36.ipynb b/doc/LectureNotes/exercisesweek36.ipynb index 816ccabf3..4c43c5de8 100644 --- a/doc/LectureNotes/exercisesweek36.ipynb +++ b/doc/LectureNotes/exercisesweek36.ipynb @@ -709,7 +709,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.3" + "version": "3.13.7" } }, "nbformat": 4, diff --git a/doc/LectureNotes/exercisesweek38.ipynb b/doc/LectureNotes/exercisesweek38.ipynb index c100028a5..48df7c2b8 100644 --- a/doc/LectureNotes/exercisesweek38.ipynb +++ b/doc/LectureNotes/exercisesweek38.ipynb @@ -117,6 +117,20 @@ "$$\n" ] }, + { + "cell_type": "markdown", + "id": "3a6e76d0", + "metadata": {}, + "source": [ + "
\n", + "\n", + "$$\n", + "\\mathbb{E}(\\boldsymbol{\\hat{\\beta}_{OLS}}) = \\mathbb{E} \\left((X^T X)^{-1} X^T y\\right) = \\mathbb E \\left((X^T X)^{-1} X^T \\tilde y + (X^T X)^{-1} X^T \\epsilon \\right) = \\beta + \\mathbb{E} \\left((X^T X)^{-1} X^T \\epsilon \\right) = \\beta\n", + "$$\n", + "\n", + "
" + ] + }, { "cell_type": "markdown", "id": "46e93394", @@ -135,6 +149,22 @@ "$$\n" ] }, + { + "cell_type": "markdown", + "id": "e39cd224", + "metadata": {}, + "source": [ + "
\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\mathbf{Var}(\\hat \\beta_{OLS}) &= \\mathbb E (\\beta_{OLS} - \\mathbb{E}(\\hat{\\beta}_{OLS}))^2 \\\\ &= \\mathbb E (\\beta_{OLS} - \\beta)^2\\\\ &= \\mathbb E (\\beta + (X^T X)^{-1} X^T \\epsilon - \\beta)^2\\\\ &= \\mathbb E ((X^T X)^{-1} X^T \\epsilon)^2\\\\ &= \\sigma^2 \\mathbb E ((X^T X)^{-1} X^T X ((X^T X)^{-1})^T)\\\\ &= \\sigma^2 ((X^T X)^{-1})^T\\\\ &= \\sigma^2 (X^T X)^{-1}\n", + "\\end{aligned}\n", + "$$\n", + "\n", + "
" + ] + }, { "cell_type": "markdown", "id": "d2143684", @@ -166,7 +196,7 @@ "metadata": {}, "source": [ "$$\n", - "\\mathbb{E} \\big[ \\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}} \\big]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\boldsymbol{\\beta}\n", + "\\mathbb{E} [ \\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}} ]=(\\mathbf{X}^{T} \\mathbf{X} + \\lambda \\mathbf{I}_{pp})^{-1} (\\mathbf{X}^{\\top} \\mathbf{X})\\boldsymbol{\\beta}\n", "$$\n" ] }, @@ -175,7 +205,42 @@ "id": "028209a1", "metadata": {}, "source": [ - "We see that $\\mathbb{E} \\big[ \\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}} \\big] \\not= \\mathbb{E} \\big[\\hat{\\boldsymbol{\\beta}}^{\\mathrm{OLS}}\\big ]$ for any $\\lambda > 0$.\n" + "We see that $\\mathbb{E} [ \\hat{\\boldsymbol{\\beta}}^{\\mathrm{Ridge}} ] \\not= \\mathbb{E} [\\hat{\\boldsymbol{\\beta}}^{\\mathrm{OLS}} ]$ for any $\\lambda > 0$.\n" + ] + }, + { + "cell_type": "markdown", + "id": "73186528", + "metadata": {}, + "source": [ + "
\n", + "\n", + "$$\n", + "\\mathbb E \\hat \\beta_\\mathrm{Ridge} = \\mathbb E \\left((X^TX + \\lambda I)^{-1} X^T \\tilde y\\right) + \\mathbb E \\left((X^T X + \\lambda I)^{-1} X^T \\epsilon\\right) = \\mathbb E \\left((X^TX + \\lambda I)^{-1} X^T \\tilde y\\right) = (X^TX + \\lambda I)^{-1} X^T X \\beta\n", + "$$\n", + "\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "65f6f914", + "metadata": {}, + "source": [ + "**b)** Why do we say that Ridge regression gives a biased estimate? Is this a problem?\n" + ] + }, + { + "cell_type": "markdown", + "id": "241e8533", + "metadata": {}, + "source": [ + "
\n", + "\n", + "The expectation value for the ridge parameters only approaches the true value $\\beta$ in the limit $\\lambda \\to 0$, this is why we call it biased. As independent of number of training samples our estimate will always differ from the true value given $\\lambda \\neq 0$ (in which case it would be OLS). This may be a problem if the goal of our analysis is the closest possible parameter estimates given near infinite number of training samples (because in this case the deviation from the true value approaches 0).\n", + "\n", + "
" ] }, { @@ -204,6 +269,26 @@ "We see that if the parameter $\\lambda$ goes to infinity then the variance of the Ridge parameters $\\boldsymbol{\\beta}$ goes to zero.\n" ] }, + { + "cell_type": "markdown", + "id": "06835e6f", + "metadata": {}, + "source": [ + "
\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\mathrm{Var}(\\hat \\beta_\\mathrm{Ridge}) &= \\mathbb E (\\hat \\beta_\\mathrm{Ridge} - \\mathbb E \\hat \\beta_\\mathrm{Ridge})^2\\\\ \n", + "&= \\mathbb E \\left[\\left((X^TX + \\lambda I)^{-1} X^T \\tilde y\\right) + \\left((X^T X + \\lambda I)^{-1} X^T \\epsilon\\right) - (X^TX + \\lambda I)^{-1} X^T X \\beta \\right]^2 \\\\\n", + "&= \\mathbb E \\left[(X^T X + \\lambda I)^{-1} X^T \\epsilon\\right]^2\\\\\n", + "&= \\sigma^2 (X^T X + \\lambda I)^{-1} X^T X ((X^T X + \\lambda I)^{-1})^T\n", + "\\end{aligned}\n", + "$$\n", + "\n", + "\n", + "
" + ] + }, { "cell_type": "markdown", "id": "74bc300b", @@ -294,6 +379,30 @@ "In order to arrive at the equation for the bias, we have to approximate the unknown function $f$ with the output/target values $y$.\n" ] }, + { + "cell_type": "markdown", + "id": "7dcfe676", + "metadata": {}, + "source": [ + "
\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\mathbb E[(y-\\tilde y)^2]\n", + "&= \\mathbb E[(f+\\varepsilon-\\tilde y)^2] \\\\\n", + "&= \\mathbb E[(f - \\mathbb E[\\tilde y] + \\mathbb E[\\tilde y]-\\tilde y + \\varepsilon)^2] \\\\\n", + "&= \\mathbb E[(f - \\mathbb E[\\tilde y])^2] + \\mathbb E[(\\tilde y - \\mathbb E[\\tilde y])^2] + \\mathbb E[\\varepsilon^2] \\\\\n", + "&\\quad + 2\\,\\mathbb E[(f - \\mathbb E[\\tilde y])(\\mathbb E[\\tilde y]-\\tilde y)]\n", + "+ 2\\,\\mathbb E[(f - \\mathbb E[\\tilde y])\\varepsilon]\n", + "+ 2\\,\\mathbb E[(\\mathbb E[\\tilde y]-\\tilde y)\\varepsilon] \\\\\n", + "&= (f - \\mathbb E[\\tilde y])^2 + \\mathrm{Var}(\\tilde y) + \\sigma^2,\n", + "\\end{aligned}\n", + "$$\n", + "\n", + "\n", + "
" + ] + }, { "cell_type": "markdown", "id": "70fbfcd7", @@ -302,6 +411,23 @@ "**b)** Explain what the terms mean and discuss their interpretations.\n" ] }, + { + "cell_type": "markdown", + "id": "6f47d057", + "metadata": {}, + "source": [ + "
\n", + "\n", + "1. Bias: Difference between true function and average prediction. Systematic error. High bias = underfitting.\n", + "\n", + "2. Variance: How much predictions change if we retrain on different data. Sensitivity to training data. High variance = overfitting.\n", + "\n", + "3. Irreducible noise $\\sigma^2$: Random noise in data we can’t model away. Sets the minimum possible error.\n", + "\n", + "\n", + "
" + ] + }, { "cell_type": "markdown", "id": "b8f8b9d1", @@ -325,7 +451,15 @@ "execution_count": null, "id": "b5bf581c", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "MSE (218.336) = Bias^2 (209.927) + Variance (8.331) = 218.258\n" + ] + } + ], "source": [ "import numpy as np\n", "\n", @@ -336,9 +470,30 @@ "# The definition of targets has been updated, and was wrong earlier in the week.\n", "targets = np.random.rand(1, n)\n", "\n", - "mse = ...\n", - "bias = ...\n", - "variance = ..." + "def calculate_key_metrics(predictions, targets, noise_free_targets=None):\n", + " y_mean = predictions.mean(axis=0)\n", + "\n", + " # --- MSE (with noisy targets) ---\n", + " mse = np.mean((predictions - targets) ** 2)\n", + "\n", + " # --- Bias^2 (wrt noise-free targets) ---\n", + " if noise_free_targets is None:\n", + " targets_nf = targets\n", + " else:\n", + " targets_nf = noise_free_targets\n", + " bias = np.mean((y_mean - targets_nf.mean(axis=0)) ** 2)\n", + "\n", + " # --- Variance (spread of predictions around their mean) ---\n", + " variance = np.mean((predictions - y_mean) ** 2)\n", + "\n", + " return mse, bias, variance\n", + "\n", + "def print_key_metrics(predictions, targets):\n", + " mse, bias, variance = calculate_key_metrics(predictions, targets)\n", + "\n", + " print(f\"MSE ({mse:.3f}) = Bias ({bias:.3f}) + Variance ({variance:.3f}) = {bias + variance:.3f}\")\n", + "\n", + "print_key_metrics(predictions, targets)" ] }, { @@ -346,11 +501,57 @@ "id": "7b1dc621", "metadata": {}, "source": [ - "**b)** Change the prediction values in some way to increase the bias while decreasing the variance.\n", + "**b)** Change the prediction values in some way to increase the bias while decreasing the variance.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "id": "c8e777a6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "MSE (731.057) = Bias^2 (728.899) + Variance (2.076) = 730.975\n" + ] + } + ], + "source": [ + "predictions = predictions * 0.5 + 20\n", + "\n", + "print_key_metrics(predictions, targets)" + ] + }, + { + "cell_type": "markdown", + "id": "dbce0baf", + "metadata": {}, + "source": [ "\n", "**c)** Change the prediction values in some way to increase the variance while decreasing the bias.\n" ] }, + { + "cell_type": "code", + "execution_count": 71, + "id": "a30ea2b6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "MSE (830.804) = Bias^2 (0.252) + Variance (830.500) = 830.751\n" + ] + } + ], + "source": [ + "predictions = (predictions - np.mean(predictions, axis=0)) * 20\n", + "print_key_metrics(predictions, targets)" + ] + }, { "cell_type": "markdown", "id": "8da63362", @@ -361,7 +562,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 72, "id": "dd5855e4", "metadata": {}, "outputs": [], @@ -379,39 +580,83 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 81, "id": "7e35fa37", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAW41JREFUeJzt3Xd4VGXCNvB7ZlImddIrqSQ0pUNiUIoSiAgKr4WyKlXsCi+6NBHEBigL6MIH6q5gY0H2FXbXRRRRioCANJGeRgKkQjKTSZ3MPN8fJIcckgnpk+Tcv+vKBXPmnDPPOZlk7jxVJYQQICIiIlIQta0LQERERNTSGICIiIhIcRiAiIiISHEYgIiIiEhxGICIiIhIcRiAiIiISHEYgIiIiEhx7GxdgNbIYrHg6tWrcHNzg0qlsnVxiIiIqA6EECgoKEBQUBDU6trreBiAanD16lWEhITYuhhERETUAOnp6ejQoUOt+zAA1cDNzQ2ouIHu7u62Lg4RERHVgcFgQEhIiPQ5XhsGoBpUNnu5u7szABEREbUxdem+wk7QREREpDgMQERERKQ4DEBERESkOOwD1Ahmsxkmk8nWxaDbsLe3h0ajsXUxiIioFWEAagAhBDIzM5Gfn2/rolAdeXh4ICAggPM6ERERwADUMJXhx8/PD87OzvxQbcWEECgqKkJ2djYAIDAw0NZFIiKiVoABqJ7MZrMUfry9vW1dHKoDJycnAEB2djb8/PzYHEZEROwEXV+VfX6cnZ1tXRSqh8rvF/tsERERGIAajs1ebQu/X0REVBUDEBERESkOAxAREREpDgMQERERKQ4DkIJMnjwZKpUKzz77bLXnXnjhBahUKkyePBkAkJOTg+eeew6hoaFwdHREQEAAEhISsH//fumY8PBwqFSqal9Lly5t0esiIqK2QwiB5Pxk5Bbn2rQcHAavMCEhIdi0aRNWrlwpDQ8vKSnBxo0bERoaKu33yCOPoKysDJ999hkiIyORlZWFXbt24dq1a7Lzvfnmm5g+fbpsm5ubWwtdDRERtXZCCKQXpONw5mEczjyMI5lHkFuci1f7vYpJd0yyWbkYgJqAEALFJnOLv66Tvabeo5v69OmDpKQkfPPNN3j88ccBAN988w1CQ0MREREBAMjPz8e+ffuwe/duDB48GAAQFhaGmJiYaudzc3NDQEBAk1wPERG1DxnGDBzKPIQjmUdwOPMwMgszZc87ahyhL9XbrHxgAGoaxSYzui38vsVf98ybCXB2qP+3cOrUqVi/fr0UgD799FNMmTIFu3fvBgC4urrC1dUV27Ztw1133QVHR8cmLzsREbUfOUU5Uu3OoYxDuGy8LHveTm2HHj49EBMYg5iAGPTw7QFHjW0/WxiAFOiJJ57AvHnzcOnSJQDA/v37sWnTJikA2dnZYcOGDZg+fTrWrVuHPn36YPDgwRg/fjx69OghO9ecOXOwYMEC2bbvvvsOAwcObMErIiKilnS95DqOZB6RanhS9Cmy5zUqDe7wvgP9A/ojJjAGvXx7wdm+dU0gzADUBJzsNTjzZoJNXrchfH19MXLkSGzYsAFCCIwcORI+Pj6yfR555BGMHDkS+/btw6+//orvvvsO7733Hv72t79JHaUB4M9//rPsMQAEBwc38IqIiKg10pfqcTTrqNSP52LeRdnzKqjQxasLYgJiEBMYgz5+feDq4Gqz8tYFA1ATUKlUDWqKsqWpU6fixRdfBACsWbOmxn20Wi2GDRuGYcOG4fXXX8dTTz2FRYsWyQKPj48PoqKiWqzcRETU/ApNhTiadVRq0jp3/RwEhGyfKI8oxAbGon9Af/Tz7wedo85m5W2ItvWpTU3m/vvvR1lZGVQqFRIS6lZ71a1bN2zbtq3Zy0ZERC2ruLwYJ7JPSDU8p3NPwyzkg3vC3cOlGp5+/v3g7dS2FwRnAFIojUaDs2fPSv+v6tq1a3jssccwdepU9OjRA25ubvjtt9/w3nvvYfTo0bJ9CwoKkJkp793v7OwMd3f3FrgKIiJqiDJzGU7mnJRqeE7lnoLJIl8sOtg1WKrhiQmIgZ+zn83K2xwYgBTMWkhxdXVFbGwsVq5ciaSkJJhMJoSEhGD69OmYP3++bN+FCxdi4cKFsm3PPPMM1q1b16xlJyKiujNZTDide1qq4TmRfQKl5lLZPv7O/lINT0xADIJcg2xW3pagEkKIOuynKAaDATqdDnq9vlpIKCkpQUpKCiIiIqDVam1WRqofft+ISEnMFjPOXT+Hw5mHcSjzEI5lHUNxebFsH2+tN2ICYtA/sD9iA2IR4hZS77nlWpvaPr9vxRogIiKiNs4iLLiYd1Gq4TmaeRQFpgLZPjpH3Y3AU9GkFamLbPOBpzEYgIiIiNoYIQRS9Cmy5SXyS/Nl+7jau6Kffz+pSSvaMxpqFZcArcQARERE1MpZW0+rKic7J/Tx73OjH09ADLp4dYGdmh/z1vDOEBERtUJXjVelsGNtPa1evr2kGp47fO6AvdreZuVtaxiAiIiIWoHK9bQOZx7G4YzDbWI9rbaMAYiIiMgG6rqeVkzgjY7Lvf16w8nOyWblbW8YgIiIiFpAe1xPqy1jACIiImoGdVlPK9ozWhqa3hbX02rLGIBIkpqaioiICBw/fhy9evVq0nMnJibi7rvvRkFBAXbu3Im77767Sc9PRGRrxeXFOJ59XGrSsraeVuXyEv0D+sNL62Wz8iodA5CCTJ48GZ999pn02MvLC/3798d7772HHj16ICQkBBkZGfDx8WnS17169SqGDRuGe+65B0FBQRg1ahT27t2L7t27S/uYTCYsWLAA27dvR3JyMnQ6HeLj47F06VIEBbXv6diJqG26dT2t33N/R7mlXLZPB9cOUqfl/gH92916Wm0ZA5DC3H///Vi/fj0AIDMzEwsWLMCoUaOQlpYGjUaDgICAJn29vLw8JCQkYODAgVi/fj00Gg1cXV2RkJCA/fv3IyIiAgBQVFSEY8eO4fXXX0fPnj2Rl5eHGTNm4KGHHsJvv/3WpGUiImqIuq6nVXUB0fa+nlZbxgCkMI6OjlLICQgIwNy5czFw4EDk5OSgsLBQ1gRmNpvx9NNP46effkJmZiZCQ0Px/PPPY8aMGdL5du/ejdmzZ+P06dOwt7fHHXfcgY0bNyIsLAxFRUUYOXIk7r77bqxdu1aacn3JkiVwdXXF8OHDsX//fvj5+UGn02Hnzp2ysq5evRoxMTFIS0tDaGhoC98pIlK6yvW0DmUewuHMw7Wup1VZy9Me1tNSCgagpiAEYCpq+de1dwYa8YNmNBrx5ZdfIioqCt7e3igsLJQ9b7FY0KFDB2zZsgXe3t44cOAAnn76aQQGBmLs2LEoLy/HmDFjMH36dPzjH/9AWVkZDh8+LP3wOzs748CBAzW+9muvvYbXXnut1vLp9XqoVCp4eHg0+BqJiOqqLutpeTh6SP13YgNiEaGLYOBpoxiAmoKpCHjXBtWc868CDi71OuTbb7+Fq+uNYZWFhYUIDAzEt99+C7W6+vow9vb2WLx4sfQ4IiICBw8exNdff42xY8fCYDBAr9dj1KhR6NixIwCga9eujb4sVKzePmfOHEyYMOG2K/oSETVE5XpahzIPSfPx3Lqelpu9G/oG9JWWl+B6Wu0HA5DC3HvvvVi7di1Q0T/n//2//4cRI0bg8OHDNe6/Zs0afPrpp0hLS0NxcTHKysqkEWJeXl6YPHkyEhISMGzYMMTHx2Ps2LEIDAxsVBlNJhPGjh0LIYRUViKixpKtp5VxGEeyal9PKzYgFl28ukCj1tiszNR8GICagr3zjdoYW7xuPbm4uCAqKkp6/Le//Q06nQ6ffPIJnnrqKdm+mzZtwquvvoq//OUviIuLg5ubG95//30cOnRI2mf9+vV4+eWXsWPHDmzevBkLFizAzp07cddddzXokirDz6VLl/DTTz+x9oeIGqXqelqHMg4hqyhL9ryjxhG9/HpJNTxcT0s5GICagkpV76ao1kKlUkGtVqO4uLjac/v378eAAQPw/PPPS9uSkpKq7de7d2/07t0b8+bNQ1xcHDZu3NigAFQZfi5evIiff/4Z3t7eDbgiIlKyuq6nVTlSi+tpKVerCEBr1qzB+++/j8zMTPTs2RN//etfERMTU+O+n3zyCT7//HP88ccfAIC+ffvi3Xffle1/63w3AJCQkIAdO3Y085W0fqWlpcjMvLGicF5eHlavXg2j0YgHH3yw2r7R0dH4/PPP8f333yMiIgJffPEFjhw5Ig1dT0lJwccff4yHHnoIQUFBOH/+PC5evIiJEyfWu1wmkwmPPvoojh07hm+//RZms1kqp5eXFxwcHBp97UTU/lRdT+tQxiGkGlJlz2tUGtzhc4dUw9PLrxfX0yKgNQSgzZs3Y9asWVi3bh1iY2OxatUqJCQk4Pz58/Dzqz5h1O7duzFhwgQMGDAAWq0Wy5Ytw/Dhw3H69GkEBwdL+1Wd7wYVw78J2LFjh9RHx83NDV26dMGWLVswZMgQpKbKf3E888wzOH78OMaNGweVSoUJEybg+eefx3fffQdUjPI6d+4cPvvsM1y7dg2BgYF44YUX8Mwzz9S7XFeuXMG///1vAKg2C/XPP/+MIUOGNOKqiai90Jfq8VvWb9Jsy9bW06qs4enr3xcu9m2zhp6al0oIIeqwX7OJjY1F//79sXr1aqBi6HVISAheeuklzJ0797bHm81meHp6YvXq1VLNw+TJk5Gfn49t27bVqQylpaUoLb05mZXBYEBISAj0en21PiglJSVISUlBREQEtFptPa+WbIXfN6K2qXI9rcMZN5q1altPKyYgBn39+3I9LQUzGAzQ6XQ1fn7fyqY1QGVlZTh69CjmzZsnbVOr1YiPj8fBgwfrdI6ioiKYTCZ4ecnXU9m9ezf8/Pzg6emJ++67D2+//bbVPiVLliyRDfcmIiLb4Hpa1FJsGoByc3NhNpvh7+8v2+7v749z587V6Rxz5sxBUFAQ4uPjpW33338/Hn74YURERCApKQnz58/HiBEjcPDgQWg01Yczzps3D7NmzZIeV9YAERFR86pcT6uy0zLX06KWYvM+QI2xdOlSbNq0Cbt375Y1a4wfP176f/fu3dGjRw907NgRu3fvxtChQ6udx9HRkX2EiIhagGw9rYzDOJFTfT2tAJcAKexwPS1qLjYNQD4+PtBoNMjKks/LkJWVddtFOZcvX46lS5fixx9/RI8ePWrdNzIyEj4+PkhMTKwxABERUfOo83paFTU8XE+LWopNA5CDgwP69u2LXbt2YcyYMUBFJ+hdu3bhxRdftHrce++9h3feeQfff/89+vXrd9vXuXz5sjRKiYiImo9sPa2MwziaxfW0qHWyeRPYrFmzMGnSJPTr1w8xMTFYtWoVCgsLMWXKFADAxIkTERwcjCVLlgAAli1bhoULF2Ljxo0IDw+X5opxdXWFq6srjEYjFi9ejEceeQQBAQFISkrC7NmzERUVhYSEBJteKxFReyOEQLI+WZptmetpUVth8wA0btw45OTkYOHChcjMzESvXr2wY8cOqWN0WlqabKHOtWvXoqysDI8++qjsPIsWLcIbb7wBjUaD33//HZ999hny8/MRFBSE4cOH46233mI/HyKiRqpcT+tQ5iEcybgxUutayTXZPpXracUGxCImIIbraVGrZPN5gFqj2uYR4HwybRO/b0QNx/W0qK1oM/MAERFR61FkKkKKIQWp+lSkGlKRqk/FH7l/3HY9rZ6+PeGg4XI11LYwAFGjqVQqbN26VerITkStl0VYkFGYIYWcFP2NwJNiSEF2UXaNx3A9LWqPGIAU4sEHH4TJZKpxQdh9+/Zh0KBBOHny5G2nFKhJRkYGPD09m6ikRNQUjGXGmwGnojYnxZCCNENatXl3qvLSeiHcPRwRugiEu4cjyjMKvf16cz0tancYgBRi2rRpeOSRR3D58mV06NBB9tz69evRr1+/eoefsrIyODg43HbOJiJqHmaLGVcLr8qarCqbsHKKc6weZ6+2R6hbKMJ1N4NOuC4c4e7hXEeLFIMBSCFGjRoFX19fbNiwAQsWLJC2G41GbNmyBXPnzsWECROwd+9e5OXloWPHjpg/fz4mTJgg7TtkyBDceeedsLOzw5dffonu3bvj559/rtYENmfOHGzduhWXL19GQEAAHn/8cSxcuBD29jc6Rb7xxhvYtm0bXnnlFbz++uvIy8vDiBEj8Mknn8DNzQ2omA9q+fLl+Pjjj5Geng5/f38888wzeO211wAA6enpeOWVV/DDDz9ArVZj4MCB+OCDDxAeHt7Cd5ao+RWUFcibrCr+TTOkocxSZvU4b633jYBTEW4idBGIcI9AoGsg7NT89U/Kxp+AJiCEqDazaUtwsnOq8+RhdnZ2mDhxIjZs2IDXXntNOm7Lli0wm8144oknsGXLFsyZMwfu7u7473//iyeffBIdO3ZETEyMdJ7PPvsMzz33HPbv32/1tdzc3LBhwwYEBQXh1KlTmD59Otzc3DB79mxpn6SkJGzbtg3ffvst8vLyMHbsWCxduhTvvPMOULE+2yeffIKVK1finnvuQUZGhrQ+nMlkQkJCAuLi4rBv3z7Y2dnh7bffxv3334/ff/8dDg7sjEltT7mlHFeNV6s3W+lTqg0zr8pB7YBQ91CpJqfy3zBdGNwdah8FQ6RkHAZfg/oOgy8yFSF2Y2yLl/PQnw7B2d65zvufO3cOXbt2xc8//4whQ4YAAAYNGoSwsDB88cUX1fYfNWoUunTpguXLlwMVNUAGgwHHjh2T7Xe7TtDLly/Hpk2b8NtvvwEVNUDvv/8+MjMzpRqf2bNnY+/evfj1119RUFAAX19frF69Gk899VS183355Zd4++23cfbsWSnIlZWVwcPDA9u2bcPw4cOrHcNh8NRa6Ev1snBT+f+0gjSYLCarx/k6+d5osnK/WaMTrgtHkEsQ59ghqsBh8FSjLl26YMCAAfj0008xZMgQJCYmYt++fXjzzTdhNpvx7rvv4uuvv8aVK1dQVlaG0tJSODvLA1bfvn1v+zqbN2/Ghx9+iKSkJBiNRpSXl1d7I4aHh0vhBwACAwORnX1jBMrZs2dRWlpqdd22kydPIjExUXY8KkJOUlJSve4JUXMot5TjivGKNMKqaq3O9ZLrVo9z1DgizD1M1icnUheJMPcwuDq4tug1ELV3DEBNwMnOCYf+dMgmr1tf06ZNw0svvYQ1a9Zg/fr16NixIwYPHoxly5bhgw8+wKpVq9C9e3e4uLhg5syZKCuT9y9wcal9JMjBgwfx+OOPY/HixUhISIBOp8OmTZvwl7/8RbZfZX+gSiqVChaL5cZ1OdV+XUajEX379sVXX31V7TlfX9/b3gOippJfkl+tX06qIRXpBekot5RbPc7P2U+qyanaCTnQJZBLRBC1EAagJqBSqerVFGVLY8eOxYwZM7Bx40Z8/vnneO6556BSqbB//36MHj0aTzzxBFDRCfnChQvo1q1bvc5/4MABhIWFSZ2VAeDSpUv1Okd0dDScnJywa9euGpvA+vTpg82bN8PPz++2VZxEjWWymHC54LJshFVls1VeaZ7V47QaLcLcw2SdkCv/5ZByIttjAFIYV1dXjBs3DvPmzYPBYMDkyZOBitDxz3/+EwcOHICnpydWrFiBrKysegeg6OhopKWlYdOmTejfvz/++9//YuvWrfU6h1arxZw5czB79mw4ODjg7rvvRk5ODk6fPo1p06bh8ccfx/vvv4/Ro0fjzTffRIcOHXDp0iV88803mD17drVh/kR1kVeSV23OnFR9Ki4XXEa5sF6bE+AScCPcVAScypFW/i7+rM0hasUYgBRo2rRp+Pvf/44HHngAQUFBAIAFCxYgOTkZCQkJcHZ2xtNPP40xY8ZAr9fX69wPPfQQ/vd//xcvvvgiSktLMXLkSLz++ut444036nWe119/HXZ2dli4cCGuXr2KwMBAPPvsswAAZ2dn7N27F3PmzMHDDz+MgoICBAcHY+jQoawRolqZzCakF6RL4UYKPIZU6Eutv9ed7JykkFO1RifMPazN1P4SkRxHgdWAi6G2P/y+KYcQAtdLrsuWeagMOZcLLsMszFaPDXQJrDYxYIQuAv7O/nWecoKIbIejwIio3SszlyHNkCaFm6prWhWUFVg9ztnOWdYnp7LJKtQ9lOtbESkIAxARtVpCCFwruYYUfYqsf06qIRVXjFdgEZYaj1NBhSDXINnEgJWhx8/Zj7U5RMQARES2V2ouxSXDJfmaVhWBx2gyWj3O1d61WnNVuC4coW6h0NqxqZOIrGMAIqIWIYRATnGObGLAys7IV41XIVBzd0S1So0glyD5elYVtTo+Tj6szSGiBmEAaiD2HW9b+P1qOSXlJbhkuFRtzpxUQyoKTYVWj3Ozd6txdfJQ91A4ahxb9BqIqP1jAKqnyhmMi4qKbjtjMbUeRUVFQA0zUFPDCCGQVZQlCzeVnZAzCjNqrc3p4NpB3gm5YkZkb603a3OIqMUwANWTRqOBh4eHtG6Vs7Mzf2m3YkIIFBUVITs7Gx4eHtBouGhkfRSXF0t9c6o2WV0yXEJReZHV49wd3GXhpvLfELcQOGgcWvQaiIhqwgDUAAEBAQAghSBq/Tw8PKTvG8lZhAXZRdlI1ifLOyEbUpBZmGn1OI1Kgw5uHWSrk1d2QvZ09OQfBkTUqjEANYBKpUJgYCD8/PxgMplsXRy6DXt7e9b8ACgyFVVvsjLcqM0pLi+2epyHo0eNI61CXENgr2GTIhG1TQxAjaDRaPjBSq2KRViQWZgp1eBUnTsnqyjL6nF2KrsbtTk6eZNVuHs4PLWeLXoNREQtgQGIqA2qHGmVrE9Gij5F+jfNkIYSc4nV47y0XtVrc9zDEewWDHs1a3OISDkYgIhaMX2pXgo4yfnJUtC5YrxidaSVndoOoW6h1RbujNBFQOeoa/FrICJqjRiAiGysckh5cn4yUgwpUtBJ1ifjesl1q8e5O7gjUheJCF2E9G+4LhzBrsGwU/NHm4ioNvwtSdRCTBYT0gvSkZKfUq3pqrYh5f7O/ojURSLSI1IKOhG6CM6bQ0TUCAxARE2syFQk1eRIzVf6ZKQb0lEuyms8xk5lhxD3EES4R0hBJ1IXiXBdOFzsXVr8GoiI2jsGIKIGEELgesl1qQanatCpbe4cJzsnqcmqavNViBuHlBMRtSQGIKJaWIQFV41XZU1WlX119KV6q8d5ab1kQacy7Pi7+EOtUrfoNRARUXUMQEQAysxl0uSAyfpkqZ9OqiEVpebSGo9RQYUg1yB5bY5HJCLcI+Ch9WjxayAiorpjACJFKSgruFmLU6Xp6rLxMizCUuMx9mp7hLmHVeuIHOYeBic7LohLRNQWMQBRuyOEQHZRtmxIeWXYySnOsXqcm70bIjwiZH1zInWRCHYNhkbNGb+JiNoTBiBqs8ot5bhccLnakPIUfQqMJqPV4/yc/KSgUzXs+Dj5cFg5EZFCMABRq1dcXoxUfao0yqoy5FwyXILJUvNitBqVBiFuIQjXhVfriOzq4Nri10BERK0LAxC1GnklefLRVhWdka8WXrV6jFajlSYGrFqbE+oeCgeNQ4uWn4iI2g4GIGpRlauV37q2VYo+BXmleVaP83T0lAWdys7IAS4BHFZORET1xgBEzcJkNiGtIK1a0Ek1pKK4vNjqcUEuQTV2RPbUerZo+YmIqH1jAKJGMZYZZR2QK/9NL0iHWZhrPMZObYcwt7Abc+ZUCTlh7mFwtndu8WsgIiLlYQCi2xJC4FrJNdkq5ZVBJ7so2+pxLvYussU7K4NOB7cOXK2ciIhsip9CJDFbzLhivFK9I7I+BQVlBVaP83HykTVZVf7r5+zHYeVERNQqMQApUEl5CS4ZLskCTrI+GZf0l1BmKavxGLVKjQ6uHeQ1OhVNWO4O7i1+DURERI3BANSO6Uv1sgU8KwPPVeNVCIgaj3HUOCLc/cbcOREeN5uuwtzD4KhxbPFrICIiag4MQG2cEAJZRVnSCuVVg871kutWj9M56mpstgp0CeSyD0RE1O4xALURJosJ6QXp0irlVefPKSovsnpcgEuAbBbkyqDjpfVi/xwiIlIsBqBWpshUVG1YebI+GemGdJSL8hqPsVPZIdQ9VF6b4xGJCPcIDisnIiKqAQOQDQghcL3kerW5c5L1ycgszLR6nLOdc7UmqwiPCIS4hcBebd+i10BERNSWMQC1oH8l/gv/vPBPpBhSoC/VW93PW+stW/Kh8v/+zv5stiIiImoCDEAtKK8kDydyTgAAVFAh2DVYaqqqXNsqQhcBnaPO1kUlIiJq1xiAWtDgkMHwd/GXhpVr7bS2LhIREZEitYpltNesWYPw8HBotVrExsbi8OHDVvf95JNPMHDgQHh6esLT0xPx8fHV9hdCYOHChQgMDISTkxPi4+Nx8eLFFriS2kXoIjAiYgQ6e3Vm+CEiIrIhmwegzZs3Y9asWVi0aBGOHTuGnj17IiEhAdnZNa8xtXv3bkyYMAE///wzDh48iJCQEAwfPhxXrlyR9nnvvffw4YcfYt26dTh06BBcXFyQkJCAkpKSFrwyIiIiaq1UQoiapwRuIbGxsejfvz9Wr14NALBYLAgJCcFLL72EuXPn3vZ4s9kMT09PrF69GhMnToQQAkFBQXjllVfw6quvAgD0ej38/f2xYcMGjB8//rbnNBgM0Ol00Ov1cHfnMg9ERERtQX0+v21aA1RWVoajR48iPj7+ZoHUasTHx+PgwYN1OkdRURFMJhO8vLwAACkpKcjMzJSdU6fTITY21uo5S0tLYTAYZF9ERETUftk0AOXm5sJsNsPf31+23d/fH5mZ1ufDqWrOnDkICgqSAk/lcfU555IlS6DT6aSvkJCQBl4RERERtQU27wPUGEuXLsWmTZuwdetWaLUN71Q8b9486PV66Ss9Pb1Jy0lERESti02Hwfv4+ECj0SArK0u2PSsrCwEBAbUeu3z5cixduhQ//vgjevToIW2vPC4rKwuBgYGyc/bq1avGczk6OsLRkSudExERKYVNa4AcHBzQt29f7Nq1S9pmsViwa9cuxMXFWT3uvffew1tvvYUdO3agX79+suciIiIQEBAgO6fBYMChQ4dqPScREREph80nQpw1axYmTZqEfv36ISYmBqtWrUJhYSGmTJkCAJg4cSKCg4OxZMkSAMCyZcuwcOFCbNy4EeHh4VK/HldXV7i6ukKlUmHmzJl4++23ER0djYiICLz++usICgrCmDFjbHqtRERE1DrYPACNGzcOOTk5WLhwITIzM9GrVy/s2LFD6sSclpYGtfpmRdXatWtRVlaGRx99VHaeRYsW4Y033gAAzJ49G4WFhXj66aeRn5+Pe+65Bzt27GhUPyEiIiJqP2w+D1BrxHmAiIiI2p42Mw8QERERkS0wABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHi2DwArVmzBuHh4dBqtYiNjcXhw4et7nv69Gk88sgjCA8Ph0qlwqpVq6rt88Ybb0ClUsm+unTp0sxXQURERG2JTQPQ5s2bMWvWLCxatAjHjh1Dz549kZCQgOzs7Br3LyoqQmRkJJYuXYqAgACr573jjjuQkZEhff3yyy/NeBVERETU1tg0AK1YsQLTp0/HlClT0K1bN6xbtw7Ozs749NNPa9y/f//+eP/99zF+/Hg4OjpaPa+dnR0CAgKkLx8fn2a8CiIiImprbBaAysrKcPToUcTHx98sjFqN+Ph4HDx4sFHnvnjxIoKCghAZGYnHH38caWlpte5fWloKg8Eg+yIiIqL2y2YBKDc3F2azGf7+/rLt/v7+yMzMbPB5Y2NjsWHDBuzYsQNr165FSkoKBg4ciIKCAqvHLFmyBDqdTvoKCQlp8OsTERFR62fzTtBNbcSIEXjsscfQo0cPJCQkYPv27cjPz8fXX39t9Zh58+ZBr9dLX+np6S1aZiIiImpZdrZ6YR8fH2g0GmRlZcm2Z2Vl1drBub48PDzQqVMnJCYmWt3H0dGx1j5FRERE1L7UuwbIZDLBzs4Of/zxR6Ne2MHBAX379sWuXbukbRaLBbt27UJcXFyjzl2V0WhEUlISAgMDm+ycRERE1LbVuwbI3t4eoaGhMJvNjX7xWbNmYdKkSejXrx9iYmKwatUqFBYWYsqUKQCAiRMnIjg4GEuWLAEqOk6fOXNG+v+VK1dw4sQJuLq6IioqCgDw6quv4sEHH0RYWBiuXr2KRYsWQaPRYMKECY0uLxEREbUPDWoCe+211zB//nx88cUX8PLyavCLjxs3Djk5OVi4cCEyMzPRq1cv7NixQ+oYnZaWBrX6ZiXV1atX0bt3b+nx8uXLsXz5cgwePBi7d+8GAFy+fBkTJkzAtWvX4Ovri3vuuQe//vorfH19G1xOIiIial9UQghR34N69+6NxMREmEwmhIWFwcXFRfb8sWPHmrKMLc5gMECn00Gv18Pd3d3WxSEiIqI6qM/nd4NqgMaMGdPQshERERHZXINqgNo71gARERG1Pc1eA1Tp6NGjOHv2LFCx/lbV/jlERERErVWDAlB2djbGjx+P3bt3w8PDAwCQn5+Pe++9F5s2bWKHYyIiImrVGjQT9EsvvYSCggKcPn0a169fx/Xr1/HHH3/AYDDg5ZdfbvpSEhERETWhBvUB0ul0+PHHH9G/f3/Z9sOHD2P48OHIz89vyjK2OPYBIiIianvq8/ndoBogi8UCe3v7atvt7e1hsVgackoiIiKiFtOgAHTfffdhxowZuHr1qrTtypUr+N///V8MHTq0KctHRERE1OQaFIBWr14Ng8GA8PBwdOzYER07dkRERAQMBgP++te/Nn0piYiIiJpQg0aBhYSE4NixY/jxxx9x7tw5AEDXrl0RHx/f1OUjIiIianL1DkAmkwlOTk44ceIEhg0bhmHDhjVPyYiIiIiaSb2bwJpyNXgiIiIiW2hQH6DK1eCvX7/e9CUiIiIiamYN6gO0evVqJCYmIigoqF2uBk9ERETtG1eDJyIiIsWpdwAqLy+HSqXC1KlT0aFDh+YpFREREVEzqncfIDs7O7z//vsoLy9vnhIRERERNbMGzwS9Z8+epi8NERERUQtoUB+gESNGYO7cuTh16hT69u1brRP0Qw891FTlIyIiImpyDVoNXq22XnGkUqna/BxBXA2eiIio7anP53eDaoC44jsRERG1ZfXqA/TAAw9Ar9dLj5cuXYr8/Hzp8bVr19CtW7emLSERERFRE6tXAPr+++9RWloqPX733Xdls0GXl5fj/PnzTVtCIiIioiZWrwB0a3ehBnQfIiIiIrK5Bg2DJyIiImrL6hWAVCoVVCpVtW1EREREbUm9RoEJITB58mQ4OjoCAEpKSvDss89K8wBV7R9ERERE1FrVKwBNmjRJ9viJJ56ots/EiRMbXyoiIiKiZlSvALR+/frmKwkRERFRC2EnaCIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcmwegNWvWIDw8HFqtFrGxsTh8+LDVfU+fPo1HHnkE4eHhUKlUWLVqVaPPSURERMpj0wC0efNmzJo1C4sWLcKxY8fQs2dPJCQkIDs7u8b9i4qKEBkZiaVLlyIgIKBJzklERETKoxJCCFu9eGxsLPr374/Vq1cDACwWC0JCQvDSSy9h7ty5tR4bHh6OmTNnYubMmY0+Z2lpKUpLS6XHBoMBISEh0Ov1cHd3b4IrJSIiouZmMBig0+nq9PltsxqgsrIyHD16FPHx8TcLo1YjPj4eBw8ebNFzLlmyBDqdTvoKCQlp0OsTERFR22CzAJSbmwuz2Qx/f3/Zdn9/f2RmZrboOefNmwe9Xi99paenN+j1iYiIqG2ws3UBWgNHR0c4OjrauhhERETUQmxWA+Tj4wONRoOsrCzZ9qysLKsdnG1xTiIiImp/bBaAHBwc0LdvX+zatUvaZrFYsGvXLsTFxbWacxIREVH7Y9MmsFmzZmHSpEno168fYmJisGrVKhQWFmLKlCkAgIkTJyI4OBhLliwBKjo5nzlzRvr/lStXcOLECbi6uiIqKqpO5yQiIiKyaQAaN24ccnJysHDhQmRmZqJXr17YsWOH1Ik5LS0NavXNSqqrV6+id+/e0uPly5dj+fLlGDx4MHbv3l2ncxIRERHZdB6g1qo+8wgQERFR69Am5gEiIiIishUGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcLoXRgk6m52PvhRxE+Log0scV4T7OcHbgt4CIiKil8dO3Be1PysVfdl6QbQvUaRHh4yJ9dfR1RYSPCzp4OsFOwwo6IiKi5sAA1IK6BLjh0b4dkJJbiJTcQlwvLEOGvgQZ+hIcSLom29dOrUKotzMiK4JRZEUwivRxga+bI1Qqlc2ug4iIqK3jRIg1aKmJEPOLypCcW4iUnEIpFCXnFiIl14gSk8XqcS4OGkT4uiDCxxWRPi6I9L1Zg+SmtW+28hIREbVm9fn8ZgCqga1ngrZYBDINJTcDUc6NUJScW4j060Ww1PId83F1RKSvi1RzFFERkEK8nOFop2nJyyAiImpRDECNZOsAVJuycgvSrhdV1BgZkZJbiKSKGqScglKrx6lVQIiX881Q5HOjBinC1wWB7lqo1WxSIyKito0BqJFacwCqTUGJCam5RUjONSK5SrNaSm4hjKXlVo/T2qsR7l21Kc1VqkXycHZo0WsgIiJqKAagRmqrAcgaIQRyjKWyUJRc0ayWdr0IJrP1t4Cns321UBTh64Jwbxdo7dmkRkRErQcDUCO1twBUm3KzBZfzimUdsCsDUoa+pNZjgz2cZP2MbjStuSLY0wkaNqkREVELYwBqJCUFoNoUlZUjNfdmf6PkimCUnGOEocR6k5qDRo1Qb2cpGEn9jXxc4OPqwCH8RETULOrz+c15gMgqZwc7dAtyR7cg+ZtICIG8ItONUJRTtUmtECnXClFWbkFithGJ2cZq53TT2lUZoeZaMSu2C8J9XODqyLcjERG1DNYA1YA1QA1nsQhc1RfL+hpVNq1dzitGbe82f3fHm/2NKkOSrwtCvZxhz1mxiYjoNtgE1kgMQM2jxGRG+vUiadh+ZX+jlNxC5BrLrB6nUasQWmUIf+Uw/khfV/i7c1ZsIiK6gU1g1Cpp7TWI9ndDtL9btef0xaaboUiqNbrxVVRmlv5/Kyd7jVRTdOuyITonzopNREQ1Yw1QDVgD1HoIIZBlKEVyZW1Rzs0lQ9KuF8Fcy7TY3i4OVUapuUqdskO9nDmEn4ioHWITWCMxALUNJrMF6dKs2IUVo9RuBKUsg/VZsVWqm0P4O1YEo8qvIA8O4SciaqsYgBqJAajtKywtl82EXRmMknMKUVDLrNgOdmqEezsjsmKE2s1lQ1zg5cIh/ERErRn7AJHiuTja4c5gHe4M1sm2CyFwrbBMCkXJVZrVLl0rQlm5BReyjLiQVX0Iv7vWDpG+8hFqlTVHzg78USIiaktYA1QD1gApk9kicDW/GEk5xltqjwpxVV/7EP5AnVY+Ss33xnD+EE8n2HEIPxFRi2ATWCMxANGtSkxmpF4rrDZCLTnHiLwik9Xj7NQqhHo7yyZ/rJwd29eNQ/iJiJoSm8CImpjWXoMuAe7oElD9ByqvsAwp1242paXkFiIpx4jUa4UoMVkqlg+pPoTfxUFT0YzmWtEh+0ZICvdxgbuWQ/iJiJoTa4BqwBogagoWi0CmoeTmQrM5hdJw/vTrRahlBD98XB0rJnuUN6uFeDnD0Y5D+ImIasImsEZiAKLmVlZuQVrFEH5phFpF7VFOgfUh/GoV0MHTWQpGkVXWVAt010LNIfxEpGAMQI3EAES2VFBiQmpuEZKrLDZbGZQKy8xWj/NxdcDAaF8M6uSDgdG+8HF1bNFyExHZGgNQIzEAUWskhECOsfSWUHRj+ZC060UwmeU/yncEuWNQJ18MivZF3zBPONhxNBoRtW8MQI3EAERtTVm5BUcv5WHvxRzsvZCD01cNsuddHDSI6+gtBaJwHxeblZWIqLkwADUSAxC1dTkFpfglMQd7L+Ri38Uc5BrLZM+HejljUCcfDIr2xYAoH7g6ckAoEbV9DECNxABE7YnFInAmwyDVDh29lCdrLrNTq9AnzBODK2qH7ghyZ2dqImqTGIAaiQGI2jNjaTl+TbomBaLUa0Wy571dHHBP9I3aoYGdfODnprVZWYmI6oMBqJEYgEhJ0q4VYU9FGDqYdA3GWxaL7RrojkGdfDA42hd9wz05DxERtVoMQI3EAERKZTJbcEzqTJ2LU1f0suedHTS4K9Ibg6J9MKiTLyJ8XLicBxG1GgxAjcQARHTDNWMpfknMxZ4LOdh3MbfaJI0dPJ2kkWUDory5hAcR2RQDUCMxABFVJ4TA2YwCqe/Qb6l5KDNbpOc1ahX6hHpgULQvBnXyRfdgHTtTE1GLYgBqJAYgotsrKivHr8nXsPdCLvZeyEFyrnzBV09ne9wT7YtB0T4Y3MkXfu7sTE1EzYsBqJEYgIjqL/16kVQ7dCDxGgpu6UzdJcBNai7rH8HO1ETU9BiAGokBiKhxTGYLTqTnY++FG4Ho9yt6VP1No7VXV3SmvtFc1tGXnamJqPEYgBqJAYioaV0vLMMviblSIMq+pTN1sMeNztSDO/lgQJQPO1MTUYMwADUSAxBR8xFC4HxWAfacz8Heizk4klK9M3XvEI8bzWUVnak17ExNRHXAANRIDEBELaeorByHkq9jz4UbgSg5R96Z2sPZHvdE+Uj9hwJ07ExNRDVjAGokBiAi27mcVySNLNuflIuCEnln6s7+bjcWcu3ki/7hXtDaszM1Ed3AANRIDEBErUN5lc7Uey7m4vfL+dU6U8dGeEv9hzr6urIzNVErVW62IO16ERKzjUjMMSIm3Av9wr2a9DUYgBqJAYiodcqr2pn6Yg6yDPLO1EE6rdR36O6OPtA5szM1UUsrMZmRlGNEYrYRSRVhJzHbiNTcIll/v+eHdMTs+7s06WszADUSAxBR6yeEwIUsoxSGDqVcR1n5zV+uahXQq0pn6p4dPNiZmqgJ6YtMSMwpuFGjU/mVY8TlvGJYSxZaezU6+roiys8VCXcE4IHugU1aJgagRmIAImp7isvMOJRSMTP1xRwkZhtlz+ucKjtT3+g/FKhzsllZidoKIQSyC0rlIaci6Ny6NmBVOid7RPvdCDpRfq7o6OeKKF9XBHs4NesSOQxAjcQARNT2Xckvxr6K2qFfLubCcEtn6mg/V6l2KDaCnalJ2cwWgfQq/XMqg05StrHarO5VBbhrq4WcKD9X+Lg62KQ/XpsLQGvWrMH777+PzMxM9OzZE3/9618RExNjdf8tW7bg9ddfR2pqKqKjo7Fs2TI88MAD0vOTJ0/GZ599JjsmISEBO3bsqFN5GICI2pdyswUnL+ul5rKT6fmwVPnN52inRkyEFwZXBKJoP3ampvapxGRGSm6hrCYnKduI5NxCWRNyVWoVEObtIjVdSYHH1wVurWzS0jYVgDZv3oyJEydi3bp1iI2NxapVq7BlyxacP38efn5+1fY/cOAABg0ahCVLlmDUqFHYuHEjli1bhmPHjuHOO+8EKgJQVlYW1q9fLx3n6OgIT0/POpWJAYiofcsvKsP+xGtSIMrQl8ieD9RpMTD6RlPZPVE+8HB2sFlZiRrCUGK60QG5SshJzDYi7XqRLPxX5WinRmRlyKkSdsJ9nNvM2n1tKgDFxsaif//+WL16NQDAYrEgJCQEL730EubOnVtt/3HjxqGwsBDffvuttO2uu+5Cr169sG7dOqAiAOXn52Pbtm0NKhMDEJFyCCGQmG2smIgxF4eSr6H0ls7UPTp4SEPte3bwgJ1GbdMyE6HivZtjLL052qpK89WtIySrctPaIcrPVdZHJ8rXDcGeTm1+oEB9Pr/tWqxUNSgrK8PRo0cxb948aZtarUZ8fDwOHjxY4zEHDx7ErFmzZNsSEhKqhZ3du3fDz88Pnp6euO+++/D222/D29u7xnOWlpaitPTmm8VgMDTyyoiorVCpVIj2d0O0vxueGhiJEpMZh1OuS7VDF7KMOJGejxPp+fhw10W4a+1wd+XM1J18EezBztTUvCwWgct5xdVHXGUbq/Vtq8rPzVHWZFVZq+Pr5sgmXlsHoNzcXJjNZvj7+8u2+/v749y5czUek5mZWeP+mZmZ0uP7778fDz/8MCIiIpCUlIT58+djxIgROHjwIDSa6tV4S5YsweLFi5vsuoio7dLaa6RwAwAZ+mLsu5CLPRWdqfXFJnz3Rya+++PG75woP9eKVe19EBvhDSeHttFUQK1PWbkFqddu9s+5WPFvco5RVitZlVoFhHg5S+Gmo9Q/xxU6p9bVP6e1sWkAai7jx4+X/t+9e3f06NEDHTt2xO7duzF06NBq+8+bN09Wq2QwGBASEtJi5SWi1itQ54Sx/UMwtn8IzBaB3y/nS0Ptj6flSR9Wn+5PgYOdGrERXhWByBed/NmZmqozlpZXa7JKyjbi0vUimK100HHQqBHp6yIbaRXl54oIHxeOYGwgmwYgHx8faDQaZGVlybZnZWUhICCgxmMCAgLqtT8AREZGwsfHB4mJiTUGIEdHRzg6Ojb4OohIGTRqFXqHeqJ3qCdmxEdDX2zCgcQbYWjvhdwbQ+8v5mLfxVy8s/0sAtzlnak9XdiZWimEELhWWCZrrqqcHfnWTvdVuTraoWPV/jkVYSfEy7nN989pbWwagBwcHNC3b1/s2rULY8aMASo6Qe/atQsvvvhijcfExcVh165dmDlzprRt586diIuLs/o6ly9fxrVr1xAY2LQzThKRsumc7DGieyBGdA+EEAJJOYU3OlNfyMGhlGvINJRgy9HL2HL0MlQVnakHVwSiXiHsTN0eWCwCV/XFNU4UmF9ksnqcj6sjovxcqoQcN0T5ucLfnf1zWorNR4Ft3rwZkyZNwkcffYSYmBisWrUKX3/9Nc6dOwd/f39MnDgRwcHBWLJkCVAxDH7w4MFYunQpRo4ciU2bNuHdd9+VhsEbjUYsXrwYjzzyCAICApCUlITZs2ejoKAAp06dqlNND0eBEVFjlZjMOJJa0Zn6Qi7OZxXInnfT2uHujpWdqX3QwdPZZmWl2zOZLbhUpX/OzeHlhSg2mWs8RqUCOng6yZqsKkdccZ265tFmRoGhYlh7Tk4OFi5ciMzMTPTq1Qs7duyQOjqnpaVBrb75V9KAAQOwceNGLFiwAPPnz0d0dDS2bdsmzQGk0Wjw+++/47PPPkN+fj6CgoIwfPhwvPXWW2zmIqIWo7XXYGC0LwZG++K1kUCmvqSiqSwHvyTmIr/IhB2nM7Hj9I3O1JG+LhgU7YvBnXwRG+kFZweb/3pWpKKyciRlF1YbcXXpWhHKrfTPsdeoEOFzszansiNypI8rO8W3YjavAWqNWANERM3JbBE4daViZuoLOTieni/r/OqgUaN/hKfUmbpLgBubRZpYXmEZEnOMuJglnyzwSn6x1WNcHDRSJ+Sq/XRCvZzZnNlKtKmJEFsjBiAiakn6YhMOJuViz4Vc7L2QU+1D2M/NEQMrhtoPjPaFFztT14kQAhn6ElmTVeWIq2uFZVaP83ZxkGpxqjZfBeq0DKKtHANQIzEAEZGtCCGQnFso1Q79mnxd1sdEpQK6B+uk2qHeoR6wV3jtQ7nZgkuVC3lWzopcUaNTWFZz/xwACPZwkvfNqQg8HK3XdjEANRIDEBG1FqXlZvyWmoe9F3Kw50IOzmXe0pna0Q5xHb0rlurwRYhX++1MXVxmRlLOzeHklV+p1wphMtf8UWanViHcx6VaR+RIXxf2s2qHGIAaiQGIiFqrbEMJ9l7MlTpTX7+lKSfCxwWDKoba3xXpDRfHtvchry8yVV/2IceIy3nFsPaJ5WSvQUe/6kEnzNtF8TVkSsIA1EgMQETUFlgsAn9c1UtD7Y+l5clGKtlrVOgX5iUNte8W6N5q+rAIIZBlKK0IOAVS/5zE7ELkGq0v5OnpbC+Fm45Vwk6QzglqThSoeAxAjcQARERtUUGJCQeSrkkLuaZfl3em9nF1lGqHBkb7wNu1+acGMVsE0qr0z6mszUnONqKg1PpCnkE67c2OyFU6I7dEmantYgBqJAYgImrrhBBIvVYkdaY+mHwNRbd0CO4erMOgTj4YFO2LPmGejWoqKjGZkZJbeMskgUYk5xSizFzzQp4atQphXs7VOiJ39HVtk013ZHsMQI3EAERE7U1puRlHL+XdWMj1Qg7OZBhkz7tW7Uwd7YtQ75o7UxtKTPLRVhVhJ/16EazMEwitvRqRPq7Vgk64twsc7Ng/h5oOA1AjMQARUXuXXVCCXyo6U++7mFttXpxwb2cM6uSLCB8XWc1OdoH1/jk6J/tqc+dE+bki2IP9c6hlMAA1EgMQESmJxSJwJsMgLeR69FKe1WUfACDAXXuzuapK4PFxdWg1naxJmRiAGokBiIiUzFhajoMVnamzDCWI9K3aP8cFblou5EmtU5taDJWIiFoXV0c7DOvmj2Hd/G1dFKJmw95nREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOAxAREREpDgMQERERKQ4DEBERESkOK0iAK1Zswbh4eHQarWIjY3F4cOHa91/y5Yt6NKlC7RaLbp3747t27fLnhdCYOHChQgMDISTkxPi4+Nx8eLFZr4KIiIiaivsbF2AzZs3Y9asWVi3bh1iY2OxatUqJCQk4Pz58/Dz86u2/4EDBzBhwgQsWbIEo0aNwsaNGzFmzBgcO3YMd955JwDgvffew4cffojPPvsMEREReP3115GQkIAzZ85Aq9Xa4Cor/PENcOKrxp1DiCYoSCPP0RrK0O7K0YRUKgCquv2rUlfZ1phjVYAKdT8WFcfLttXn9W89tp6vX63s9Xj9Go+tzzlqKntdr6G2Y+twjlrLXss1QADCcvO9LsQt20QdtqGO+zXlscLKNksdtzX0WNTjNapug/X9Wuz+1OVY1PP7V+W+VL3uvpOAu55rgV+KNVMJYdvf3rGxsejfvz9Wr14NALBYLAgJCcFLL72EuXPnVtt/3LhxKCwsxLfffittu+uuu9CrVy+sW7cOQggEBQXhlVdewauvvgoA0Ov18Pf3x4YNGzB+/Phq5ywtLUVpaan02GAwICQkBHq9Hu7u7k13sftWALsWN935iIiI2qp7ZgHxi5r0lAaDATqdrk6f3zatASorK8PRo0cxb948aZtarUZ8fDwOHjxY4zEHDx7ErFmzZNsSEhKwbds2AEBKSgoyMzMRHx8vPa/T6RAbG4uDBw/WGICWLFmCxYtbIJh0SgDcAht/nsq/0hp3krZfhiYrRxNoFffjlr8ia/wLso5/Fd7uHFb/Qq3hL706vd5tXv+2r2ft9ev6enX4y785j5Wuz9p9q8v3rx7H1lRjUdsxNdUc1al2ykqNVb2PRT1fQ11DbVZzlq+Rx9a1tq6mWtc67WeltrdFj0X1bR4hjfyd1zg2DUC5ubkwm83w9/eXbff398e5c+dqPCYzM7PG/TMzM6XnK7dZ2+dW8+bNk4WqyhqgJud/x40vIiIisimb9wFqDRwdHeHo6GjrYhAREVELsekoMB8fH2g0GmRlZcm2Z2VlISAgoMZjAgICat2/8t/6nJOIiIiUxaYByMHBAX379sWuXbukbRaLBbt27UJcXFyNx8TFxcn2B4CdO3dK+0dERCAgIEC2j8FgwKFDh6yek4iIiJTF5k1gs2bNwqRJk9CvXz/ExMRg1apVKCwsxJQpUwAAEydORHBwMJYsWQIAmDFjBgYPHoy//OUvGDlyJDZt2oTffvsNH3/8MQBApVJh5syZePvttxEdHS0Ngw8KCsKYMWNseq1ERETUOtg8AI0bNw45OTlYuHAhMjMz0atXL+zYsUPqxJyWlga1+mZF1YABA7Bx40YsWLAA8+fPR3R0NLZt2ybNAQQAs2fPRmFhIZ5++mnk5+fjnnvuwY4dO2w7BxARERG1GjafB6g1qs88AkRERNQ61Ofzu1UshUFERETUkhiAiIiISHEYgIiIiEhxGICIiIhIcRiAiIiISHEYgIiIiEhxGICIiIhIcRiAiIiISHFsPhN0a1Q5N6TBYLB1UYiIiKiOKj+36zLHMwNQDQoKCgAAISEhti4KERER1VNBQQF0Ol2t+3ApjBpYLBZcvXoVbm5uUKlUTXpug8GAkJAQpKenc5mN2+C9qjveq7rjvao73qu6472qu+a8V0IIFBQUICgoSLaOaE1YA1QDtVqNDh06NOtruLu784ekjniv6o73qu54r+qO96rueK/qrrnu1e1qfiqxEzQREREpDgMQERERKQ4DUAtzdHTEokWL4OjoaOuitHq8V3XHe1V3vFd1x3tVd7xXddda7hU7QRMREZHisAaIiIiIFIcBiIiIiBSHAYiIiIgUhwGIiIiIFIcBqAnt3bsXDz74IIKCgqBSqbBt27bbHrN792706dMHjo6OiIqKwoYNG1qkrLZW33u1e/duqFSqal+ZmZktVmZbWbJkCfr37w83Nzf4+flhzJgxOH/+/G2P27JlC7p06QKtVovu3btj+/btLVJeW2rIvdqwYUO195VWq22xMtvK2rVr0aNHD2kyuri4OHz33Xe1HqPE9xQacK+U+p6qydKlS6FSqTBz5sxa97PFe4sBqAkVFhaiZ8+eWLNmTZ32T0lJwciRI3HvvffixIkTmDlzJp566il8//33zV5WW6vvvap0/vx5ZGRkSF9+fn7NVsbWYs+ePXjhhRfw66+/YufOnTCZTBg+fDgKCwutHnPgwAFMmDAB06ZNw/HjxzFmzBiMGTMGf/zxR4uWvaU15F6hYkbaqu+rS5cutViZbaVDhw5YunQpjh49it9++w333XcfRo8ejdOnT9e4v1LfU2jAvYJC31O3OnLkCD766CP06NGj1v1s9t4S1CwAiK1bt9a6z+zZs8Udd9wh2zZu3DiRkJDQzKVrXepyr37++WcBQOTl5bVYuVqr7OxsAUDs2bPH6j5jx44VI0eOlG2LjY0VzzzzTAuUsPWoy71av3690Ol0LVqu1srT01P87W9/q/E5vqfkartXfE8JUVBQIKKjo8XOnTvF4MGDxYwZM6zua6v3FmuAbOjgwYOIj4+XbUtISMDBgwdtVqbWrlevXggMDMSwYcOwf/9+WxfHJvR6PQDAy8vL6j58b91Ql3sFAEajEWFhYQgJCbntX/btkdlsxqZNm1BYWIi4uLga9+F76oa63CvwPYUXXngBI0eOrPaeqYmt3ltcDNWGMjMz4e/vL9vm7+8Pg8GA4uJiODk52axsrU1gYCDWrVuHfv36obS0FH/7298wZMgQHDp0CH369LF18VqMxWLBzJkzcffdd+POO++0up+195YS+kxVquu96ty5Mz799FP06NEDer0ey5cvx4ABA3D69OlmXxTZ1k6dOoW4uDiUlJTA1dUVW7duRbdu3WrcV+nvqfrcKyW/pwBg06ZNOHbsGI4cOVKn/W313mIAojahc+fO6Ny5s/R4wIABSEpKwsqVK/HFF1/YtGwt6YUXXsAff/yBX375xdZFafXqeq/i4uJkf8kPGDAAXbt2xUcffYS33nqrBUpqO507d8aJEyeg1+vxz3/+E5MmTcKePXusfrArWX3ulZLfU+np6ZgxYwZ27tzZ6jt+MwDZUEBAALKysmTbsrKy4O7uztqfOoiJiVFUEHjxxRfx7bffYu/evbf9K9LaeysgIKCZS9k61Ode3cre3h69e/dGYmJis5WvtXBwcEBUVBQAoG/fvjhy5Ag++OADfPTRR9X2Vfp7qj736lZKek8dPXoU2dnZspp5s9mMvXv3YvXq1SgtLYVGo5EdY6v3FvsA2VBcXBx27dol27Zz585a25XpphMnTiAwMNDWxWh2Qgi8+OKL2Lp1K3766SdERETc9hilvrcacq9uZTabcerUKUW8t25lsVhQWlpa43NKfU9ZU9u9upWS3lNDhw7FqVOncOLECemrX79+ePzxx3HixIlq4Qe2fG81axdrhSkoKBDHjx8Xx48fFwDEihUrxPHjx8WlS5eEEELMnTtXPPnkk9L+ycnJwtnZWfz5z38WZ8+eFWvWrBEajUbs2LHDhlfRMup7r1auXCm2bdsmLl68KE6dOiVmzJgh1Gq1+PHHH214FS3jueeeEzqdTuzevVtkZGRIX0VFRdI+Tz75pJg7d670eP/+/cLOzk4sX75cnD17VixatEjY29uLU6dO2egqWkZD7tXixYvF999/L5KSksTRo0fF+PHjhVarFadPn7bRVbSMuXPnij179oiUlBTx+++/i7lz5wqVSiV++OEHIfiekqnvvVLqe8qaW0eBtZb3FgNQE6ocqn3r16RJk4QQQkyaNEkMHjy42jG9evUSDg4OIjIyUqxfv95GpW9Z9b1Xy5YtEx07dhRarVZ4eXmJIUOGiJ9++smGV9ByarpPAGTvlcGDB0v3rtLXX38tOnXqJBwcHMQdd9wh/vvf/9qg9C2rIfdq5syZIjQ0VDg4OAh/f3/xwAMPiGPHjtnoClrO1KlTRVhYmHBwcBC+vr5i6NCh0ge64HtKpr73SqnvKWtuDUCt5b2lEjd+aRAREREpBvsAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARERGR4jAAERERkeIwABEREZHiMAARNdLu3buhUqmQn59f52PCw8OxatWqZi1Xc1OpVNi2bVuTnW/y5MkYM2ZMk52vqrKyMkRFReHAgQNNds4hQ4Zg5syZTXa+5lTf99uGDRvg4eHRrGVqLdatW4cHH3zQ1sUgG2AAonZt8uTJUKlUePbZZ6s998ILL0ClUmHy5Mk2KdvtGAwGvPbaa+jSpQu0Wi0CAgIQHx+Pb775Bu1xAvcPPvgAGzZskB43ZcBYt24dIiIiMGDAAGmbSqWSvnQ6He6++2789NNPTfJ6rc2RI0fw9NNPN+k5q94/FxcXREdHY/LkyTh69GiTvk5zmzp1Ko4dO4Z9+/bZuijUwhiAqN0LCQnBpk2bUFxcLG0rKSnBxo0bERoaatOyWZOfn48BAwbg888/x7x583Ds2DHs3bsX48aNw+zZs6HX621dxCan0+mapdZBCIHVq1dj2rRp1Z5bv349MjIysH//fvj4+GDUqFFITk5u8jLYmq+vL5ydnZv8vJX37/Tp01izZg2MRiNiY2Px+eefN/lr3cpkMjXJeRwcHPCnP/0JH374YZOcj9oOBiBq9/r06YOQkBB888030rZvvvkGoaGh6N27t2zf0tJSvPzyy/Dz84NWq8U999yDI0eOyPbZvn07OnXqBCcnJ9x7771ITU2t9pq//PILBg4cCCcnJ4SEhODll19GYWFhncs8f/58pKam4tChQ5g0aRK6deuGTp06Yfr06Thx4gRcXV0BAHl5eZg4cSI8PT3h7OyMESNG4OLFi9J5Kpsyvv32W3Tu3BnOzs549NFHUVRUhM8++wzh4eHw9PTEyy+/DLPZLB0XHh6Ot956CxMmTICLiwuCg4OxZs2aWsucnp6OsWPHwsPDA15eXhg9erR0b86dOwdnZ2ds3LhR2v/rr7+Gk5MTzpw5A9zSBDZ58mTs2bMHH3zwgVTLkJKSgqioKCxfvlz2uidOnIBKpUJiYmKN5Tp69CiSkpIwcuTIas95eHggICAAd955J9auXYvi4mLs3LkTALBnzx7ExMTA0dERgYGBmDt3LsrLy2t8jTfffBN33nlnte29evXC66+/Lru+5cuXIzAwEN7e3njhhRdkH+TN+f2s2gS2YsUKdO/eHS4uLggJCcHzzz8Po9FY47XVpvL+hYeHY/jw4fjnP/+Jxx9/HC+++CLy8vKk/W7385CRkYGRI0fCyckJERER2LhxY7Uyq1QqrF27Fg899BBcXFzwzjvvAAD+9a9/oU+fPtBqtYiMjMTixYtl36f8/Hw89dRT8PX1hbu7O+677z6cPHlSdh0PPvgg/v3vf8v+SCIFaPblVolsaNKkSWL06NFixYoVYujQodL2oUOHipUrV4rRo0fLViV++eWXRVBQkNi+fbs4ffq0mDRpkvD09BTXrl0TQgiRlpYmHB0dxaxZs8S5c+fEl19+Kfz9/QUAkZeXJ4QQIjExUbi4uIiVK1eKCxcuiP3794vevXuLyZMnS68TFhYmVq5cWWOZzWaz8PT0FE8//fRtr++hhx4SXbt2FXv37hUnTpwQCQkJIioqSpSVlQkhhFi/fr2wt7cXw4YNE8eOHRN79uwR3t7eYvjw4WLs2LHi9OnT4j//+Y9wcHAQmzZtkpXPzc1NLFmyRJw/f158+OGHQqPRyFbABiC2bt0qhBCirKxMdO3aVUydOlX8/vvv4syZM+JPf/qT6Ny5sygtLRVCCLFmzRqh0+nEpUuXRHp6uvD09BQffPBBte+VEELk5+eLuLg4MX36dJGRkSEyMjJEeXm5eOedd0S3bt1k9+Dll18WgwYNsnqPVqxYIbp06VJte9XyCyHE9evXBQDx4YcfisuXLwtnZ2fx/PPPi7Nnz4qtW7cKHx8fsWjRImn/qitcp6enC7VaLQ4fPiw9f+zYMaFSqURSUpJ0fe7u7uLZZ58VZ8+eFf/5z3+Es7Oz+Pjjj1vk+1n1/bZy5Urx008/iZSUFLFr1y7RuXNn8dxzz0nPr1+/Xuh0Oqv3tKb7V+n48eMCgNi8ebMQdfx5iI+PF7169RK//vqrOHr0qBg8eLBwcnKSlRmA8PPzE59++qlISkoSly5dEnv37hXu7u5iw4YNIikpSfzwww8iPDxcvPHGG7JzP/jgg+LIkSPiwoUL4pVXXhHe3t7Sz7QQQhQWFgq1Wi1+/vnnWq+Z2hcGIGrXKj9Us7OzhaOjo0hNTRWpqalCq9WKnJwcWQAyGo3C3t5efPXVV9LxZWVlIigoSLz33ntCCCHmzZtX7QN4zpw5sgA0bdq0auFl3759Qq1Wi+LiYiFuE4CysrIEALFixYpar+3ChQsCgNi/f7+0LTc3Vzg5OYmvv/5aiIoPMgAiMTFR2ueZZ54Rzs7OoqCgQNqWkJAgnnnmGelxWFiYuP/++2WvN27cODFixAjpcdUPwC+++EJ07txZWCwW6fnS0lLh5OQkvv/+e2nbyJEjxcCBA8XQoUPF8OHDZftXDUDiloBR6cqVK0Kj0YhDhw4JUfH98fHxERs2bLB6n2bMmCHuu+++aturlr+wsFA8//zzQqPRiJMnT4r58+dXu541a9YIV1dXYTabayzfiBEjZCHipZdeEkOGDJFdX1hYmCgvL5e2PfbYY2LcuHFCtMD309r7TQghtmzZIry9vaXHjQlAxcXFAoBYtmyZEHX4eTh79qwAII4cOSI9f/HiRQGgWgCaOXOm7DxDhw4V7777rmzbF198IQIDA6XXcXd3FyUlJbJ9OnbsKD766CPZNk9Pz1rfR9T+2Nm6BoqoJfj6+mLkyJHYsGEDhBAYOXIkfHx8ZPskJSXBZDLh7rvvlrbZ29sjJiYGZ8+eBQCcPXsWsbGxsuPi4uJkj0+ePInff/8dX331lbRNCAGLxYKUlBR07dq11rLWtYPz2bNnYWdnJyuPt7c3OnfuLJUXAJydndGxY0fpsb+/P8LDw6VmtMpt2dnZtV5XXFyc1ZFEJ0+eRGJiItzc3GTbS0pKkJSUJD3+9NNP0alTJ6jVapw+fRoqlapO11opKCgII0eOxKeffoqYmBj85z//QWlpKR577DGrxxQXF0Or1db43IQJE6DRaFBcXAxfX1/8/e9/R48ePfDGG28gLi5OVr67774bRqMRly9frrHv2PTp0zF16lSsWLECarUaGzduxMqVK2X73HHHHdBoNNLjwMBAnDp1CmiB72dVP/74I5YsWYJz587BYDCgvLwcJSUlKCoqanRfocr3b+W9u93Pw4ULF2BnZ4c+ffpIz0dFRcHT07Paufv16yd7fPLkSezfv19qDgMAs9ksXcvJkydhNBrh7e0tO664uFj2vgQAJycnFBUVNeraqW1hACLFmDp1Kl588UUAuG1/lsYwGo145pln8PLLL1d7ri6drn19feHh4YFz5841SXns7e1lj1UqVY3bLBZLg1/DaDSib9++sg+5Sr6+vtL/T548icLCQqjVamRkZCAwMLDer/XUU0/hySefxMqVK7F+/XqMGzeu1g9tHx8fKWTcauXKlYiPj4dOp5OVsyEefPBBODo6YuvWrXBwcIDJZMKjjz4q26cp7ntjv5+pqakYNWoUnnvuObzzzjvw8vLCL7/8gmnTpqGsrKzRAagyrEVERAB1+Hm4cOFCnc/t4uIie2w0GrF48WI8/PDD1fbVarUwGo0IDAzE7t27qz1/a4f769evN/o9QG0LAxApxv3334+ysjKoVCokJCRUe75jx45wcHDA/v37ERYWBlSMNDly5Ig0HLtr167497//LTvu119/lT3u06cPzpw5g6ioqAaVU61WY/z48fjiiy+waNEiBAUFyZ43Go3QarXo2rUrysvLcejQIWl497Vr13D+/Hl069atQa9d1a3X9euvv1qtverTpw82b94MPz8/uLu717jP9evXMXnyZLz22mvIyMjA448/jmPHjsHJyanG/R0cHGQdeSs98MADcHFxwdq1a7Fjxw7s3bu31uvo3bs31q5dCyFEtRqngICAGr9PXbt2xf/93//Jjtm/fz/c3NzQoUOHGl/Hzs4OkyZNwvr16+Hg4IDx48dbvbaaNPf3s9LRo0dhsVjwl7/8BWr1jXEwX3/9dZOdf9WqVXB3d0d8fDxQh5+Hzp07o7y8HMePH0ffvn0BAImJibJO1Nb06dMH58+ft3ruPn36IDMzE3Z2dggPD7d6nqSkJJSUlFQbFEHtG0eBkWJoNBqcPXsWZ86ckTVDVHJxccFzzz2HP//5z9ixYwfOnDmD6dOno6ioSBpC/eyzz+LixYv485//jPPnz2Pjxo2yuWsAYM6cOThw4ABefPFFnDhxAhcvXsS//vUvqfapLt555x2EhIRIQ4rPnDmDixcv4tNPP0Xv3r1hNBoRHR2N0aNHY/r06fjll19w8uRJPPHEEwgODsbo0aMbfb/279+P9957DxcuXMCaNWuwZcsWzJgxo8Z9H3/8cfj4+GD06NHYt28fUlJSsHv3brz88su4fPkyUHHvQkJCsGDBAqxYsQJmsxmvvvqq1dcPDw/HoUOHkJqaitzcXKlGQ6PRYPLkyZg3bx6io6OrNdXd6t5774XRaMTp06frfO3PP/880tPT8dJLL+HcuXP417/+hUWLFmHWrFlSaKjJU089hZ9++gk7duzA1KlT6/x6AJr9+1kpKioKJpMJf/3rX5GcnIwvvvgC69ata9C58vPzkZmZiUuXLmHnzp149NFHsXHjRqxdu1aqYbndz0OXLl0QHx+Pp59+GocPH8bx48fx9NNPw8nJ6bZNpAsXLsTnn3+OxYsX4/Tp0zh79iw2bdqEBQsWAADi4+MRFxeHMWPG4IcffkBqaioOHDiA1157Db/99pt0nn379iEyMlLWtEjtHwMQKYq7u7vVGgoAWLp0KR555BE8+eST6NOnDxITE/H9999L/RFCQ0Pxf//3f9i2bRt69uyJdevW4d1335Wdo0ePHtizZw8uXLiAgQMHonfv3li4cGG1mpzaeHl54ddff8UTTzyBt99+G71798bAgQPxj3/8A++//z50Oh1QMQ9L3759MWrUKMTFxUEIge3bt1drEmmIV155Bb/99ht69+6Nt99+GytWrKix5gwV/VL27t2L0NBQPPzww+jatSumTZuGkpISuLu74/PPP8f27dvxxRdfwM7ODi4uLvjyyy/xySef4LvvvqvxnK+++io0Gg26desGX19fpKWlSc9VNtdMmTLlttfh7e2N//mf/6mxec6a4OBgbN++HYcPH0bPnj3x7LPPYtq0adIHqzXR0dEYMGAAunTpUq2vWF005/ezUs+ePbFixQosW7YMd955J7766issWbKkQeeaMmUKAgMD0aVLFzz33HNwdXXF4cOH8ac//Unapy4/D59//jn8/f0xaNAg/M///A+mT58ONzc3q323KiUkJODbb7/FDz/8gP79++Ouu+7CypUrpRpclUqF7du3Y9CgQZgyZQo6deqE8ePH49KlS/D395fO849//APTp09v0D2gtksl2uOUskTUKOHh4Zg5c2arXeph3759GDp0KNLT02UfZNb8/vvvGDZsGJKSkmSdhZuaEALR0dF4/vnnMWvWrGZ7nfbu8uXLCAkJwY8//oihQ4c262udPn0a9913Hy5cuCD9YUHKwD5ARNRmlJaWIicnB2+88QYee+yxOoUfVNRCLFu2DCkpKejevXuzlC0nJwebNm1CZmZmnWqm6KaffvoJRqMR3bt3R0ZGBmbPno3w8HAMGjSo2V87IyMDn3/+OcOPAjEAEVGb8Y9//APTpk1Dr1696r3cQnOv+ebn5wcfHx98/PHHNQ7hJutMJhPmz5+P5ORkuLm5YcCAAfjqq6+atOnPmsrO2qQ8bAIjIiIixWEnaCIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSHAYgIiIiUhwGICIiIlIcBiAiIiJSnP8Psv17mYtOMDcAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "n = 100\n", - "bootstraps = 1000\n", + "n = 40\n", + "bootstraps = 100\n", "\n", "x = np.linspace(-3, 3, n)\n", - "y = np.exp(-(x**2)) + 1.5 * np.exp(-((x - 2) ** 2)) + np.random.normal(0, 0.1)\n", + "y_nf = np.exp(-(x**2)) + 1.5 * np.exp(-((x - 2) ** 2))\n", "\n", "biases = []\n", "variances = []\n", "mses = []\n", "\n", - "# for p in range(1, 5):\n", - "# predictions = ...\n", - "# targets = ...\n", - "#\n", - "# X = ...\n", - "# X_train, X_test, y_train, y_test = ...\n", - "# for b in range(bootstraps):\n", - "# X_train_re, y_train_re = ...\n", - "#\n", - "# # fit your model on the sampled data\n", - "#\n", - "# # make predictions on the test data\n", - "# predictions[b, :] =\n", - "# targets[b, :] =\n", - "#\n", - "# biases.append(...)\n", - "# variances.append(...)\n", - "# mses.append(...)" + "p_degrees = list(range(1, 5))\n", + "for p in p_degrees:\n", + " predictions = np.zeros((bootstraps, int(n*0.2)), dtype=float)\n", + " targets = np.zeros((bootstraps, int(n*0.2)), dtype=float)\n", + " targets_nf = np.zeros((bootstraps, int(n*0.2)), dtype=float)\n", + " for b in range(bootstraps):\n", + " x_sample, y_sample = resample(x, y_nf)\n", + " X = PolynomialFeatures(degree=p).fit_transform(x_sample.reshape(-1, 1))\n", + " X_train, X_test, y_train, y_test_nf = train_test_split(X, y_sample, test_size=0.2, shuffle=False)\n", + " y_train = y_train + np.random.normal(0, 0.1, size=y_train.shape)\n", + " y_test = y_test_nf + np.random.normal(0, 0.1, size=y_test_nf.shape)\n", + " model = LinearRegression().fit(X_train, y_train)\n", + "\n", + " predictions[b, :] = model.predict(X_test)\n", + " targets[b, :] = y_test\n", + " targets_nf[b, :] = y_test_nf\n", + "\n", + " mse, bias, variance = calculate_key_metrics(predictions=predictions, targets=targets, noise_free_targets=targets_nf)\n", + " mses.append(mse)\n", + " biases.append(bias)\n", + " variances.append(variance)\n", + "\n", + "plt.plot(p_degrees, mses, label=\"MSE\")\n", + "plt.plot(p_degrees, biases, label=\"Bias^2\")\n", + "plt.plot(p_degrees, variances, label=\"Variance\")\n", + "#plt.plot(range(1, 5), np.array(biases) + np.array(variances), label=\"Bias^2 + Variance\", linestyle=\"dashed\")\n", + "plt.xlabel(\"Model Complexity (Polynomial Degree)\")\n", + "plt.ylabel(\"Error\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "id": "18d104bd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "np.float64(0.5624956986315015)" + ] + }, + "execution_count": 79, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.var((predictions))" ] }, { @@ -463,7 +708,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "lecture-materials", "language": "python", "name": "python3" }, @@ -477,7 +722,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.15" + "version": "3.13.7" } }, "nbformat": 4,