{
 "metadata": {
  "name": "",
  "signature": "sha256:2a5a2d0a90558796b39210f4feb2bf1f0cb288979f67139864fdf50f7331add0"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Covariate Shift"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Most of the probabilistic prediction analysis assume that the distribution of the training set is not different than the distribution of the test. However, this may not be the case in real life. This might be due to various reasons. For example, if you collect data in an engineering lab, your samples might be biased towards male subjects because of the underrepresented ratio of females in engineering departments. When you collect your test data from outside the lab, you will probably have equal number of females and males. This phenomenon is called covariance shift in statistics. This might cause a problem when you validate your model. Shimodaira et. al, 2000 proposed to use weighted maximum likelihood to overcome this problem. Below, I will illustrate their toy example."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's generate a training data set $\\left(x^{(n)}, y^{(n)} \\right)$ of size $n=100$ where $$y = -x+x^3+\\epsilon$$ $$\\epsilon \\sim \\mathcal{N}\\left( 0, 0.3^2 \\right)$$ and density $q_0 \\left( x \\right)$ of covariate $x$: $$\\mathcal{N} \\left(0.5, 0.5^2 \\right)$$ Let's generate test data set from normal distribution $$q_1 = \\mathcal{N}\\left(0, 0.5^2\\right)$$\n",
      "\n",
      "Shimodaira showed that the optimum weights are the ratio of densities of training to test data sets is: $$w(x) = \\frac{q_1(x)}{q_0(x)}$$\n",
      "\n",
      "The first figure compares the regression analysis results of orginary least square (OLS) and weighted least square (WLS) using the weights defined as above. It is not easy to tell which model is better. We know that OLS is the optimum unbiased estimator but we don't know if tis model will decrese the error in the test set. \n",
      "\n",
      "The second figure shows OLS fit to the test data. Clearly, this regression curve is much more similar to the regression curve computed from WLS. In other words, WLS gave almost equivalent result to the future OLS by using past data. You can read the original paper for details: \"Improving predictive inference under covariate shift by weighting the log-likelihood function\" by Shimodaira, Hidetoshi."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "n = 100\n",
      "mu_x0, sigma_x0 = 0.5, 0.5\n",
      "x = np.random.normal(mu_x0, sigma_x0, n)\n",
      "mu_e, sigma_e = 0., 0.3\n",
      "e = np.random.normal(mu_e, sigma_e, n)\n",
      "y = -x + x**3 + e\n",
      "plt.figure()\n",
      "plt.plot(x,y,'*r')\n",
      "plt.ylabel('x')\n",
      "plt.ylabel('y')\n",
      "min_x, max_x = np.min(x), np.max(x)\n",
      "x_range = np.linspace(min_x, max_x, n)\n",
      "y_true = -x_range + x_range**3\n",
      "plt.plot(x_range,y_true,'b',label='true')\n",
      "\n",
      "x_add_ones =np.column_stack((x, ones(shape(x)[0])))\n",
      "x_range_add_ones =np.column_stack((x_range, ones(shape(x_range)[0])))\n",
      "\n",
      "\n",
      "# OLS for d=1\n",
      "theta_0 = np.dot(numpy.linalg.inv(np.dot(np.transpose(x_add_ones),x_add_ones)),\n",
      "                 np.dot(y,x_add_ones))\n",
      "y_0 = np.dot(x_range_add_ones,theta_0)\n",
      "plt.plot(x_range,y_0,'g',label='OLS')\n",
      "\n",
      "# compute weights from q1(x)/q0(x)\n",
      "mu_w, sigma_w = -0.28, 0.38\n",
      "w = exp(-((x-mu_w)**2)/float(2*sigma_w**2))\n",
      "x_w = x * w\n",
      "x_w_add_ones  =np.column_stack((x_w, ones(shape(x_w)[0])))\n",
      "\n",
      "# weighted log-likelihood\n",
      "theta_1 = np.dot(numpy.linalg.inv(np.dot(np.transpose(x_w_add_ones),x_w_add_ones)),\n",
      "                 np.dot(y,x_w_add_ones))\n",
      "y_1 = np.dot(x_range_add_ones,theta_1)\n",
      "plt.plot(x_range,y_1,'m',label='WLS')\n",
      "plt.legend()\n",
      "plt.title('Training data set')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "WARNING: pylab import has clobbered these variables: ['e']\n",
        "`%matplotlib` prevents importing * from pylab and numpy\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 121,
       "text": [
        "<matplotlib.text.Text at 0x10d199910>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4U2X7wPFvOiiFQumCUpqyWqaIIAIqo76iojgYivqK\nIir4vqhsZOhPwEEZMkQ2CiqIytBXVIaKlCGKIFOWLcOWUqC0UEpLZ57fHycN3Q1t0pOU+3Ndvcg4\nOblzEp77PPOAEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBClctU7ACHKYT3gAhy08bblZQKWA5cq\n4L2EEMJpXQVSzH8mIC3P/ad1jMuWTEAjK7ZrYN7Wxa7RFBYOxFbwewon4KZ3AOKm45Xn9ingReCX\nIrZzA7IrJCL9GfQOQAio+DMTIYoTDpwBXgfigY+BWsD3wAUgCfgOqJfnNZFoCQXgeWAHMN287Umg\nexm3bQhsA64APwHz0JqAijMaOGuO/4UCz/UA9gHJQAwwIc9z28z/XkarIXUAGqMlyItAArAC8C7h\nvWcB5837Pwi0ND/uAbwP/AOcAxYAVYHqwAYgyPyeV4DAEvYvhBAV4hTwL/PtcCALiADc0QovX6CX\n+bYXsAr4Js/rt3C9AH4eyEQr9A3Af4C4Mm77GzANrZZyN1ph+1kxn6E7WoHbAqgGrCR/k1FXrhfS\nrczbPma+X5/CTUaNgXvNx8Af2IpW6BflAWAPUNN8vynXC/dZwP/QkqoXsA6YnCcmaTISQjiUggkh\nA6hSwva3oZ3R5ypYyEflea4aWmFb+wa3DUFLTFXzPL+c4msIS7le0AKEUXIfwmxgpvl2A0rvQ+gJ\n7C3muXuA42g1i7z7MKD11eSN4U60mhBIH4IohjQZCUeSgHbmnqsasAg4jXaWvhWt+aS4NvdzeW6n\nmf/1KmrDErYNQks66XmeL6nwrFvg+ZgCz3dAS0YX0JqGXgb8SthfHeBLtOanZLREVNz2W4C5aE1a\n59GOVQ0gAO3Y/Yk20ukSWjORfwnvK4QkBOFQVIH7I4EmQHu0RNAVLRnYsxM2Hq2pyjPPYyGlbJ/3\n+YLbrkRruglGa75ZyPX/dwU/L2i1jRzgFrTP/Cwl/z/9EGiH1mTVBK0/IwG4Zn7Mx/xXi+tNS0W9\nrxCSEIRD80Ir2JLRCukJJW9uE/+gtctPRGvHvxN4mOIL0VVoTVDN0c7KC8bohXaGnomW2P6dZ18J\naE1GjQtsn4rW2VsPrYAvTju0Gog7Wi0nHS2ZKGAJWvNUgHnbesD95tvn0WodNREiD0kIwpEULHRn\no52pXwR2ojV7FFcwqyKeK+u2z6AlgkTgHeAr8jdl5bXRHOcvwN/A5gL7Ggy8jVbA/595X7nSgPeA\nX9GaqdoDk4C2aEnwO2BtCZ+jJrDY/NrTaMdpuvm5MUA08Lt5Xz+h1SAAjgFfoPUpJCGjjIQDcUUb\nlved3oEIUYyvqJjaiRC6coQawlDgCNKuKRxHO7RmHBfgQeBRtH4AIYQdBQM/ow2fkxqCcBQPo40W\nSkVrXumvbzhC3BxWA23QRo9IQhBCCB3p2WT0MNrY7H3IWi5CCKE7PQviyWhjrLPRZoXWRBtR8Vzu\nBo0bN1YnTpzQJzohhHBeJ4BQvYMoq+KajJQzmzBhgt4hlIszx+/MsSsl8evN2eOnjIN0HGGUUS4Z\nZSSEEDpylISwFW1onxBC3BSUUkwbOxbthN4xOEpCqJTCw8P1DqFcnDl+Z44dJH69VUT8m9auJX7+\nfH78+mu7v5e1HH10j3Kk7CmEEOW1YtEivpwzh9ZZWbwbFcWbYWEccHfnqSFD6PfyyzZ5D4PBAGUo\n3yUhCCGchq+vL5cuXdI7DIfh4+NDUlJSocclIQghKj2DweBQbe56K+54lDUhSB+CEEIIQBKCEEII\nM0kIQgghAEkIQgghzCQhCCGEDTRo0IBffvlF7zDKRRKCEELYQEkjoLKzsys4mrKRhCCEEOX07LPP\nEhMTwyOPPEKNGjWYPn06Li4uLF26lPr169OtWze2bt2K0WjM97oGDRqwefNmQFvKYsqUKYSGhuLv\n78+TTz5Z4XMuJCEIIUQ5LV++nJCQEL7//ntSUlLo27cvANu2bePYsWNs3Lix2PkC5jkDzJkzh3Xr\n1rFt2zbi4+Px8fHhlVdeqdDP4Vah7yaEEHZksNFU2/LOfcst/CdOnIinp6dVr1m0aBFz584lKCgI\ngAkTJlC/fn1WrFiBi0vFnLtLQhBCVBqONom5YBNRSU6fPk2vXr3yFf5ubm6cP3+eunXr2iO8QiQh\nCCGEDRiKqJ7kfax69eqkpaVZ7ufk5JCQkGC5HxISwrJly7jzzjvtG2gJpA9BCCFsoE6dOpR0yd8m\nTZqQnp7O+vXrycrK4t133yUjI8Py/H/+8x/Gjx9PTEwMAAkJCaxbt87uceclCUEIIWxg3LhxvPvu\nu/j6+rJ27dpCNQZvb2/mz5/PSy+9RHBwMF5eXvmalIYOHcqjjz7K/fffT82aNbnzzjv5448/KvQz\nyGqnQginIaud5iernQohhLALSQhCCCEASQhCCCHMJCEIIYQAJCEIIYQw0zMhVAV2AfuBI0CEjrEI\nIcRNT8+ZyunAPUCaOY4dQCfzv0IIISqY3k1GufO4qwCuQJKOsQghxE1N74TggtZkdB7YgtZ0JIQQ\nQgd6L25nAm4DvIFNQDgQmXeDiRMnWm6Hh4cTHh5eUbEJIcQN+eSTT5gxYwYnT56kZs2a9OrVi4iI\nCLy9vZk4cSInTpxg+fLlhV63Y8cOXn/9dY4cOYKrqyvNmzdn9uzZtGvXzqr3jYyMJDIystzxO9LS\nFf8HXAPez/OYLF0hhLBw5KUrZsyYwfTp0/nss8+49957OXPmDIMHDyYhIYFff/2V9957r8iEcOXK\nFUJCQli0aBF9+/YlIyOD7du3ExgYSKtWrUp8T1svXaFnQvAHsoHLgCdaDWESsDnPNpIQhBAWjpoQ\nrly5Qr169Vi2bBmPP/645fHU1FQaNmzI1KlTiYmJITo6ulBC2LNnD/fdd1+ZLpdZmdYyqgv8gtaH\nsAv4jvzJQAghnMLOnTtJT0+nd+/e+R6vXr06Dz30ED/99FOxr23atCmurq48//zzbNy4scKvo5yX\nnn0Ih4C2Or6/EKKSMUyyTaOHmnBjtZCLFy/i7+9f5KUu69aty59//knTpk2LfG2NGjXYsWMHU6dO\nZeDAgZw7d46HHnqIJUuWULt27TLFX1Z6dyoLIYTN3GhBbiv+/v5cvHgRk8lUKCmcPXsWf3//El/f\nrFkzli1bBsDx48fp168fw4YNY+XKlXaLuSh6DzsVQgind+edd+Lh4cHatWvzPX716lU2btxIt27d\nrN5X06ZN6d+/P3/99ZetwyyVJAQhhCgnb29vJkyYwGuvvcamTZvIysri9OnT9O3bF6PRSL9+/VBK\nYTKZyMjIID09nfT0dDIyMjh+/DgzZ84kLi4OgNjYWL744gtdrq0sCUEIIWxg9OjRTJ48mVGjRuHt\n7U3Hjh2pX78+mzdvpkqVKhgMBr744gs8PT2pVq0a1apVIywsjBo1arBr1y46dOiAl5cXd955J7fe\neiszZsyo8M/gSPMQiiLDToUQFo467FQvlWnYqRBCCAciCUEIIQQgCUEIIYSZJAQhhBCAJAQhhBBm\nkhCEEEIAkhCEEEKYSUIQQggBSEIQQghhJglBCCHKKSIigoceeijfY2FhYUU+9tVXX+Hi4sLJkycL\n7SczM5ORI0diNBqpUaMGDRs2ZPjw4XaNPS9JCEIIUU5du3Zl586dlmUk4uPjyc7OZv/+/ZhMJstj\nJ06coEuXLsXuJyIigr1797J7925SUlKIjIzk9ttvr5DPAJIQhBCi3Nq1a0dWVhb79+8HYPv27dxz\nzz00adIk32OhoaHUrVu32P3s2bOHnj17EhgYCED9+vXp16+f/T+AmSQEIYQopypVqtChQwe2bt0K\nwLZt2+jcuTOdOnVi27ZtlsdKqh0AdOzYkZkzZ7JgwQIOHTpU4Qv5yRXThBCVRqQh0ib7CVfhN/ya\nrl27sm3bNoYNG8aOHTsYPnw4QUFBLFq0iGHDhrF9+3ZGjRpV4j7GjRuHj48Pn3/+OcOHD8fPz4+I\niAiee+65Mn6SGyPLXwshnIYjL3+9ZcsWnnzySY4fP84tt9xCXFwcV65coUmTJhw9epSAgABOnDhB\n/fr1cXFxITo6mkaNGhW7v4yMDD7++GOGDBnCX3/9RbNmzQptI8tfCyGEA+rYsSPJycksWbKEu+++\nG4CaNWsSFBTE4sWLCQoKon79+lbvz8PDg8GDB+Pj48PRo0ftFXY+khCEEMIGPD09adeuHTNnzszX\nV9CpUydmzpxJ165d822f91Ka6enp5OTk8MEHH7B161auXbtGdnY2n376KVevXqVNmzYV8hkkIQgh\nhI107dqVhIQEOnXqZHmsc+fOXLx4sVCHcsuWLS2X0qxWrRqffPIJ1apVY+TIkdStW5eAgAAWLFjA\n2rVradCgQYXEr2cfghH4DKgNKGAxMKfANtKHIISwcOQ+BD3Yug9Bz4QQaP7bD3gBfwI9gbyNZZIQ\nhBAWkhDyq0ydyufQkgHAVbREEKRfOEIIcXNzlD6EBkAbYJfOcQghxE3LESameQFrgKFoNYV8Jk6c\naLkdHh5OeHh4RcUlhBBOITIyksjIyHLvR++Jae7A98AGYHYRz0sfghDCQvoQ8qtMncoG4FMgEShu\nfVdJCEIIC0kI+dk6IejZZHQ30A84COwzPzYO2KhbREIIh+bj45Nb2Am042FLjn5kpYYghLCrFYsW\n8eWcObTOyuLdqChCq23Fs+Y2xk4MoN/LL+sdXpk4Yw1BCCF098ygQfj5+rJt5Ei20pXEjAYsn5nA\nw0/11ju0Cucow06FEEIXBoMBg8HAtUuX+bfnbDq6v0OVKtyUTVOSEIQQN73YqChqvbIR/9DWDPm0\nO7FRUXqHpAtHT4HShyCEsLucHGjdGiIi4JFH9I6m/Jxx6QohhHAIy5aBnx88/LDekehLaghCiJva\n1avQpAl8+y3ccYfe0diG1BCEEKIM3n8f7rmn8iSD8pAaghDipnX2LLRqBX/+CRV0DZoK4YxLV1hD\nEoIQwm5efFHrO5g2Te9IbEsmpgkhxA3YswfWr4djx/SOxHFIH4IQ4qajFAwZAu++C97eekfjOCQh\nWClxQyJ/9fmLK39c0TsUIUQ5ff45ZGbCgAF6R+JYpA/BSjmpOcR/HE/sjFg8G3kSMjYEn/tl5UUh\nnM3Vq9CsGaxaBXfdpXc09iGdyhXElGXiwpcXiJkag4u7C8YxRgIeD8DFTSpbQjiDceMgNhZWrNA7\nEvuRhFDBlEmRuD6RmCkxZMZnYhxlJPD5QFw9XfUOTQhRjCNHoGtXOHgQ6tbVOxr7kYSgo+Rfk4mZ\nGsOVP64QPCSYoMFBuNdy1zssIUQeSmkT0Pr0gdde0zsa+5KE4ABSD6cSMy2GxO8TqftCXYKHBeNR\nz0PvsIQQwPLlMHs2/PEHuFbyirwkBAeSHpNO7MxYzn92Hv9e/oS8HkK1ptX0DkuIm9alS9CihbZe\nUfv2ekdjf7KWkQOpGlKVsNlhdIjqQNWQquzrvI+/esuQVSHsSSnFtLFji7zo/OuvQ69eN0cyKA9J\nCHbk7udOgwkN6HiqI7XuqcXhJw6z/579JG1KKvJHK4Qou01r1xI/fz4/fv11vse3bIFNm2DKFJ0C\ncyLSZFSBcoesxk6LxeBmkCGrQtjAikWL+HLOHFpnZfFuVBRvhoVxwN2dp4YMofezL3PrrVrfwc10\nrQPpQ3AiSimS1icRMzWGjLgMGbIqRDkopdi4Zg3bRo4kIjaWcUYjXWfO5IE+fRgzxkBsLHzxhd5R\nVizpQ3AiBoMBvx5+tNnWhuafNSdpQxK/N/ydf977h6xLWXqHJ4RTMRgMGAwG0i9fZkSLFly7fBmD\nwcDu3QY+/RQ++EDvCJ2H3glhKXAeOKRzHLrxvtubVutacdvm20j7O41dobs4MfoEGXEZeocmhNOI\njYqi+7JlzPjrLx5ctowTR07Rvz98+CHUrq13dM5D7yajzsBV4DOgVRHPV8omo5IUHLJqHG2kerPq\neoclhMNTSjF93DhGR0QwapSBuDj48ku9o9KHM/chNAC+QxJCPlmJWcTNiyNubhzenbwJGRNCzQ41\n9Q5LCIe1cc0aNr3wAsGj1zFjQTiHDmkXv7kZSR9CJePu506Dt8xDVsNrcbivDFkVoigrFi3i4ZYt\n2T5+PBNT4N13GlHP7VU2rFmkd2hOx+FrCBMmTLDcCQ8PJzw8vGKicjCmLBMXvrpA7FQZsipEXnlH\nGcXFvsuh6q5EfOLBA3363DTL00dGRhIZGWm5P2nSJJAmo8rPMmR1SgwZZ2XIqhCgNRd98Ox6dpne\n5Okqd/HoJ/N4oE8fvcPSjVxT+SaRO2TVr4efZZXV05NOyyqr4qa2748L/FZlPr9s8SDh1Dxio6L0\nDskp6V1D+ALoCvgBF4C3gGV5npcaghUKrbI6PBiPIFllVdwcsrKgSxdtWetRo/SOxjE48yijkkhC\nuAHpMenEzojl/PLz+Pf2J2S0rLIqKr+RI+Hvv7WVTF2kSw2QhCDyyErMIm5uHHHzzENWx4ZQs70M\nWRWVzzffwPDhsHcv+PrqHY3jkIQgCslJzSF+aTyx78fi2diTkDEh+Nzvc9OMvBCV28mT0LEjfPcd\ndOigdzSORRKCKFahVVZfNxLwhAxZFc4rNRXuugtefBGGDNE7GscjCUGUKnfI6j8R/5B5NhPjaBmy\nKpyPUvD00+DhAZ98AlLhLUxmKotS5Q5ZbbujLc2X51lldfI/ZF0uepXVkq5CJYQe3n8foqNh4UJJ\nBrYmCeEmlbvKauufW5N2PI1djXcRPSq60CqrxV2Fyl7smYAkuTm/DRtg1iytM9nTU+9oKh9JCDc5\nr1u8aP5pc9rtawc5sLvVbo69dIwvJy21rA8zMyWFbePG8XDLlqxYZN/1YeyZgCo6uQnbOngQ+veH\nNWvAaNQ7GqEHJSpW5sVMdWriKbUjYIf6uf1GFVG7q1KgxhqNasPq1cpkMtnlfZcvXKh6tGihxoeF\nKROo8WFhqkeLFmr5woUOvW9RMeLjlQoJUWrlSu2+yWRSU8eMsdvv0dkBZaoGSw3hJqNKaTZx93On\nwQRtlVXTLVncmjCYNdUWUeNiQwwY7DZk9ZlBg3hl4kRM6ekYAFN6Oq9OmsQzgwY59L6F/aWmwqOP\nwksvaZ3JILU9e7EmIQwBfOwdiKgY1v5Hcq3uyunQQ7h+4UqnhY/SMeD/YEhVzn95HlO2yeZxFXcZ\nRFskIHvuW9hXVhY88QS0bAlvvpl/qeuKbMoU170HRAOrgO5U7FBVvWtelUZ5m01MJpNK+C5B7e20\nV/3W6Dd1Zv4ZlZ2WbdMYF0+erDauWaNMJpPauGaNWhIR4RT7FjfG2uaenBylnn1WqR49lMrMvP7a\n9atWqbFGY4U0ZTorythkZC0XtGTwJVpymAw0tucbmul9XCsNW/5Hurzjsjr46EG1o84Odfrd0yoz\nKdMOEYvKasPq1WpYjRpq45o1xW5jMik1apRSHTsqlZpa9OuHt2ihhpayH0eUkpGi9sfvV2sOr7Fb\nIsPOfQgm4BxwHshBa0JaA0wvy5uKimfLZhPvu71p9W0rbtt8G2l/a0NWT4w+UWjIqjNQDjQU1ZFi\nsYcbae55+23YtAm+/x6qFVifMTYqiu7LljHjr794cNkyh1zqOiUjhf3n9rP68Goitkfwwrcv0GVZ\nF+rOqEvA9ACe+foZVhxaQWpWqt6h5mNNaTAUeA5IBD4CvgGy0JJJFPatKajK+p9DD0siIghp0oT7\ne/fmx6+/JjYqipfGji33ftNj0omdGcv5z87j38sf42gj1ZtVt0HE9pd7Hd7uy5bpfkEVR4rFHlSe\nK5tFxMYyzmik68yZha5sNmUKfPopREZCnTr6xVualIwUopOiiUqKyv9vYhRXMq7Q2LcxYb5hhPmG\nWW6H+oZSr2Y9XAz2Hc9jz6UrJgFLgX+KeK4FcORG3/QGSEIog/R0OHNG+zt/HhITtb+rVyEtDa5d\nA5NJ+zMYtCUAPDzAywu8vbW/gACoXRsCA6FePesmAWUlZhE3L464ueZVVseEULODY66yumLRIr6c\nM4fWWVm8GxXFm2FhHHB356khQ+j38stFvkYpxfRx4xgdEWHTDumyxOKscpOewWjEFBvLgwWS36xZ\nMG8ebN2q/e4Kstd3UJwrGVeISsxf4Of+pWSmEOobSqhvqKWwb+zTmDC/MIJqBNm90C+JPa+YNqGE\n5+yZDEQpMjLg0CFt6d8DB+DYMTh+HBIStP9MRqN2huXvry0N7OsLwcFQtSq4umrJwGSCzEwtiaSm\naq+NioKLF7Vkcu4cxMVpSaJ+fQgNhbAwaNpUG/nRrJm2PzAPWX2rAcaRRuI/judw38N4NvLEOMaI\n7wO+Nv0PXN6C4ZlBg/Dz9WXbyJHXh6JOnlzimbllhNYdd9j0DL4ssTir3OaevLXUXBERsHQp/PJL\n0ckA7PMdJKcn5yvs857pp2Wl5Tu77xTSiQG3DSDUN5SgGkGVbqSao38ah6khVPSZSVFSU7Uzp61b\nYft2LQmEhkLbttC6NTRvrhXURqNW4NuKyaQlh9OntTVkoqO15HP4MJw4oSWKNm3gttvgjjugXTuo\nWdO8yupXF4idGguuEDImxGarrNqieSXv2WpOTAwZ3bqx4OuvC32/FXEGX9qZc2WmFEycCKtWwebN\nEBRUeJvyfgeXrl0qdJafezs9O/362b1vGGF+YZaz/kCvQKcs9GW1UzvTq303Lg6+/lpb8/2337TC\n/1//gs6dtTXgq+vcVJ+VpSWHffu0v927Yf9+CAmBu+/O/VP4HE8idmoMGXEZGEcZCRxQtlVWbVk4\n5+1TmTJ6NGfmzuXRzz8v9P1a2/ZdHvbq33F0JpN2xbPNm+Hnn7VmyqJY8x0kXUuynNlHJ0UTfen6\n7YycDEshn3u2n1vw16lexykL/ZJIQrATPdp3k5Lgiy9g5Uo4ehQeeQQeewy6ddPOvK2lV60mK0tr\nytq5E3bs0GozBgOEh0P34MtU+fIX6qXXIfi1YIIGB+Hu4271vm1dOFv7/d7MZ/DFKe/vKzMTnn8e\nYmNh3TrwKWX664bVq/n2lRdIburPpcxz+Pd/BAKrWM70s3Ky8p3d523fr129dqUr9Etizz6Em1pF\nte8qBdu2wYIFsHEjdO8O48fDffdBlSpl26e92rxL4+6u1WTatoVXX9U+24kTsGULfPzpFf6I6UQL\nXwODl1+k/tRd1H0+kEZjjHjU8yh13wWHz5piY8s169ja77ektu+bVXl+X1euQJ8+2kCGH3+8PmhB\nKUXitcTrzTqJUZYC//DZQ/CaC83q+FIjLYjLSRd54q7+DL5jMI19Gt90hb49SEIoha0LoIIyMmDF\nCpgzRztjeuUVLSmUdrZUkrxnvTNTUnhz3Dg+fOst3UatGAzw++ZFfDtnDndlZfGziuY/1Xsw7VIX\n/Gs/TruF2dy/YDepbf1pPtFIywdLbgezZeFs7fc7cNw4y21bJVdH6Jcqi/L+vk6cUDz0xEVu6RzN\nI89HEfF7/o5cIN+Z/oOhD1qaePw8/ZzqWAnbssssvhtlj2UPUlOVmj1bqXr1lOreXakff9RmZ9qC\nI07vLymmlBSl1i3PVB/cfkr9z2WHmlH9kHrn38lq61alsrLsH9uNfr+2WmmzpBm7tl7N05b7s+b3\nZTKZ1Pmr59WvMb+qT/Z9ot7Y/Ia69c1mqumM25VhnLeqNslH3bH4DvXvtf9WE7ZMUJ/t/0z9Fvub\nSkhNkGUobAA7L12hF72Pq81lZCg1d65SgYFK9e6t1J499nkfR5zev37VKnWfu7sa1rx5sTFlpWSr\n30bHqg3eO9Vir33q3pqJasDzJvXtt0qlpWnb6L30sTVLL5TEmnWlyvseBY9RefdX0IbVq9XQGl5q\nUJsm6vFmnmrEklfU+J/Hq76r+6q2i9qqmhE1le9UX9VhSQf1zNpn1DPz+qrmbR9U1ZtGqlXfXbRJ\nDKJ4OGlC6A4cQ5vxPKaI5/U+rjZjMin15ZdKNWqk1AMPKLV3r33fzxEXcxvx9NPqVU9PtWH16lJj\nysnMUfGfxasdTf9Q6+vtVq82P6dq1chRTzyh1Nhhv6nBXnUrPMnZ6roKJZ1h2+o9chPA2EGDyr2o\nYXxKvNp2eptaunepGvfzOPX4qseVcVJd5fl2VeU31U81mxqmOr59m5oUOUmtOLBC7TqzSyWmJVqO\n2QPNblctanyrbmW/ern+vXItigqAEyYEV7SF8hoA7sB+oHmBbfQ+rjZx4IBSXbsq1bq1Ups36x1N\nxStPIWcymdTF7y+qvZ33qk2+P6r/1lyomnluUTVIVk2q/6RaB72uFs5aWiE1Bls2xRVXgyvvexQ8\n1uNCQ1XnevXUaz4+JTbvxF2JU1tPb1Uf/fmRGvvTWNXnqz6q9YLWqvp71ZX/NH9150d3qn5f91OT\nIieplQdXqj/O/KGS0pJKjWf/fpMKCb6smnt8oVLxVEP9/HRvvrwZUMaEoGencnu0hHDafP9L4DHg\nqF4B2VpqKrz1Fixfrk28efll204Ysxdl487O8ozUMhgM+PXww6+HH5d/vYzrUBOP7btMfUbwhrsv\ncXWHMXx8AIGZP5OQuofxb7djyVT7dNSW1AF9o8esqI7x3H20uv32Mg9iKHisVUYGDz75JOeWLGbQ\n7WEkpZ3h2qXNRG7eY5mcdSLpBNWrVM83XPPxFo9b1uCpVbXWDR8rk0kbKPHW/6VT130CXTLm8AaK\n9KR0IoYP52JiYqVblqMy0DMh1ANi89w/A3TQKRab++UXGDgQOnbUZvQGBOgdkfVsPVzVViO1at1d\ni6yxKezsP4Hz7s/y8uWW7D4+m24Bf9Asph7vLBvArHkpBBu6kJC6hzfeuYNaN16Wlai4EU43esyK\nGrW0cc0a4ufP58Ijj5RpFJVJmTibcpZDaUfYY7zAHXf7ccklnmSPpaSMysK7WjL+qj5nz/1Bj6a9\neaLFE5YeWHgyAAAeHUlEQVTROzU9ip/gcqPJ7vRp7befkgJ791Vl4VsXOP+jLwsSE3nO15f2XbrI\n1eoclJ7jt/qg9SEMNN/vh5YQXsuzjZow4fpSSuHh4YSHh1dUfGVy7Rq8/jp8+602fLRHD70jsp49\nJ+HZaiZu3v1sWvQNZ2fHYYxqSKZpG/Nc11Df341O59vxdvUBnLzWkVtaJjFsZH0ee4x8ycEWtSCl\nFM89+CCXYmPLdcyKO+4+RiOfbdiQLz6TMhF3Ja7IJRhOJJ3Au6o3XqkeNPJuRJc295IWnYTruSxG\nj3yPGh41yvQ5rZ2lbzJpv/mJE2H0aBgxAtzcZFJfRYiMjCQyMtJyf9KkSeD4E4/z6QhszHN/HIU7\nlvVshrth+/Yp1by5Uk89pVRS6c2rDscRh6uWZv2qVepRN18133+k+oFv1OdVpqhkmquxRqNa8+k3\nasUKk+rZU6maNZV68EGlli5VKjHRNqNuNqxerYZ6eanJI0aU65jl5OSoQY89psaY9zHGGKxeH/mS\n6tWiqnpt8SA1atMo9dgXj6mW81oqz3c9VeD7garz0s7qhf+9oCZvm6xWH16t9sXvUykZKWX+LEW5\nkb6f3buV6tBBqbvuUuro0fzPOeIAh8oOJ+xD2AOEoXUqnwWeBJ7WMZ4yUwoWLtT6C2bNgmee0SZj\nORt7T8Kzh5+/+YYQ92s0XNCRA9t/o8r8U+xwf4e2Z+Oo8bcLvd/Rvo+UFPjhB5gx7SSDXgrAWLUO\nY9KeYsPo6Tc8aS/vGf2sq1fpv3Ilnhcu8EJgIDWtvPBQjimH2CuxRCdF8+3Pq4jM+p6sLopFfm6k\n1DhDlcyl3NHexMltX5KRVpX7u/Rk0pOTaOzbGK8qXrY4dKUqqu/nlffe4+Cff6KUwmAwcOGC9rv/\n3/+01Ur79weXAmsX2mNSn7APPRNCNvAqsAltxNHHOGGH8tWrMGiQ1k+wc6e2NLQzc6QlGlQJzTp5\nC+X3r13jzfHjibx8mRotg3ht12NseX0Dpo9c2fP9HkLGhhDweABPPeXCk0825JsV37B8+E9sSfsX\n/zv9PmHNMrh4zZ+YGG1RvtIULChTr12j2bBhuLq50fqOOyzHLLfQ//vi3yxYNoMGnVtw4tIJopKi\nOH35NJ6mqhjOZRJyxY2WyTkkXapFvRM+BNWoT/Yf+1ibnMw4o7fVazWVdLyKer60+0WdIOzftYsL\nixezrlVHDp7syQcfaAn36NHyza4Xwhp61rpKdeKEUi1aKDVgwPVJU0Jji8ljpc3kLdi8NXnECDXU\nyyvfEM6L319UezvtVb81+k2dmX9GZadl5xvy+YpXgJo49lc1YIBSfn5KtWmj1KRJWvNfSaHn7mNo\ny+bq+eDqasDrfdQ9nTxUrw97qB6f91BNP2yqPN7xUMaZRnXrtJbqlt7u6sUFz6q1R9aqg+cOqtTM\n1GKb6NavWlWmSYWlNYMVfL60+0pdb+75bOFC1TU4WA30C1HTGaG8XM+poJo/qPffWZnvO5k6ZozK\nycnRdeKgcM55CNbQ+7gWa8sWperUUWrePL0jcUzlaaO3tu069z0GBAaqgS4u6tnAwGK3v7zjsjr4\nyEG1o84Otfq+T9XGT9YWatPOylIqMlKpYcO0CYRGo1KDByv17XeZavjIwWrD3xvU3F1z1dANQ1Wr\nSU1VvclByn2im/IY5arqveiuBj6M+tej/ur2bsHq4Z53qY8WzCnxc5hMJvVyz55qqJdXvsL/Rtvc\nSzteBZ9/tXZt1apKFdUnIECZQPUJCFCtqlRRr9auXezxS0gwqf5PHVTVXC6ox1mlBtS5v1BfSe73\nMXnkSJvOihY3DkkIFeejj5SqXVupn37SOxLHY4uZttZ2bucWnDk5OWryiBGqt7d3qR27V/+6qo48\nd0Rt992uokZGqfQz6UoppTKzM1VUYpRa//d6Nef3Oeq19UNUl4UPKb9JYcrwfx6KYSGq5iud1N0R\n/1Fjvn1ffXP0f+rwhcMqNSO12JpK7tl+cZ9j/apV6mVXVzV55MhydbiWdrwKPj/GaFTvjRihxgQH\na/eDg7X7Rbz+wAGlBg1SqlYtpe6756R6ttodhWouud95n4AA1QPUIHd3S+KRWcn6wAk7lZ2OUvD2\n2/DZZ9oa/02a6B2R47HFcuFlWYG0zZ13cmHJkhK3z8rJ4mzgWU69eYp/ev/D0U+OEtw0mD9b/ckn\n7T8hp1GOZXJWY5/GuKVnwYbjdLwUzOvHE+lWvQn/uN3HER5glW917r1Xu0aFulqV9MuXeSEwELe4\nOOJWruTTq1d5c/x4tl+7xq0JCfni+nzxYr6cM4fQS5dYkJPD8JUreWTDBp4aMqRMQ3GtmTBXcLKb\nt8FARnJy/vvm59NirrL+x0b831QD587BSy9pF0Fat/RLQl4ZU6h/Kfc73zpiBF2A97OzAaju4cGr\nkyZJR7ITkYRgpexs+O9/tauC7dypXatYFGarkUo32rmdu334Yw+zYtVifjj2Hcd+P5NvnH7slVjq\n1ahnuXBKbIvjDBozjO7/60740nBqdapFyJgQanbQJmmpOxUbg7SL8ezmMh3TlhEywocx03px5Iji\njWHf8fnnj/DzT/cR4H+GB3vU4FLM51zZMRU4hyk9nQ5dunBvr1480KeP5XN41KrFleRkMhISMAAZ\nCQlccXVFWXExKFVEx7FSis8+/JDnli61vE/M338zbexYWt1+O/Hz55NQYLLb2iVL6JPn/op5X1Gn\n/2aOnmxH5Olsmv8ezbvTtOtx5M6uL260UO73m5GczCxfX5okJdHL15cGVo64EsJa+ta78vjhB21R\nuhTbDvWulOw97jwjO0MdTTiqvjv+nZq5c6Ya/P1gdf/y+1WjDxopj3c8VKMPGqkHlj+gXvnhFTXr\nt1nq++Pfq2MJx1RGdoZlHwX7OLKvZqvYObFqZ8hOtS98n0rcmKhMJpMaM2iQut3FRb1UpYoygRpa\nt67q0aKFGjNokOX12dnaqrUzZih1V/szqprhvPJyi1ehbt+oF545qDZvzj8vxWQyqfdGjFD/MRjU\ncFAvGwyWZqO8iuqYL6pvpqjHcuMeWrdukU13164ptXWr1oHesaM2T6NXL6VWrlTqypWSYyjKoMce\nU12Dg9VgPz9lAjXYz091DQ5Wg3r2tP6LFTZDGZuMHD11mz+bYzCZCo+xFvaRkZ3BqcunLNfEtVxA\nJSmK+JR4jN7GIi+KXr9Wfaq4Xr/EnCpwRl3SbOxnBg3i/dfH82yrocROj8XgZiCqyV5OB/1N8tpV\nTImNZVDVqhzz96dT1aq8Fx1daGbykogIjGFNaNy6Nx/P282fewxc4w4OHtSGZbZuDS1bwoGtC8j6\n83NCG2STFHOAkD59mL5iRb5jkHeGcMLFi4Xi3pySgrvBQBcvL8tj265eJUsp/uXlRafoaL5zc2N2\ntgtPu7fkcqt+tO4ynD17DOzbB82bQ9eu2tX5OncGD4/Cx2zT2rVseuEFHli6lIN79pQ4pNXe150W\n1pNLaFYASQa2lZ6dzslLJ/NdLjH3wujnrp7D6G20XBS9uX9zHm36KI19GtOgVgPcXa27DnPBNYZK\n6uPYtHYtZxfN48Cydtx/sDdJ65NwneJKncjGHEqOZVQzXzzjTvJg375cWb26yD6SvM0qU2a3t9w2\nmeDkSThwAI4cgYupXUkKeprfz3pjcMvEZ8slonqC0QjnYv7g+O4faEIC/025nQkD5uIa4Enz9o+S\nvO1nzlOH5FRvBr81k4xMF3ZGfMgO6nAsqSlt+/yHS8kebPj+ID/iR3R2IxYTQrWs09Ted5iGzf9i\nwoRWtG9f/PW5N61dS8zs2YR//jmdPD2ZmZJC/9dew/PCBd64dInJixYVeo0zTmoUhTn6t+VQNQRx\n465lXePkpZNFrr1z/up56teqX+RF0UO8Q6wu9PPKPbut26ABX334Yb4z6v1ubviFhPD0iy/y44sv\nWtbWqfb00xzcsaPYWsPINs/Sw/1ZqsZ6ce1fyWzOWk36prXlXptHKcW0sePoPzyCmBgDZ85oF5w/\nd06x949TnNwZTbV0d87ijfIx4upWjdSLV1Au7hhMWXj6eOPubiIr4RTuHtlUzYwnrGtzatXMIGHD\nIs4ZYkhNP4qfyym6mNKZBTxfpQon/fx4ecKEQjOzVyxaxBcffADnzvHdpUuMCAxk0/nz+CtFczc3\nFmVn80ZYGAeLWavJVutVifKTGoLQTVpWmuVMv2ATz4XUC9SvVd9ypn9L7Vvo1awXob6h1K9VHzcX\n2/4Ec2sEty5dyisTJ+arCXR64gnOL1nCL35++TpYY6Ki6NytW7G1BnVyHdnLetG66d1sfPpbup34\nN6pbX+6eF87wIf8mpoyzuTetXcu5BfM50F6rvbS3VCgMjH15KomRH9GmihuLMzMZXrUuu3Ny6DKw\nJ5MXLsy3ZPb1QvgksVFfaY/168R9vXoxZdQo9i5ciOc1rXTwqVGD8R98QPfHHy8Uz78HDmTNxx8T\neOkSPwFkZdGqShXw98fl7FltOe0SRo3JEhXOT2oIwippWWmcSDpRqD0/KjGKi2kXaejTsMgzfaO3\n0eaFflGK6hv49do1aiQkYPDx4eyFC3SqU4dZ8fFFrkhacEVOz6ef5lCBWkPuENJHZn9Bo2O3cOaj\nWGIyN1P3/XbcP6RXqTFaai8NG/JVKavKLp48mcTERJJXr2ZKbCyveXvT46OPbrhNfuOaNSx79ll8\nMzNJVgrfqlV5ZPnyQgX2ikWLWPzuu7SMi+MxpVhbpQp7MzNpdt999Bs4kPX9+3O1Rg28r12T1Uqd\ngNQQRLmlZqZy4tIJy5l+3maexGuJNKzVUOvA9QmlTWAb+rbsS6hvKMaaRlxd9L3yT1F9A7d37sx9\nvXtzf+/eTB09mrMff1zs3IiCw1xj/v6bLuZaw+fAj6dOcXdAALMyMhgy4T9svXyZW2vU453UzkSN\n6MDiN2dzvMUB3v9tabEFdkm1l4LxDBo/XktS5rkVruUYvhvSuzf39tIS1uZvvik0hDc3mbbNzmaW\nUrwJ7DeZaNatG13vuYfP586l32ef5Rs6KyonSQg3mauZV/N34iZFWzpyL6VfopFPI+0M3yeUdkHt\nePqWpwn1DSW4ZrDuhX5JiurUfLBPH0shW9rEteIuWJN++TJ/Nm9Ow1OnIDMTA9qEq6defZUrq1bR\nkM/4OPBXGjcZyr2R3dnS6mdav98W3wd8LfvPW3uZmZJS7IS1goW9LRYazPu5gCKbinKT6ab//hcD\nkOjmRiM3N54111YC8ix2JzWDyk0SQiWUkpFiKfTzNvFEJ0VzOf0yjXwaWc7029drzzO3PmMp9F0M\nzjuUqqQCtCyFa97XTBk9mjNz5xaa6ftoYCBn40/RyRTBDHWB+RefIrHnBbxr+9F6ahsCngi4PpPX\nXCPIKWbCWkF5C/P7e/dm+rhxlmWn8yo4tPZG5Rb2WVevMig4GM/kZDw6dODdwYPpWquWlsTGjbvh\nZcKF85E+BCd1JeNKkYV+VGIUKZkpNPZpTGPfxpbO3NxLJQbVCHLqQl8vBUfQrF2yhD4DB3Jfr16W\n5qgPk5O18fczZtKuajix02LJiMvAOMrIIe/f+GHgcySmpxfbjl9SwV7SVcs2rlnDxgEDSO/WjQVf\nf12mpFDw88VERRHcuLHMK3BSZe1DcPRv9qZOCMnpyUUX+klRXM28qk3MMp/p5xb4Yb5h1K1RVwr9\nClTSJSKTdyaz87/byT6Uzd/VN/Hq1RUMDPTiZE4O/3nnnXxn23kL/dwaQVHDZ3M7oAFLU1SnqCjW\nGQwcqVePgW++aZOzeLn0pfOSTmUndTn9cpHt+dFJ0aRlpeUbsdM5pDMDbhtAmF8Ydb3qypmagyiq\nOSrv2X73/T3YNOsbAv8viN2s4J7UX/n39I48MEjr6C3UxzBuHO8NG8atCQncumJFiR3Qh/78kx+X\nLQNgvlIMz8nhyzlzAMqdFBzpYkmiYjh6iVIpagiXrl0q8iw/Oima9Ox0y5m+pWnH/G+gV6AU+g6q\ntHb7gk08uferBd5Kw9N30NjjEer2DcI42ki1ptUsyz60jI3lAzc3bUSTeYhs7vDZsMaNC52pb1i9\nmmXPPYd3RgZLlGKYnx/dFy6Upp2bnNQQdJZ0LanIiVlRSVFk5mTmK+zvaXAPA9sOJMwvjDrV68h/\nXCdUcEmMXEWd7X/41lvUCwujd56z7VP799LQvRH7u+zHu5M3Lne4FzmiKe/w2aI6oM9ER1OrWzdc\nv/+e5/z88EpJkSUjRJk5+q/GoWoIiWmJRZ7lRyVGkW3KzneWn1v4h/mFEVAtQP6DVhIlLY7X7+WX\nb3iRt5zUHOI/jufohMO41nOh9fS2LP45gjPz5uJRRI2gqFiqnj/PS4mJbAsN5df0dFq2a8f8b74p\n9bOUd3SScFxSQ7Cz1YdXM/C7gZaVNUN9Qnmg8QMMbjdYCv2bSGkXALrRRd5cq7sSPCSYoP8GceGr\nC5x8/SRtE++h3X+7cu+0Hvy0rvBEsqJi6Z6YyNaMDMbNmlVix2/BlUyLquWIm5ckBCv1adGHx1s8\nLoX+Tc6aAr8snbEu7i4E9gukzjN1SNqQRMyUGP5o+ge3jOpIt6E9yxxLQUWtZCpzDEQuSQhWkmGc\nIldpBX55FnkzGAz4PeSH30N+JO9MJmZqDKffPk3wa8EEDQ7C3Sf/CrClxVLU+kkfZmQwPDubH0+d\nojllu8ypqJwc/XTXofoQhCivsrTbpx5JJWZaDInrEgl8IRDjcCMe9Tys2l/u6KYHli5FKWXp2xjm\n50d8aipBDRuizpyROQaVTFn7EPQ67X0COAzkAG11ikGICmdpt//6a6tfU71FdZp/0px2+9tBDuxu\ntZtjLx4j9VhqsftbsWgRD7dsyfbx45mZksL28eOJGDGClAsXGNGiBZkpKbR59VVmHj7Mg8uWyRwD\nAehXQ2gGmIBFwEhgbzHbSQ1BVAqljU66EVmJWawf8C1u66uQWvUED6QuYVpYdqmjnbI7dy60fpJc\nwKZycrZRRsd0el8hdFHa6KQb4e7nzqPf9mHj52tIfO0QR3iL22OS6DqyHt0GPgIUs/qreS4DyAVs\nRNGkp1SIClCwgL52+XK5JpAZDAYMVQ38nbOGNc3f46RhI24rqvNn2z85/8V5TNkmS4fzjL/+kmYh\nYRV71hB+AgKLeHw88J21O5k4caLldnh4OOHh4eWNSwhd2HptoIL7O/n3HnrdOpCYqTGceuMUD416\njsCHAuU6BjeByMhIIiMjy70fvUcZbUH6EISwudwhq1d2XSl2yKqovJxtlFFeeiclISod77u8afVt\nK2775TbSotLY1XgX0SOjST+TrndowoHpVRj3AuYA/kAysA94sIjtpIYghA2kx6ZzZtYZzn1yDv+e\n/hhfN1K9WXW9wxJ2IhfIEUKUKisxi7h5ccTNjcP7bm9CxoZQs0NNvcMSNubMTUZC3JSUUkwbO5aK\nPOlx93OnwVsN6HiqI7X+VYvDTx5m/z37SdyYWKFxCMckNQQhdFLSdZIriinLRMKqBGKmxIArhIwJ\nIeCJAFzc5FzRmUmTkRBOwpazlm1FKUXSem2V1Yy4DIyjjAQOCMTV01WXeET5SEIQwknc6EV0KpoM\nWXV+0ocghJOw9axlWytyyOqoaDLiMvQOTdiZJAQhdOAMy0qUtMqqqJwc45SkeNJkJISDsAxZnReH\n910yZNWRSR+CEKJC5KTlEP9xPLEzYqnaoCohY0PwfcDXYZq8hCQEIUQFswxZnRoDLjJk1ZFIQhBC\n6EIpRdKGJGKmxpARK0NWHYEkBCGE7mTIqmOQhCCEcBiph1OJmR5D4rpEAgcEEjw8mKrBVfUO66Yh\nCUEI4XDSY8yrrH56Dv9e/hhHyyqrFUESghDCYeVbZbWTNyFjZMiqPUlCEEI4vJzUHOKXypBVe5OE\nIIRwGqYsExe+ukDs1FhZZdUOJCEIIZyOZZXVqTFknMnAONI8ZLWaDFktD0kIQginZhmy+vsVgofI\nkNXykIQghKgUUo+kEjPNPGT1hUCMw4141PPQOyynIglBCFGppMemc2amechqT3+Mr8uQVWtJQhBC\nVEpZSVnEzZVVVm+EJAQhRKUmQ1atJwlBCHFTsAxZnRYrq6wWw9kSwnTgYSATOAEMAJKL2E4SghCi\nSJZVVqeYh6zKKqsWzpYQ7gM2AyZgivmxsUVsJwlBCFGqvKus1nu1HvVeqXdTD1l1toSQVy+gD9Cv\niOckIQghrCZDVjVlTQiO0Oj2ArBe7yCEEM6veovqNP+kOe0OtIMc2N1qN8deOEbqsVS9Q3MKbnbc\n909AYBGPjwe+M99+A60fYWVxO5k4caLldnh4OOHh4TYLUAhROVU1ViV0Vij136xP3Lw49nfdj/dd\n3hjHGPHu6K13eDYXGRlJZGRkufejZ5PR88BA4F4gvZhtpMlICFFuOWk5xH988wxZdbY+hO7ADKAr\ncLGE7SQhCCFsptAqq6+HENC38g1ZdbaEEAVUAZLM938DBhexnSQEIYTN5RuyGptnyGolWWXV2RKC\ntSQhCCHsKvm3ZGKmVK4hq5IQhBCiHCrTkFVJCEIIYQOVYZVVSQhCCGFDWUlZxM2LI25unNMNWZWE\nIIQQdpCTZl5l9X3nGbIqCUEIIezIlGUiYVUCMVNjHH6VVUkIQghRAZRSJK1PImaq466yKglBCCEq\nWN5VVoNfCyZocJBDDFmVhCCEEDrJN2R1QCDBw4OpGlxVt3gkIQghhM7yDVnt5Y9xtD5DViUhCCGE\ng8g3ZPVub0LGhFCzQ80Ke39JCEII4WD0GrIqCUEIIRyUZcjqlBhtlVU7D1mVhCCEEA4u3yqrZzK4\nfe/tuNey/agkSQhCCOFErv51Fa9bvOyyb0kIQgghgLInBMebcy2EEEIXkhCEEEIAkhCEEEKYSUIQ\nQggBSEIQQghhJglBCCEEIAlBCCGEmV4J4R3gALAf2AwYdYpDCCGEmV4JYRrQGrgN+B8wQac47Coy\nMlLvEMrFmeN35thB4tebs8dfVnolhJQ8t72AizrFYVfO/qNy5vidOXaQ+PXm7PGXlZuO7/0e8CyQ\nBnTUMQ4hhBDYt4bwE3CoiL9HzM+/AYQAnwCz7BiHEEIIKzjC4nYhwHrgliKeiwYaV2w4Qgjh9E4A\noTf6Ir2ajMKAKPPtx4B9xWx3wx9ICCGEc1mD1ny0H1gL1NY3HCGEEEIIIYRDegI4DOQAbUvY7jRw\nEK2p6Q/7h2U1a+PvDhxDazYbUwFxWcMXbSDA38CPQK1itjuNYx17a47lHPPzB4A2FRSXtUqLPxxI\nRjve+4A3Kyyy0i0FzqPV9ovjyMe+tPjDcdxjD9qE3i1oZc5fwJBitnPk76BEzYAmaB+ypAL1FFoB\n5misid8VrbO8AeCO1mzWvCKCK8U04HXz7THAlGK2c6Rjb82xfAht0AJAB+D3igrOCtbEHw6sq9Co\nrNcZrYAprkB15GMPpccfjuMee4BAtMm9oM3nOk45f/+OtpbRMbQzVGs4wgipgqyJvz1aIXAayAK+\nROtY19ujwKfm258CPUvY1lGOvTXHMu/n2oVW86lTQfGVxtrfgqMc74K2A5dKeN6Rjz2UHj847rEH\nOId2EgFwFTgKBBXY5oa+A0dLCNZSwM/AHmCgzrHcqHpAbJ77Z8yP6a0OWvUZ87/F/Wgc6dhbcyyL\n2ibYznFZy5r4FXAXWnV/PdCiYkKzCUc+9tZwpmPfAK22s6vA4zf0Hegx7PQntKpOQeOB76zcx91A\nPBBg3t8xtGxfEcobv7JtODekuNjfKHBfUXyceh77gqw9lgXP8vT8DvKyJo69aG3FacCDaGt/NbFn\nUDbmqMfeGs5y7L3QRm4ORaspFGT1d6BHQrjPBvuIN/+bAHyDVvWuqEKpvPHHkX91VyNa1q4IJcV+\nHi1ZnAPqAheK2U7PY1+QNcey4DbB5sccgTXx5133awMwH60PJ8m+odmEIx97azjDsXdHG7q/Ai1h\nFXRD34EjNxkV13ZXDahhvl0duJ+SRznopbj496BNzGsAVAGexDE6rtYB/c23+1P0j8vRjr01x3Id\n8Jz5dkfgMtebxvRmTfx1uP5bam++7UgFUkkc+dhbw9GPvQH4GDgCzC5mG6f+DnqhtXddQztT3WB+\nPAj4wXy7EVpHyn60oVbjKjjGklgTP2jVz+NoHYqOEr8vWt9AwWGnjn7sizqWL5v/cs01P3+Akkev\n6aG0+F9BO9b7gZ041kKQXwBngUy03/0LONexLy1+Rz72AJ0AE1p8uUNjH8S5vgMhhBBCCCGEEEII\nIYQQQgghhBBCCCGEEEIIIYQQQgghhBDO7P8BrIJ2zVwdbfYAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10cf76210>"
       ]
      }
     ],
     "prompt_number": 121
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# test set\n",
      "mu_x0, sigma_x0 = 0., 0.3\n",
      "x1 = np.random.normal(mu_x0, sigma_x0, n)\n",
      "y1 = -x1 + x1**3 + e\n",
      "\n",
      "plt.figure()\n",
      "plt.plot(x1,y1,'*r')\n",
      "plt.ylabel('x')\n",
      "plt.ylabel('y')\n",
      "min_x, max_x = np.min(x1), np.max(x1)\n",
      "x_range = np.linspace(min_x, max_x, n)\n",
      "y_true = -x_range + x_range**3\n",
      "plt.plot(x_range,y_true,'b',label='true')\n",
      "\n",
      "x1_add_ones =np.column_stack((x1, ones(shape(x1)[0])))\n",
      "x1_range_add_ones =np.column_stack((x_range, ones(shape(x1)[0])))\n",
      "\n",
      "\n",
      "# OLS for d=1\n",
      "theta_0 = np.dot(numpy.linalg.inv(np.dot(np.transpose(x1_add_ones),x1_add_ones)),\n",
      "                 np.dot(y1,x1_add_ones))\n",
      "y_0 = np.dot(x1_range_add_ones,theta_0)\n",
      "plt.plot(x_range,y_0,'g',label='OLS')\n",
      "plt.legend()\n",
      "plt.title('Test data set')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 122,
       "text": [
        "<matplotlib.text.Text at 0x10d5497d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX6wPHPsLgjgooLoKagCW6QS2reqCzNrFxKy8xK\nQS1NU+sq1u9qXXNpMTUzFXLPTNHSXHDJSLulpqi55K6J+4riAgLz/f1xBkQcYAZm53m/XvNyZjhz\nzjPjzHnOdwchhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQwq5OAE/YOwghbMHN3gEIYQE3gBTD\nTQ/cyvH45ULsLwHoY+ZrlOFmCj1Q28z9F1Utw3HlNy+KzMPeAQhhAeVy3D+OdtLfWIT9mZoAikJn\ng2M40nGFEMJhHQceN9x3A0YAR4BLwPeAj+FvpYAFhuevAtsAP+BjIAO4jVZimZLHcV4F/jG8fmSu\n4zYH/jDs9wzwJeBp+NsmtCv/rFLSi0AFYCVwAbgC/AT45/MehwOngOvAgRzH1eXzfk8ajptVEmuR\nz/6FEKJYyXkCHwz8DlRHO3FPBxYa/tYPWIGWQHRAGOBl+NsvQO98jhGCdvJ9BCgBfA6k5zhuOFry\ncANqAvsNsWTJXVXlC3Q2xFIOWAz8kMex66ElgaqGxzVy7Cu/91sTqaoSQgijciaO/TnuA1QD7gDu\nwBvA/4CGRvbxC/m3cfyHuydkgDJAWq5j5fQOsCzH44LaOJqglTyMCQLOozXEe+b6W17v1w1p4xAW\nJG0cwpXVQrty1+d4LgOtSmo+EAgsQqsqWgC8b/g75N/OUQ2tqijLLeByjsd1gYnAQ2hJxQPYns/+\nygBfAO24W7VUDq0klDuOI2iJaDQQCqwFhgJnyfv9Vsnn2EKYTa4+hCs7CbRHOxln3cqgnWQzgI/Q\nTr6tgI5AL8PrCmocP4uWdLKUASrmePw12tV/EOCNlpDy+60NQ0s2zQ3bP4qWNPJqyP4OaINW/aSA\nCYbn83u/tmjwF8WEJA7hyqYDY9HaAQAqA88Z7kegVVO5o7VXpAOZhr+dB+rks984tETTGq2N4yPu\n/S2VM+zzFvAg8Gau1+fefzm0xvhraO0do/I5dl206qiSaNVjqTnizu/9XkQrieT3voQQoljK2cah\nA4ag9Ty6jlbNM8bwt5cMz98AzgGTuHvyfxg4iNbOMCmP4/Ti3l5Vx3Ictw3wN1ry2AR8aPg3Sz+0\n3lZXgRfQqr5+MWx/AOiLlgyMXdg1BLYa3s9ltAb+rIby/N4vhjguGI7bPI/3JYTDm4V29bUnn22m\nAIeB3Wg9X4QQQhRjbdCSQV6JowOw2nC/BbDFFkEJIYRwbLXIO3FMB7rneHwA6SEihBB25eiN4/5A\nUo7Hp4AAO8UihBACx08ccH+XROlWKIQQduToAwBPc29/+QDDc/eoU6eOOnr0qM2CEkIIF3EUbbyR\nWRy9xLGCu4OyHgaS0Xph3ePo0aMopRzqNmrUKLvHIDG5VlwSk8Rk6RuFHNdj7xLHd2ijZCuhtWWM\n4u78OzPQelR1QOuPfhNtfiEhhBB2ZO/EYcoiOwOtHoUQQgiTOXpVldOKiIiwdwj3kZhM54hxSUym\nkZisz1VWA1OG+jphoJTi0+ho3hs3Dp3OVf6b7UM+S+GqDN9ns7/UUuJwUWuXLuXstGmsW7as4I1F\nvuSzdE6+vr7odDq56XT4+vpa9LN1lcsnKXEYLJgxg0VTptA4PZ0xhw/zQXAwuz09eWnQIHr262fv\n8JyKfJbOTafTIecFTV6fRWFLHPZuHBcW9krfvlT09WXTsGHoAH1qKgPHjqVd1672Ds3pyGcphHFS\nVeVisoqmqcnJDA0J4XZycvZzwjzyWQphnJQ4XFDS4cO0nz2bp7p0Yd2yZSQdPmzvkJyWfJZC3M9V\nLp2kjUMIcQ9p47jL0m0cUlUlhBA2VqtWLTZu3GjvMApNEocQQthYfqWhjIwMG0djPkkcQghhQ6++\n+ionT57k2WefxcvLi08//RQ3NzdmzZpFzZo1adu2Lb/++iuBgYH3vK5WrVr8/PPPgDYodfz48QQF\nBVGpUiW6d+/O1atXbfYeJHEIq1FK8cmIEVLPLEQO8+fPp0aNGqxcuZKUlBS6desGwKZNmzhw4ADx\n8fF5tkdk9eibMmUKK1asYNOmTZw9exYfHx8GDBhgs/cgiUNYTXEbcS2J0rnodJa5FVXW92X06NGU\nLl2aUqVKFfiaGTNmMGbMGKpXr46npyejRo0iLi4OvV5f9IBMIIlDWNyCGTPoGBrK5pEjmZiSwqbo\naDqGhrJgxgx7h2ZVxS1ROjulLHOzlNxVU/k5ceIEnTt3xsfHBx8fH0JCQvDw8OD8+fuWK7IKSRzC\n4l7p25cBo0ejT029O+L6ww95pW9fm8Vgy6v/4pooReEZG0Sa87myZcty69at7MeZmZlcvHgx+3GN\nGjWIj4/n6tWr2bdbt25RrVo16wZuIIlDWJwjjLi25dW/IyRK4VyqVKlCfstd161bl9TUVFavXk16\nejpjxowhLS0t++/9+/dn5MiRnDx5EoCLFy+yYsUKq8edRRKHsIqsEdef793L07Nn22zEtT2u/h0h\nUQrnEh0dzZgxY/D19WXp0qX3fVe8vb2ZNm0akZGRBAQEUK5cuXuqsgYPHsxzzz3HU089Rfny5WnZ\nsiXbtm2zWfyu8s2WkeMC0Kqo4uPi2DRsGOOSkogODOTRiRNp17WrWSdyc9fgiBk3jhp1694zNUnk\niBFFeSuiiGTk+F0yO64Q+ch99a9PSirU1X92VVezZibNhhsVHZ19X2bPFa6u2FVVbfpnExduXrB3\nGMKKilJNJg3dQhSs2FVVDV8/nJmJM3nigSeICo+ibe22uLu5Wzk84SwsVdUl7E+qqu6SSQ6LaMKT\nEzgx+ARta7fl/Y3vU2dKHT769SOSriXZOzThAKShW4iCFbvEAeBdypv+Tfuzve92fuj+A+dunKPJ\njCY8s/AZfvj7B9Iz0+0dorAje/UIE8JZuMplVJF7Vd1Kv0Xc/jhiEmM4cuUIrzV+jT5hfQiuGGyh\nEIUQtiRVVXdZuqpKEocRBy4dIDYxlnm75xHqF0pkWCRd6nehtGdpix1DCGFdkjjukjYOG3iw0oN8\n9tRnnBp6igHNBjDvr3kEfhHIoDWD+Ov8X/YOT1hJQdOUyCSGQmgkceSjhHsJXgh5gbU917K973Yq\nlKpAh2870CK2BbGJsaSkpdg7RGFBBU1TIpMYCuFalK2kZ6arlQdXqk6LOqkK4yuoyOWRakvSFqXX\n620WQ1Hp9Xo1Yfhwp4rZmuZPn66eCQlRI4ODlR7UyOBg9UxIiJo/fbpJfxeOyZbnhcKYPXu2atCg\ngSpTpoyqWrWqevPNN1VycrJSSqlRo0apnj17Gn3d5s2bVcuWLZW3t7fy9fVVrVu3Vn/++We+x8rr\nswCKdfG5aP+DhXTm+hk1bvM4VWdyHdVwWkM1ectkdfnWZbvEYo41S5aod7y8VHxcnL1DcQh6vV6t\nXrxYjQgMVArUiMBAtWbJkuzEWtDfhWOy13nBFJ999pmqUqWKWrt2rcrIyFAnTpxQHTp0UM2aNVN3\n7tzJM3Fcu3ZNeXt7q0WLFim9Xq9u376t1q1bp/766698j5fXZ0EhE4dUVRVBNa9qjHhkBIfePsTk\n9pPZcmoLtSfX5pVlr/DL8V8cri5cRkUbV9DYDRnbISzp+vXrjB49mqlTp/LUU0/h7u5OzZo1Wbx4\nMSdOnGDBggV5frcOHTqETqeje/fu6HQ6SpUqxZNPPknDhg1t+h4kcViAm86Nxx54jIVdF3Js8DGa\nV2/OoPhB1J1al/G/jedsyll7hwjI9N/5KWjshoztEJby+++/k5qaSpcuXe55vmzZsnTo0IH169fn\n+dp69erh7u7O66+/nr0ehz3IJIcW5lval8EPD2ZQi0FsO72NmMQYQqaFEFErgsiwSNoHtbfbFCeW\nmgDQFRU0SaFMYuh6dB9a5nuvRplXs3Dp0iUqVaqEm9v91+3VqlVjx44d1KtXz+hrvby8+O2335gw\nYQJRUVGcO3eODh06EBMTg5+fX6HiLwxXOWMoR6sWyiklLYVFexcRuzOWMyln6N2kN73DelOzQk2b\nx2Kv6b+VmdOUuzL5LGzDUcdxxMfH8+yzz5KWlnZf8njttddIT0+nXr16HDlyhPnz5+e7r4MHD9Kz\nZ0+Cg4NZuHBhnttZehyHqyhaS5UN7T63W729+m1VcUJF9dT8p9SSfUtUWkaavcOyOmmQv0s+C9tw\n1PNCcnKyKlu2rFq8ePE9z6ekpCg/Pz/1zTffqNGjR+fZqyq3L7/8UjVs2DDfbfL6LJDGcefQqEoj\npjw9haQhSfRq1Iuv/vyKwC8CeW/dexy4dMDe4VmcNMjfJZ+FAG11v1GjRvH222+zdu1a0tPTOXHi\nBN26dSMwMJCePXuilEKv15OWlkZqaiqpqamkpaVx8OBBJk6cyOnTpwFISkriu+++o2XLlnZ+V87J\npMzsqA5dOqSGrx+uqnxaRT0y6xE1d9dcdfPOTXuHZRHSlfUu+Sxsy9HPC998841q0KCBKl26tKpS\npYrq379/9jiO0aNHK51Od88tMDBQnT59WnXr1k35+/ursmXLKn9/f9W/f3+VkpKS77Hy+iwoZInD\n3nVb7YFJgDsQC0zI9fcIYDlwzPB4KTDGyH4Mn4FzS89MZ+WhlcTujOWPpD94qcFLRIVHEVYtzN6h\nFUl8XBzxvXtzpUQJfO7cocPs2cW2gTk+Lo61vXujCwxEn5TE08X4s7A2R23jsAdXmqvKHZiKljxC\ngJeB+ka2+xUIM9yMJQ2X4enuSef6nVnVYxW7+++mStkqdPq+Ew/NfIjp26dzPe26vUMslKTDh6nS\nty++aWlU69evWHdllW69whXYs8TREhiFljgAsrr2jM+xTQQwDHi2gH25RInDmEx9JhuObSAmMYYN\nxzbQuX5nosKjaBnQ0il65CyYMYNFU6bQOD2dMYcP80FwMLs9PXlp0CB69utn7/CEC5MSx12WLnHY\ncxyHP5Bz2b1TQItc2yigFbAbOA28C+y3SXQOwt3NnXZB7WgX1I4LNy8wd9dc3lj+Bu46dyLDI+nV\nuBeVylSyd5h5eqVvXyr6+rJp2LC7gw7HjpXqGSGcmD0ThymXAolAIHALeBr4EahrbMPRo0dn34+I\niCAiIqLIAToav7J+vNf6Pd5t9S6bT24mJjGGj379iHZB7YgKj+LxBx7HTedYHeVk0KEQjiMhIYGE\nhIQi78eev96HgdHcraqKBvTc30Ce03HgIeBKruddtqqqIFdvX2XhnoXEJMZwPe06fcL68HqT1/Ev\n72/v0LLZa9ChKN6kquouV1oB0AM4CDwBnAG2oTWQ/51jmyrABbTSSXNgMVDLyL6KbeLIopRix9kd\nxCbGsnjfYlrXaE1UeBQdgjvg4SYzy4jiRxLHXa6UOECrfsrqjvsNMA7IajGdAQwA3gQy0KqrhgJb\njOyn2CeOnG7cucGSfUuISYzhn2v/8Hrj1+kT3ofaPrXtHZoQNuPr62u3SQAdjY+PD1eu5K6ocd7E\nYSmSOPKw78I+YhNjWbBnAY2rNCYyPJLOD3ampEdJe4cmhLAzZxzHIWwg1C+UL9p/wakhp4gMjyQ2\nMZaALwIYEj+E/ReLTwc15ULrhbvSexHOSRJHMVHSoyQvNXiJDb02sKXPFsp4lqHtvLa0ntWa2Ttn\nc/POTXuHaFWWXi/c2MnbVid0WftcCMsozFQxxV56ZrpafmC56riwo/IZ76P6/dRP/Xn6T5eaO8la\n64Ubm+HW2rPeytrnwtKQNcdFUZy6dkqN+XWMqjWplmoyvYmaunWqunr7qr3DKjJLTyxo7OTdompV\n9Ui1alY/ocskicLSkGnVRVH4l/fn/X+9z9FBR/mk7SdsPrmZWpNq0euHXmz6Z5PT1qdbar1wZaiG\n6hEVdd/yu6OmTCF60iSzluRVhajWkrXPhaOQDv5OTim4eBGOHYOkJDh9Gs6cgQsX4PJluHIFrl2D\nmze1W2oqZGRoN70e3N3Bw0O7lSkDZcq4Ua7ck1So8CQtK13iwNl5dNnVHzePTNr6RtK93ms0DvLD\n3x88Pa31niy7Ql7WxII5ByCaK6tdYX2zZveNhHdzc0MpZdbo+Ox2imbNzJp+xRLvRYiicpVLFeWs\nV8SmUgpOnoSdO2HvXti3D/7+G44ehRIloHZtqFED/P2henXw84OKFbVb+fJQrhyULQulSt1NFG5u\nkJmp3dLT4datuwkmORmuXtUSz/nzil2X/+DPzBiSyv1AiVNtubMliuqpbakb5E5QENStCyEhEBoK\nAQFQlPN91tTj7R1gynFjkzQmXLvGvzp1Yuz06dknb6WUSaPjZdJH4UhkHIeLJY6bN2HLFvjtN+3f\n7du1k314ODRooJ2gQ0IgKAgqVLBdXNdSr7Fo7yJm7ojhXMpFHvfpTd0bvblwJJD9+2H/fi32xo2h\nSRMIC4OmTbV43d3z37cjnlSVUsTHxbFp2DDGJSURHRjIoxMn0q5r10KVhiy9PyGKwhlnxxU5pKfD\n1q2wfr12271bO/G2aQP9+0OzZlpJwt68S3nTr2k/+jXtx86zO4lNjGXSviY83PJhBr4VSce6Hbl2\n1ZO//oJdu2DjRvjkE60KLTwcHn4YWreGVq2gUq5JfR1xJl1LT9Iokz4KVyCJw46uXYM1a2D5coiP\nhwcegCefhP/+Vzuxli5t7wjzF1YtjK+e+YpPn/qUuP1xTNwykbdWv8VrjV+jT1gfHn88OHvb5GT4\n80/4/XeYOhVefVWr0oqIuHurXNkxT6qWbleQdgrh7FzlMsdpqqpu3ICffoLvvoOEBK1E8fzz8Oyz\nUK2avaMrugOXDhCbGMv8v+YTUjmEyLBIuoZ0pZRHqXu2y8zUSiQJCdpt82atncbP6w8iItJ5J7oN\nv66SmXSFsCZp43DgxKHXw6ZNMGsWrFihVdW89JKWMMqXt3d01nEn8w4rDq4gJjGGHWd20KNhD6LC\no2hYpaHR7dPTYds2rZouPl5r+H/8cejQAZ55xjGq6YRwNZI4HDBxXLoE33wDM2dqPZp694ZXXoHK\nle0dmW2dSD7B7J2zmbVrFtW9qhMZFslLDV7Cq6RXnq+5dAnWroVVq7REEhSkJdrOnbVOAaJwLN3V\nWTg3SRwOlDh27oRJk7TSRefO8OabWs+i4v47zdRnEn8kntidsSScSKBr/a5EhkfSwr9Fviex9HSt\nKmv5cli2TOta3KULvPii1nuruH+u5nCkrs7C/iRx2DlxKAXr1sGnn8KBAzBoEPTpo42jEPc7m3KW\nubvnEpsYS2nP0kSFR9GzUU98S/vm+zqltEb2pUth8WJtDEu3blrVX2iojYJ3ArlLFo7Y1VnYX2ET\nh6sweW4WvV6vJgwfbrH5ffR6pX78UamwMKVCQ5WaO1eptDSL7LpYyNRnql+O/6J6LO2hvMd5qx5L\ne6iNxzaqTH1mga/V65Xatk2pYcOUCghQqlEjpSZMUOrkydzbWfb/3BnknnBR5rkSxiCTHJrGUjOY\n6vVKLV+uJYwmTbTkIb/Borl085KavGWyajCtgaozuY4at3mcOnP9zH3bGUsEmZlKJSQo1bevUr6+\nSj3+uJbEU1KsP2utI8lvBt2sz2FISIgaXEw+D5E/JHHkz5JTUv/2m1KtWyvVoIEkDGvQ6/VqS9IW\nFbk8UlUYX0F1WtRJrTy4UmVkZiilCk4Et28rFRenVFij48rDLVk18opTCbRR0UGuPw15fiWLmWPH\nqvi4OKXX61V8XJyKGTfO3uEKO0MSR/4sUVQ/fFip559XqkYN7Wo2I6Mo/2XCFNdTr6uZ22eq5jHN\nle9HPiroxUrqrSa1TEr+er1efTtzuXrM+yMVwl5VweO4er3HX+r0adfO9FKyEKZCplXPX1GmpE5J\ngeHDtekyHn4YDh6EXr0KnntJFJ1XSS+iHopia+RWNkb9Qkibpsxuf5L2PWF/1Sv0G/VBntOX63Q6\nfH3u0Fj/KU/V70a7Er05d74coaE6OnfWuvnq9TZ+QzaQNTL98717eXr2bBmZLizOVVrTDckzfzHj\nxpk0g+ndnWo9d4YOhbZtYdw4GYhmb/Fxcazq+wZHWnizt/o5bjxQjj7NI4kMj+TBSg/et72x//Pu\nA0bw3XcwY4Y2A3DfvtoYGz8/O7whIexIuuNauDvu8ePw1ltw6pR2gmnVyqK7F4WUOxFsP/gH15u5\nM3f3XIIrBhMVHsULIS9QxrOMSfvbvh2mT9e693boAAMGQMuWMjYkNyUDB12SJA4LJY7MTJg8GcaO\nhffe00ob1lqwSFhOemY6Kw+tJHZnLFtObeGl0JeIDI8krFqYSa+/ehXmzIFp07QBhgMHQo8ejj/R\npK3IwEHXJInDAonj0CF44w1t3YtZs6BOHQtEJmwu6VoSs3bOYtauWVQuU5nI8Eh6NOxB+ZIFTwym\n12sDOb/8Ups7q08frRQSGGiDwB2QDBx0bTIAsAj0eqUmT1aqYkWlpkzRxgQI55eRmaHWHF6jun7f\nVXmP81av//i6+t/J/5nck+7wYaUGD1bKx0epF1/UumEXt67XMnDQtSG9qgrn3Dmtbvvbb7WV9t5+\nW1tSVTg/dzd32ge1J65bHIfePkRIpRDeWP4GodNCmfjHRC7dupTv64OCtDnHTpyARx6B117TetUt\nWqTNn1UcFKU3onBdxfoUuXq1trRps2baEq1BQfaOSFiLX1k/3mv9HgcGHGB6x+nsPLeToClBdI/r\nzoZjG9CrvPvlli+vzT128CCMHAlff61VY372mbYYlytRSvHJiBGoHFW/0r1X5OYqlw1KmdHGkZEB\nH3wACxdqJY02bawYmXBYV29fZeGehcQkxnA97Tp9wvrwepPX8S/vX+Brd+yAzz/Xpn5//XUYPBhq\n1LB+zNYmjeDFizSOm5g4Tp+Gl1+GMmVg/vzitzaGuJ9Sih1ndxCbGMv3+77nkRqPEBUeRYfgDni4\n5b+68smTWi+8OXPg6ae1nniNG9smbkuSRvDiSRrHTdShg1JjxkgDuCuyxCy4KWkpalbiLNUytqWq\n/nl1NXLDSHX0ytECX3f1qjYzb/XqSj35pFIbNjhXQ7o0ghdPSOO4aX78Ed5/XxrAXdHapUs5O20a\n65YtK/Q+ypUoxxthb/B7n99Z13Mdt9Jv0SK2BW3ntWXR3kWkZaQZfV2FCvDvf8OxY1qJduBAre1s\n8WJtbJCjk0ZwYQ5X+VYYkqcojqxdzZKWkcaPB34kJjGG3ed307NhT6IeiiKkct5r2Or18NNPMGEC\nXLyoVWH16gWlShU5HKsxd0oe4fykqkoUW7asZjl65agauWGkqvZZNdXqm1ZqVuIsdSPtRj6xKbVp\nk1ZFWq2aVp117Vrhj18cF6US1oNUVYniypbVLLV9avPxEx9zcshJhrcezg8HfiDwi0D6r+zP9jPb\n7+nGqsWm9dpbtQrWrNHWo69dW+vVd/Gi8e6v+clZHWfua4WwFEkcwiXYeqyBh5sHz9V7jhUvr2DP\nm3sILB9ItyXdCJ8ZzlfbviI5Nfm+1zRuDN99pw00vXgR6tWDzh2PcGjq8gLbZRbMmEHH0FA2jxzJ\nxJQUNkVH82iNGpycNKlIbTpCOKP2wAHgMDA8j22mGP6+G8hrxjq7FveEUEpbP33dkXWq25Juynuc\nt+r1Qy+16cQmo9VK86dPV08EP6paVIhVPlxWjcovVo0rNFPz8lmUKqs6bj6oph4eanC1akVezVIU\nbzhhVZU7MBUteYQALwP1c23TAQgCgoG+wNe2DFC4PmXB6h43nRtP1nmS71/4nsNvH6Zxlcb0W9mP\nB796kE//9ynnb5zP3vaVvn0Z9vEAHvP6kMMEUzr9OAeSVxE7rwO7d9+/75zVcTvq1+cBDw+4cwcd\noE9NZeCHH+a5oJU13qso3uyZOJoDR4ATQDqwCHg+1zbPAXMN97cCFYAqNopPmMBWJyNrHccSXXiN\nqVy2MkNbDmXfW/uY9dws9l/aT72p9Xhh8QvEH4lHr/TodDpSLlygX4kbNEkdwUJq43FoIS2aXiCs\n0T9s2XLvPrOq4ybu20fYgAGk37hhVpuOtd6rELb0AhCT43FP4Mtc2/wE5FxCaQPwkJF92a+sV8xl\nrW9t7XWtLX2c+dOnq2dCQtTI4GCbVfck305WX//5tXpoxkOqxhc11LP/fULNW/i1WvX996pnxYoq\nxtAjbPm3S9XUqXpVs6ZSjz1mfDDhzLFjVXxcnNLr9So+Lk7FjBtn1nvtUL++erVdO+mdVcxRyKoq\ne+qKaYmjdY7HG4BwI/uy9+df7NjqxGvOcbK6qmZmZhbYZdXULrx5dX8tarfYxDOJ6q2Vbymf8T6q\n6WdhqmNYKTUotL4a7OWl1ixZoiYMH67S0vRqzhyl6tVTqkULpVasKNxodGPvdezQoWqwDRK+cGwU\nMnG4F+ZFFuKNVjW1wPD4GeAa8L8c2zwKXAX2Gh6PRGssv5lrX6MBEhISSEhIAKBWrVqWj1hka/jQ\nQ/j6+XFi9WraXr/OhrJl6TZ+PJ1ffdWi3WDNOc7apUvZ9cEH7LpyhfPTp6OrV4+gkJDsZU9bPfFE\n9mt0Oh1H//6bg99/z7qgINIvXCC4UyeCQkKM7jNrXwU9b6pqXtV4pu4zvN3ibXYk/Mq+JnfYHp5M\nxfb/4kjCTtLnxeFRvx5dXw7hzTfBxwc+/FBb5tbXF+rXN332g5zvdZqvL7vPnaPi8eN8cfkys3fu\n5Ivp09F7eNCoadP7XmvssxPOKyEhgTlz5mSfK3/99VeAD+0dlzk8gKNALaAEsAvjjeOrDfcfBnLV\n+mazd+IulrKqj4aEhFj16rWg42SVSgb6+Sk9qL6enuoZUF0rV1bPhISo4X37Gq3myq+6J3dJJzo4\nWLWoWFG91aOH1UpaE6aOVg+87KvKDHdXEa+hnnuyimrf8MHsfev1Sq1cqVTLlkoFByv1zTdKpaWZ\ntu+s95qZmanGDh2qBnp751nSylmaslVVpLAPnLCqCuBp4CBaI3m04bl+hluWqYa/78Z4NRVI4rAL\nc+rZrXn5ihdvAAAdUklEQVScrKqY4YaqmP7u7moNqIEVKqg2AQEqOijI7JN87uqdwRUrqm6lSqnV\nixdbbZR61jHfqxmgloSgHogspcr/10sNXD1Q/XXurxzbKbVxo1Jt2yoVGKitWnnzpunHKSgRr1my\nRHUsWVI9GhCgRgYHq0xQHXx8VIf69aXLrwU5wiwAOGnisBS7ffDCMWSdDKP8/VU/UN0DAtQgLy81\nduhQoyd5U360a5YsUQNKllRdS5RQ/XQ6tcaQfB4NCFAdS5a0Skkr90l9zsJp6v82/p/y/9xfNY9p\nrmZun6mup17P3n7rVqWef16pKlWUGjtWqeTkgo+RVyLOWcrKBNWzalXVx8NDrQHVX6dTY4cNk8Z0\nC3KE0hyFTByuUmlp+AxEcZU1Qd8/Bw9y5epVKlasSI3gYDYsW0bGTz+hCwxEn5TE04YFikxZsChm\n3DgCg4PR6/V899ZbPHr5MkcDA7nzyCM82aUL7bp2tfhkgDknGly7dCkLpk5l/i+/oFd64o/EE7sz\nloQTCXSt35XI8Eha+LdAp9Oxdy+MHw/x8dCvn7awlJ+fecdWShEfF8emYcMYl5TEY2XKcP32bcI9\nPZl55w5DqlXjiI+PrNFRRI609olMciiEEbmvrt/s1MnsNgpbteXkddzcxztz/Ywau2msqjO5jmow\nrYGavGWyunzrslJKqSNHlOrXTykfH6Xeflupf/4p3DGHhISojp6eKurZZ7OrAAd6e8saHRbgSGuf\n4IQjx4WwGJXHAMGo6Gjade2KTqejXdeufLVsGQNGj0afmmryqGtbz4NlbF6qjqGhLJgxA9B6ZEW3\niebQ24eY0n4KW05tofbk2ryy7BVOuv3CtK/17N0LJUtCWJi2tO3+/aYdO+d7Hfjdd5T38iLNMHmk\nu15v1uSRef2fFHey9onjsHmmFo7FnPpie5UgTJWZman6Pv989pW+KVekl25eUpP+mKRCvwpVdSbX\nUeM2j1NnU86qK1eU+u9/lfLz09pC/vjDvFhmjB2r+nXqpDIzM83uAOEIdfiOylYdSwqCNI6L4qgw\nAxEd5UeblzVLlqhupUqpgaVLm53c9Hq92pK0RfVZ3kdVGF9BdVrUSa06tEpdT8lQU6cqVbOmUv/6\nl1KrV5s2mLAwJ397jMp3JI7QW8pUSOIQxZEj1RcXVc4T7gxQr1atqjrUr6+i+/YtVHK7nnpdxeyI\nUc1jmquAiQHq/zb+nzp04bhasECpRo2UathQqXnzlLpzJ/9YitqV2Zn/TwqjoGTrSIkFaeMQxZE1\n6ouVlermC9rvK337Zre/9AX8PT15+6OP+Hj6dKO9tgran1dJLyLDI9kauZVVPVZxLfUaLWc3ZZ5q\nx/sLlzB2wh1mz4Y6dWDiREhJMR5LVlvQgNGjOX38+H3Hyx1Hca3DL6htKosrTDYpiUM4PUs3Xlvr\nh13Qfs094ZoTZ6MqjZj89GSShiTxWuPX+Hr7NPrsCyR8xLtMnHeAbdvggQdg+HA4fdp4LLu2buWc\nkePFx8WxYeJE1i5dmv2crTsUOAJjyTZnxwtTE4uwHbsV9YRzMlZdYK26eXP2a0r7i6XiPHTpkBqx\nfoSq+llV1WZWG/XZ+rnqzUE3lY+PUj17KvV/A2NVfFycmjd9uno0IEC9VanSPccbYJh+5dWqVdVg\nQ9VacWrLMCa/jheOWIWHtHEIYTpj9dDW+mEXZr/51YNbOs47GXfUsv3L1NMLnlZlRpVSvZe+qd4e\nm6j8/bVp3Zcv16uVi+4/3tyvv1aPBgSoKA8PpQcV5eGhHg0IyHMVw+KgoMTvaD36KGTi8LDwCVwI\nh5Zz1O7ElBQ+iI7my//8J3vUbs7qGX1SkkXq5nNX+5iy3+xqqGbN7hvZXpj95cfT3ZPO9TtTel8m\n1WM2ceeDKywv3Ykq/1eZWpmRjBrbg4sn2xN8ZSsD6/2Kx5mD6HQ6Xu3Xj1OHD3Pyiy/QAW6ZmbTr\n3p2eBaxE6MqioqOz7xubkSCrCu+pLl2yZx1wRpI4RLHySt++VPT1ZdOwYXfroceOzf6RW+uHbcp+\nlVL0evppriYl5ZnYzNmfqXIm05gzN/ngi0QalChLaN9WHKu6geOdRhB4tRFX/hrKwm2fEPGvo+zY\n8jPtuuq4fPYsnqVKMeSBB/A4fpwrZ8+6fCN4URSUWIRt2bW4J5yLo1UXZFmzZIkaXK5cnhMzWktB\nVV/nb5xXn/z2iar7ZV0VPClEPTZyovINuKieeUapQa8vUqsXG6+acaRup8I4pDuuEKZxtB4/OXvb\nfHHjBn8vXMiV06fpXbWqTbqyFtSby6+sH++1fo8DAw7wTafpBDTdSeaAIM7/qzsrT1Vk0PudmTRJ\nR/PHu97TbdgVup0K41ylTGlInkI4H5VrVtqu3t40jYxk+CefsP6HHyw6+25ecs7Ka8qMv8mpyXz7\n17fEJMZw4do1/E734fiyN+j6lD81/H7gzxUf0CTD/rO/ivwVdnZcSRzC4SnD8qXvjRvnsvXnWdO8\n557+3dEppdhxdgexibF8v3cx1dIf4cqGPpQ52ZrgaxNZfPUrxgd68+jEidmTTToTV//uFTZxmFJV\nNQjwMXfHQlhKYas8lBPNzupo1Wem0ul0NK3elOkdp3NqaBLvdexM7V4TuPTag/z50CH8fDex/Ox4\n9v1dCWe8TpXqtsL7GG3p1sVAexzzf99ebUvCioo60M3VZmfN2djs6A3Po8YMUg3fqqN8x1dS1f/d\nQlVoNVPVqZeqPv5YqZMn7R1dwZx9okZTvx9YeQCgG1rSWISWRMYCdax5QDPZ6L9D2FJhB7o5+48+\nLzkToaMnxaz4Viz+Ti3as0g9MfcJVeHjSqrhu++o8nX2qrZtlZo7V6mUFHtHapwjjvI2h6nfD6zc\nq0oPnAPOA5loVVdxwKeFOagQpijsZHkFzRnkbHLPcbSwRw/+06MH4Q4431HuWLe8/x/md/+I12+/\nyI43t/Js+zKU6f8kJ59sxecbZ+Nf6yY9e8Lq1ZCebu/o73LWiRptNR+WKYljMLAD+AT4H9AAeBN4\nCOhi0WiEyKUwdf/2/tErC7et5E6EPl5e1C1fnp44XlLML2nX9qnNx098TNLQk3z6/AhqtvsBt3cD\nOdGgPyOmbKdadUX//rBxI2Rk2PudOGe7k60umkwZOe6LliD+yfW8HnjWotEIkYspI22VkZ4vRR1Z\nbWyfpspvupDCyJ0IM44fJ10phllwWhRLMWU6FA83D56r9xzP1XuO09dPM2fXHL4p3Q2/DuU5lxLF\n0JGvcPZ4Bbp0gS5dICICPD1t/16ccZS3paejcXU2qDUUjsoa9f2OtvJdzsnzhvXood595RWHXcGw\nMCssZuoz1boj61S3Jd2U9zhv1Wnuq6rvx7+qZs31ytdXqVdfVWrJEqWuX897H47eYcBWzPn8kdlx\nRXFjjRO1rHxnfxdvXlQTf5+o6k+tr+p+WVe9v+oTNW7KedWunVLlyin11FNKTZqk1MGD9y5/6+gd\nBhwRMuWIKG6sUZ9blH3au23FUpSdx79UKlOJIS2HsO+tfcx+fjan0/cz/mZdyvV+gXm/r6VPVCZ7\n9sBjj0FQEDzx6D4eCnyHdcMnyAJJwLFj8MUXcOeO9Y4hiUM4LWucqIu6T2dsUM3N1oPe8kpUOp2O\nVoGtmP38bP555x/a1m7Lx3+8z79P1aHGqx/xx74kfvgBOjwXglvl4Xx1fCNN2c76s+/RtONMnn7B\nMToMWFtmJmzbBqNHQ+PG0LIl/P33vUsBW5pzXQrlTdnr6kjYl7lzLNlrn84g5/TqtpxjKmu6lXaz\nZvHX9u0FdkjYeXYnsYmxLNq3iIcDHqZPWB/WRM+h5M+/c7HS0ySdCSbtwX4cOlaFwEBo1Uo7mTZr\nBiEh4OHki0koBYcPw6+/aj3Q1q+HqlWhfXvo1El7r+7upu1L5qqSxCFEkahcky1GBwZadY6p3Ilq\nSLVq/Hb+PE9GRjKugGompRQfR79LwEsN+Wz9BI5cOMDDJZozc8h8jv+6h6TDh3n93RHs2QO//w5b\ntsCff8KpU9pVeZMm2r+NGmnJpHx5i789i7l2DRITtVLFtm3a+/H0hEcf1W7t2kFgYOH2XdjE4eS5\nVwhhKbbuypm1qNaafv3QAWkXL/JfvZ7NGzfSMTQ035LO2qVLOTfpK9Z/W5nnS5em51V44/FDNPw8\nhDrl6/D+M/8hXd0mLKw0YWEwYID2umvXYOdO2LUL/vgDpk+HAwfA1xfq14e6daFOHa3tpFYt7YRc\noQJYu5lKr4dz5+DIETh0SCtR7NsHe/bA5ctakmvRAl58ET77TIvNnk1nkjiEENlsubRpVlJyy8gg\nyt8f9zNncANUWhoDx40zOnYiZynly7Q0hmRksO74cepnwGN/ePF+12mkhnrwzc5vGBw/mB4NexAV\nHkXDKg0B8PbWxoVERNzdp14P//yjtQscOaLdfv4ZTpyApCTt7/7+UKUK+PlptwoVtJu3N5QpA6VL\nazcPD3Bz026ZmVoD9Z07cPu21uaQkqIlr0uX4OJFuHBBO8bZs9r+goIgOFi79emjlYgeeEDbnyOR\nqiohRL6UYTDku2PH8tnIkRadYjyrPUmv1zOnVy/SvLyolZqa57TyuavT3qlYkbM3b1L9gQdQp07d\n87oTySeYtXMWs3bOwr+8P1HhUXQP7Y5XSS+zYrx2DU6f1k7y589rJ/zkZO355GQtKWTdMjO1RJOZ\nqSURT08oUQJKlQIvL+3m7Q2VK2s3Pz8ICNASU8mSpsWjLDjVe2GrqlyFHXpAC1E8ZI2PGDtsmNXG\nSZgzaC3n0r9vliihxr37br6vy8jMUCsPrlSdFnVSFcZXUH2W91FbkrY47fgaS45XoZDjOFwl0xg+\nAyGEpWRVC5W6eJHUixfx9/Rkeno6g/z8OF6pkt1W9CtKr7dzN84xZ9ccYhNjKeNZhsjwSHo26olv\naV8rR1101uj1JiUOIYRFZY2EHx4QoFaD6u/urhSo4RYYEW/v6UEy9Zlq47GN6uW4l5X3OG/VY2kP\ntfHYRpWpz7RLPKawxswEyMhxIYQlZTVep127xpyAAMjMpG9AAKkWGGiZe5ChsvFodTedG4898BgL\nuy7k6KCjtPBvwaD4QdT9si7jfxvPuRvnbBKHORxpZgJXKaIoW33hhChOsqqF/jl0iCtXrlDR15ca\ndesWelBkXtUtDVu3JnXRItrbca11pRTbTm8jJjGGpX8vJaJWBFHhUbSr0w53NxNH1FmZpQenOtsA\nQF/ge6AmcALoBiQb2e4EcB1t8ah0oHke+5PEIYQTULl6RXXx8eFK2bK0Kl2aj204Wr0gKWkpfL/v\ne2ISYziTcoY3mrxB77De1KpQy24xWUNhE4e9qqpGAOuBusDPhsfGKCACCCPvpCGEcBK5q1sC09Np\n160bysFWa/Qq6UVkeCRbI7eyqscqklOTaTqzKe0WtGPJviXcybTiDIJOwF6J4zlgruH+XKBTPtu6\nSnWaEIJ7J4LsMGcOl8+dc4h6+7w0qtKIKU9PIWlIEr0a9eKrP78iYGIA7657lwOXDtg7PLuw1//O\nVbR1y7NiuJLjcU7HgGtoVVUzgJg89idVVUI4KWecVPLw5cN8s/Mb5uyaQ92KdYkMj+SFkBco41nG\n3qGZxRHbONYDVY08/z5aKSNnoriC1u6RWzXgLFDZsL+3gc1GtpPEIYQoMmXmqOz0zHRWHlpJTGIM\nW05t4aUGLxEVHkVYtTAbRFt0jjjJ4ZP5/O08WlI5h5YcLuSx3VnDvxeBH9DaOYwlDkaPHp19PyIi\ngoick9EIIYQJzF0v3tPdk871O9O5fmdOXjvJnF1z6PR9JyqXqUxkeCQ9GvagfEnHmXo3ISGBhISE\nIu/HXlVVnwCXgQloDeMVuL+BvAzgDqQAZYF1wIeGf3OTEocQTsjcK3xrseSo7Ex9JuuPrSc2MZYN\nxzbQuX5nIsMiaRXYyqHabsD5elWNRyuRHAIeNzwGqA6sMtyvila62AVsBVZiPGkIIZyUrVcbzIsl\nlyF2d3OnfVB74rrFcejtQ4RUCqH3it6ETgtl4h8TuXTrUoH7sPWASHPZK3FcAdqidcd9irtjOM4A\nzxjuHwOaGG4NgHE2jlEIYSULZsygY2gom0eOdIh1wq01KtuvrB/vtX6PAwMOMKPjDHad20XQlCBe\ninuJDcc2oFd6o69zlISaF5lyRAhhc5a8wrcUa64Xr9PpaFOzDfM6z+P44OO0qdGGd9e9S9CUIMZs\nGsPp66cBx0uoeXGsCrfCkzYOIZxM1lrjusBA9ElJea7B4aqUUmw/s53YxFgW719Mmxpt6BPWB/e/\nbvO/d/9tk+V7na2NQwhRzFniCt/R2wLyo9PpaObfjBnPziBpSBKdHuzEhP9N4PXDb/Fz6Hl6Nwty\nyAGRICUOIYQTyyq12HNyREv78ON32O19hM23txKgq0aztFC+HD6Hkh4mLhFoBkccAGhLkjiEKEas\nsaiRo0nLSOPHAz8SkxjD7vO76dmwJ5HhkYT6hVrsGJI4JHEIUWzknmXX2m0B9nbs6jFm7ZzF7F2z\nqeldk6jwKLqFdqNsibJF2q+0cQghio3CdJ915vaQ2j61GfP4GP555x+iH4nmhwM/EPhFIP1X9mf7\nme02f0+SOIQQ+XLUE665jeuOPjbCFB5uHjxb71lWvLyCPW/uIaB8AN2WdCN8ZjhfbfuK5FRjyxpZ\nnquU6aSqSggrcfYGaFdvD9ErPb8c/4WYxBjij8TzXL3niAyPpE2NNgVW20kbhyQOISzKVU64xak9\n5NKtS8zfPZ85u+ewtudaqpYzNkH5XY44O64Qwom90rcvFX192TRs2N3R3WPHOl2pI3d7iD4pySHH\nRlhCpTKVGNJyCENaDrHqcaSNQwhhlLXmb7IHa04nUhw53zfAOKmqEsIKcq/Od/LQIa5eu2b3adCF\nZUh3XCGExUVFR2e3BbTr2pXA4GCn7plUlB5ijtq7zB4kcQghCmSJWVsd4cRblC65rtCdV9xLCSGs\nR6/Xq9WLF6sRgYFKgRoRGKjWLFmi9Hq9yftYs2SJesfLS8XHxVkxUuPmT5+ungkJUSODg5Ue1Mjg\nYPVMSIiaP326VV/r6IBCZXEpcQghClSUhnJHWGOiKOt/OOLaIfYmiUMIYZLC9kxyhBNvURJffq9V\nDlD9Zg8yjkMIYZKo6Ojs++aM5XCUcRRZiS+rh5g5XXLzem12u0ezZk43vqUoXKU/nSpuGV8IZ5K7\nW2/S4cNEjhhh77AKzVVG1cuUI5I4hBA2olxkGhMZxyGEcHmO0qbgSqPqC0MShxDCaTjSWIriPI2J\nq6RHqaoSwoW5SpuCo5HZcYUQLstVZup1FVJVJYRweMW9TcHRSIlDCOEUijIOQ1iWq6RraeMQQggz\nSXdcIYQQNiGJQwghhFkkcQghhDCLJA4hhBBmkcQhhBDCLJI4hBBCmEUShxBCCLNI4hBCCGEWeyWO\nF4F9QCYQns927YEDwGFguA3iEkIIUQB7JY49QGdgUz7buANT0ZJHCPAyUN/6oQkhhMiPveaqOmDC\nNs2BI8AJw+NFwPPA31aKSQghhAkcuY3DH0jK8fiU4TkhhBB2ZM0Sx3qgqpHnRwI/mfB6s2YtHD16\ndPb9iIgIIiIizHm5EEK4vISEBBISEoq8H3vPjvsLMAxINPK3h4HRaG0cANGAHphgZFuZHVcIIczk\nzLPj5hX0diAYqAWUALoDK2wUkxBCiDzYK3F0Rmu/eBhYBawxPF/d8BggAxgIrAX2A98jDeNCCGF3\n9q6qshSpqhJCCDM5c1WVEEIIJyKJQwghhFkkcQghhDCLJA4hhBBmkcQhhLAqpRSfjBiBdGBxHZI4\nhBBWtXbpUs5Om8a6ZcvsHYqwEEkcQgirWDBjBh1DQ9k8ciQTU1LYFB1Nx9BQFsyYYe/QRBHZa3Zc\nIYSLe6VvXyr6+rJp2DB0gD41lYFjx9Kua1d7hyaKSEocQgir0Ol06HQ6UpOTGRoSwu3k5OznhHOT\nEocQwmqSDh+m/ezZPNWlC+uWLSPp8GF7hyQswFVSv0w5IoQQZpIpR4QQQtiEJA4hhBBmkcQhhBDC\nLJI4hBBCmEUShxBCCLNI4hBCCGEWSRxCCCHMIolDCCGEWSRxCCGEMIskDiGEEGaRxCGEEMIskjiE\nEEKYRRKHEEIIs0jiEEIIYRZJHEIIIcwiiUMIIYRZJHEIIYQwiyQOIYQQZpHEIYQQwiySOIQQQphF\nEocQQgizSOIQQghhFkkcQgghzCKJQwghhFnslTheBPYBmUB4PtudAP4CdgLbrB+WEEKIgtgrcewB\nOgObCthOARFAGNDcyjFZVEJCgr1DuI/EZDpHjEtiMo3EZH32ShwHgEMmbquzZiDW4ohfFInJdI4Y\nl8RkGonJ+hy9jUMBG4DtQJSdYxFCCAF4WHHf64GqRp4fCfxk4j5aA2eByob9HQA2WyQ6IYQQhWLv\naqBfgGFAognbjgJuAJ8b+dsRoI4F4xJCiOLgKBBk7ousWeIwVV7JqwzgDqQAZYGngA/z2NbsNy6E\nEMK5dAaSgNvAOWCN4fnqwCrD/drALsNtLxBt4xiFEEIIIYQQxZGpgwfbozWmHwaG2yAuX7QG/EPA\nOqBCHttFo8W/B1gIlHSAmCoAccDfwH7gYQeICbSqyp2Y3pnCmjEForXJ7UMrAQ+yUiymfG+nGP6+\nG22Mky0UFNcrhnj+Av4HNHKAmLI0AzKALg4SUwTa93ovkOAAMVUC4rlbu/O6DWKyuQeBumg/4rwS\nhztag3ktwBPtA6lv5bg+Af5tuD8cGG9km1rAMe4mi++B1+wcE8BcoLfhvgfg7QAxAQwFvgVWWDEe\nU2OqCjQx3C8HHMTy3ylTvrcdgNWG+y2ALRaOobBxteTu96a9DeIy9TfuDmwEVgJdHSCmCmgXHwGG\nx5UcIKbRwLgc8VzGMdq/rSK/xNESLYNmGWG4WdMBoIrhflXD49x80U44Pmj/MT8Bbe0ckzdaMrMV\nU2IC7Ye1AXgM65c4TI0ppx+BJywchynf2+lA9xyPc8ZuLeb+nnyAU1aNyPSY3gHeAmZj/cRhSkxv\nAR9ZOY6cTImpH/CV4X5tTBic7egDAAvLH63xPcspw3PWVAU4b7h/HuM/5ito3YlPAmeAZLSToz1j\negC4iPbDSgRi0Hq02TMmgC+A9wC9FWMxN6YstdCqiLZaOA5TvrfGtgnAusz9PfXhbqnIWkz9rJ4H\nvjY8Vg4QUzDaBeQvaAObX3WAmGKAULRz0m5gcEE7ddTiSFEHD1rrC5JXXO8bOb6xGOqgXQHVAq4B\nS9Dqhr+1Y0weaCW3gcCfwCS0K5L/2DGmjsAFtHrgiCLEYcmYspRDaw8ajDauyJJM/d7m7sJu7ROi\nOft/DK3as7WVYsliSkxZ32WF9plZe9yaKTF5ov3enkC7QPsDrVrvsB1jGolWhRWBdo5aDzRGGwph\nlKMmjieL+PrTaI2ZWQKxTNE5v7jOo52YzgHV0E58uTUFfkerQwRYBrSiaImjqDGdMtz+NDyOo+jV\nekWNqRXwHFp9fimgPDAP6GXHmED70S8FFqBVVVmaKd/b3NsEGJ6zJlN/T43Qrl7bA1cdIKaHgEWG\n+5WAp4F0rNdmZkpMScAltKEIt9Emem2M9RKHKTG1Aj423D8KHAfqoZWIXM4vaF8MYzzQPoBaQAls\n1zie1WNhBMYbWBuj9VoojXb1MxcYYOeYQPvy1jXcHw1McICYsjyK9ds4TIlJh5a8vrBiHKZ8b3M2\njj+MbRrHTYmrBlojrDV75JkbU06zsX6vKlNiehCtetodrcSxBwixc0wT0WbmAK2a9hRadZpLMWXw\nIGhXFwfRvsy2GDzoi/aFyN2lM3dc/+Zud9y5aFex9o6pMVqJYzdaKciavapMjSnLo1i/V5UpMT2C\n1t6yC60KbSfalbWlGfve9jPcskw1/H03+XdJt2VcsWgl6azPxhbr55jyWWWxReIwNaZ3uXsOsFa3\nbnNiqoR2cbbbEFMPG8QkhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCEv7f0L7VCa9\n8ivkAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10cf9a8d0>"
       ]
      }
     ],
     "prompt_number": 122
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Robust models to outliers"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is important to consider the outliers in your data set when building a model for prediction or classification. Obviously, one way is to identify and reject the outliers. Another way is to use a model which is more robust to outliers. In prediction analysis, using absolute error instead of absolute squared error will make your model more robust. When dealing with categorical data, using tree-based models will yield more robust results."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "MSE versus MAE"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Both mean squared error (MSE) and mean absolute error (MAE) are used in predictive modeling. MSE has nice mathematical properties which makes it easier to compute the gradient. However, MAE requires more complicated tools such as linear programming to conpute the gradient. Because of the square, large errors have relatively greater influence on MSE than do the smaller error. Therefore, MAE is more robust to outliers since it does not make use of square. On the other hand, MSE is more useful if we are concerned about large errors whose consequences are much bigger than equivalent smaller ones. MSE also correspons to maximizing the likelihood of Gaussian randomv variables."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Error metric for a binary classifier"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Binary: http://www.sergulaydore.com/evaluation-metrics-for-binary-classification/\n",
      "\n",
      "Multiclass: http://www.quora.com/Classification-Machine-Learning/What-are-some-good-error-metrics-for-multi-class-labeling\n",
      "http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.104.8244&rep=rep1&type=pdf\n",
      "(F-measure (major and minor))\n",
      "\n",
      "A micro-average is generated in a traditional manner: pool all your results into one big contingency table and calculate the F-score from that.\n",
      "\n",
      "A macro-average, on the other hand, is obtained by calculating the F-score for each class, and then averaging those to get a single number.\n",
      "\n",
      "Because the micro-average gives equal weight to each individual decision, large classes will dominate small ones. This makes it a measure of how effective your classifier is on the large classes in the collection.\n",
      "\n",
      "On the other hand, macro-averaging gives equal weight to every class, and so will give you a sense of how effective you are on the small classes.\n",
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Various ways to predict a binary response variables"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Decision Trees"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A decision tree (DT) is a flowchart-like tree structure, where each internal node denotes a test on an attribute, each branch represents an outcome of the test, and each leaf node (or terminal node) holds a class label [20]. The topmost node in a tree is the root node.\n",
      "During tree construction, attribute selection measures are used to select the attribute which best partitions the tuples into distinct classes. Three popular attribute selection measures are Information Gain, Gain Ratio, and Gini Index. When DTs are built, many of the branches may reflect noise or outliers in the training data. Tree pruning attempts to identify and remove such branches, with the goal of improving classification accuracy on unseen data."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Regression"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Linear regression (LR) is used to model continuous- valued functions. It is widely used, owing largely to its simplicity. Generalized linear models represent the theoretical foundation on which LR can be applied to the modeling of categorical response variables. Common types of generalized linear models include logistic regression (LogR) and Poisson regression. LogR models the probability of some event occurring as a linear function of a set of predictor variables. Count data frequently exhibit a Poisson distribution and are commonly modeled using Poisson regression."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Naive Bayes"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bayesian classifiers are statistical classifiers. They can predict class membership probabilities [20]. Na\u00efve Bayes (NB) probabilistic classifiers are commonly studied in machine learning. The basic idea in NB approaches is to use the joint probabilities of words and categories to estimate the probabilities of categories given a document. The na\u00efve part of NB methods is the assumption of word independence, i.e. the conditional probability of a word given a category is assumed to be independent from the conditional probabilities of other words given that category. This assumption makes the computation of the NB classifiers far more efficient than the exponential complexity of non-na\u00efve Bayes approaches because it does not use word combinations as predictors"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "SVM"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A support vector machine (SVM) is an algorithm that uses a nonlinear mapping to transform the original training data into a higher dimension. Within this new dimension, it searches for the linear optimal separating hyperplane. A hyperplane is a \u201cdecision boundary\u201d separating the tuples of one class from another. With an appropriate nonlinear mapping to a sufficiently high dimension, data from two classes can always be separated by a hyperplane. The SVM finds this hyperplane using support vectors (\u201cessential\u201d training tuples) and margins (defined by the support vectors).\n",
      "\n",
      "Although the training time of even the fastest SVMs can be extremely slow, they are highly accurate, owing to their ability to model complex nonlinear decision boundaries. They are much less prone to over fitting than other methods"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Regularization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "While optimizing the fit to the data, we also optimize the simplicity of the model. Models which fit data and also simmpler are better models. Reducing the complexity helps us aviod overfitting so that our model can generalize to different data sets. This general methodology is called regularization. Ridge regression and lasso are examples of using regularization. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "How to estimate the number of retweets?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You can run GLM using the features described in this paper: http://www-users.cs.umn.edu/~echi/papers/2010-socialcom/2010-06-25-retweetability-cameraready-v3.pdf\n",
      "\n",
      "Basically, they make use of both content and contextual features. Content features include the number of URLs and hastags and contextual features include the number of followers and followees and the age of the account. The number and frequency tweets and the number of favorited tweets seem not to matter much."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "How to predict the weather using social media?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "http://www.quora.com/How-could-you-collect-and-analyze-data-to-use-social-media-to-predict-the-weather you could grab data from the Twitter, Facebook, Flickr et al public APIs filtered by location, then run some natural language processing routines over the content to extract a relevant metric from the text features i.e. words, 2-grams. Then you can compare that time-series to rainfall data or some such using e.g. @Granger causality or even go crazy and try to learn a more complex relationship using a @Kalman filter . But, excuse me if I'm missing the point, any lead time that you find in a social media signal for predicting the weather will be coincidental. Unless you are focusing on the small subset of Twitter users that have enough money to use technology that can affect the weather"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "How to construct a feed to show relavent content for a site that involves interactions with items?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I would use a similar algorithm to EdgeRank which is used by facebook http://techcrunch.com/2010/04/22/facebook-edgerank/. Basially, I will assign a score to each object (post) and rank them in a decreasing order and show only some items which have higher scores. This score can be a function of the relationship between users, comments and likes that the object got from other users and also time."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "How would you design the people you may know feature on linkedin or Facebook?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "I would assign probability of connections based on common features such as location, skills and etc."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}