线性回归/05. 线性回归的扩展.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在现实生活中，纯粹的线性模型是很难遇到的，我们不妨再来看看下面的数据集：\n",
    "\n",
    "| x | y |\n",
    "|---|---|\n",
    "|40 |216|\n",
    "|50 |399|\n",
    "|60 |496|\n",
    "|70 |507|\n",
    "|80 |432|\n",
    "|90 |271|\n",
    "\n",
    "如果我们按照线性回归的解法，也能得到解，如下所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = [[315.33333333]\n",
      " [  1.1       ]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "X = np.matrix([\n",
    "    [1, 1, 1, 1, 1, 1], \n",
    "    [40, 50, 60, 70, 80, 90]]).T\n",
    "Y = np.matrix([216, 399, 496, 507, 432, 271]).T\n",
    "\n",
    "a = (X.T * X).I * X.T * Y\n",
    "print(\"a = {0}\".format(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "得到一元线性回归：\n",
    "\n",
    "$$\n",
    "y = 1.1x + 315.33\n",
    "$$\n",
    "\n",
    "我们画出回归线和散点数据："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFQlJREFUeJzt3X2MXNd53/HvU75pRS+1eiEpcilDMkJQTpNIlBcqXQGGIyWh5RoWYdiN0sBmDCH8I0Jq94Wp2KIvRhNEBovaMlqoYaMmlBFYtlWGEhTDjEDFaFPUMkhTEmXThBhHtbhkuJQsipK9oiny6R9zhhySy93Z3VnODs/3Awzm3jN3Vs+Ohuc399wzZyMzkSTV5+91uwBJUncYAJJUKQNAkiplAEhSpQwASaqUASBJlWorACJiICIej4gfRMS+iHh/RFwTEU9HxEvl/upybETElyLiQES8EBG3zeyvIEmainbPAB4CvpmZNwO3APuAB4CdmbkS2Fn2Ae4GVpbbBuDhjlYsSeqImOiLYBGxCHgeeE+2HBwR+4EPZubhiFgGfCszV0XEH5Xtr5x/3Iz9FpKkSZvbxjHvAY4CfxIRtwC7gc8AS5udegmBJeX4QeCVlucfLG3nBEBEbKBxhsDChQvfd/PNN0/n95Ck6uzevfvVzFw81ee3EwBzgduA383MZyPiIc4O94wlxmi74DQjM7cAWwCGhoZy165dbZQiSWqKiP83nee3cw3gIHAwM58t+4/TCIQjZeiHcj/ScvwNLc9fARyaTpGSpM6bMAAy8++AVyJiVWm6C/g+8CSwvrStB54o208CnyqzgdYAbzj+L0mzTztDQAC/C/xZRMwHfgh8mkZ4fC0i7gN+BHyiHPsN4MPAAeCn5VhJ0izTVgBk5nPA0BgP3TXGsQncP826JEkzzG8CS1KlDABJqpQBIEmVMgAkqVIGgCRVygCQpEoZAJJUKQNAkiplAEhSpQwASaqUASBJlTIAJKlSBoAkVcoAkKRKGQCSVCkDQJIqZQBIUqUMAEmqlAEgSZUyACSpUgaAJFXKAJCkShkAklQpA0CSKmUASFKlDABJqpQBIEmVMgAkqVIGgCRVygCQpEoZAJJUKQNAkio1t9sFSDXavmeYzTv2c+jYKMsH+ti4dhXrVg92uyxVxgCQLrHte4bZtG0voydPATB8bJRN2/YCGAK6pBwCki6xzTv2n+n8m0ZPnmLzjv1dqki1aisAIuLliNgbEc9FxK7Sdk1EPB0RL5X7q0t7RMSXIuJARLwQEbfN5C8g9ZpDx0Yn1S7NlMmcAfxyZt6amUNl/wFgZ2auBHaWfYC7gZXltgF4uFPFSpeD5QN9k2qXZsp0hoDuAbaW7a3Aupb2R7Ph28BARCybxn9HGtP2PcPc8eAz3PTAX3DHg8+wfc9wt0tqy8a1q+ibN+ectr55c9i4dlWXKlKt2g2ABP4yInZHxIbStjQzDwOU+yWlfRB4peW5B0vbOSJiQ0TsiohdR48enVr1qlbzQurwsVGSsxdSeyEE1q0e5A8/9osMDvQRwOBAH3/4sV/0ArAuuXZnAd2RmYciYgnwdET8YJxjY4y2vKAhcwuwBWBoaOiCx6XxjHchtRc60nWrB3uiTl3e2joDyMxD5X4E+HPgduBIc2in3I+Uww8CN7Q8fQVwqFMFS+CFVKkTJgyAiFgYEf3NbeDXgBeBJ4H15bD1wBNl+0ngU2U20BrgjeZQkdQpXkiVpq+dM4ClwF9HxPPAd4C/yMxvAg8CvxoRLwG/WvYBvgH8EDgA/HfgdzpetarnhVRp+ia8BpCZPwRuGaP9NeCuMdoTuL8j1UkX0Rw/dzkFaepcCkI9ywup0vS4FIQkVcoAkKRKGQCSVCkDQJIqZQBIUqUMAEmqlAEgSZUyACSpUgaAJFXKAJCkShkAklQpA0CSKmUASFKlDABJqpQBIEmVMgAkqVIGgCRVygCQpEoZAJJUKf8msKRJ275nmM079nPo2CjLB/rYuHaVf5+5BxkAkiZl+55hNm3by+jJUwAMHxtl07a9AIZAj3EISNKkbN6x/0zn3zR68hSbd+zvUkWaKgNA0qQcOjY6qXbNXgaApElZPtA3qXbNXl4DkDQpG9euOucaAEDfvDlsXLuqi1Vd/k6dTl576wRHjp9g5M23OXL8xLR/pgEgaVKaF3qdBdQZY3XsZ+6Pv83Imyc4cvxtXn3rBKezs//tyOzwT5yCoaGh3LVrV7fLkKSOmW7Hft275rO4/wqW9C/g+kVXsGTRApYsauwvLfeDV1+5OzOHplqjZwCSNAnNjr3Zgbd27EfLfTsd+9JFC/j5ZYsaHXt/o3NfuqjRft27FjBvzsxfojUAJInOduzvXdbf+JR+3if2xf2XpmNvlwEg6bJ2fsfe2sEfbRmaOfpmex37kv4rWHrV7O7Y22UASOpJF+vYR95sjLE3O/ZX3/oZp8bo2a9dOJ/F/Qu4/qorzn5i7z93nP26dy1g/tze69jbZQBImlXaGYoZ7xP7tQvnn+nEz3xiX3R2jH1J/4LLvmNvlwEg6ZIYb1bMZMbYmx17jZ/YO63tAIiIOcAuYDgzPxIRNwGPAdcA3wU+mZk/i4gFwKPA+4DXgF/PzJc7XrmkWaFT0x3PmRXT7Nzt2GfUZM4APgPsAxaV/c8DX8jMxyLivwH3AQ+X+9cz8+ci4t5y3K93sGZ1iEv6ajzjDcW007Ffu3B+mQlzYcd+qac7amxtBUBErAD+EfAHwD+PiADuBP5JOWQr8B9oBMA9ZRvgceC/RETkbPjGmc5wSd96zfR0Rzv23tHuGcAXgd8D+sv+tcCxzHyn7B8Emr3GIPAKQGa+ExFvlONfbf2BEbEB2ADw7ne/e6r1a4rGW9LXAOhNduyarAkDICI+Aoxk5u6I+GCzeYxDs43HzjZkbgG2QGMpiLaqVce4pG/v6NQ89nMunp73BSXH2OvUzhnAHcBHI+LDwBU0rgF8ERiIiLnlLGAFcKgcfxC4ATgYEXOBq4Afd7xyTcvygT6Gx+jsXdL30jl1OnntJycYOT7+PPaJpjue84m93+mOat+EAZCZm4BNAOUM4F9m5m9GxNeBj9OYCbQeeKI85cmy/3/L4884/j/7uKTvzGm3Yx/vC0qt89jt2DVTpvM9gH8FPBYRvw/sAR4p7Y8AX46IAzQ++d87vRI1E1zSd/JaO/aRljH1TnTsTndUN7gctKo3Ucd+drrj+EsKNC+Utn7ztNnBL7Zj1wyICJeDlsYyVsc+cvwER958u+2OvdmJ33x9/wUdu5/Y1esMAPWciaY7TuYTe2vHvrilg/cTu2pgAGjWmO489taOfdXSfpYsavwlpWbH7id26VwGgGbcdOexX7Nw/pkhl5uv7z8z1r64v7HMwFI/sUtTYgBoymZyHrtj7NLMMwB0ATt2qQ4GQEU6OY/95uv7uf4qO3aplxkAl4FOzGMf6y8oefFUurwZALOY89glzSQDoAvG79ibbVOfx27HLqkdBkAHdapjX2LHLukSMADacPp08mpLx96YHTO1T+yrljbXYz93zRjnsUu61KoOgNOnk9d+8rNywXTsjn3k+AmOvnVizI699QtKduySes1lGQAX69hHWr51ascuqXY9FQDjT3dsv2NfMkbH3rhvbNuxS6rBrAmAkXPG0yc3j/38jr25PswS14qRpIuaFQGwd/gNbv+DnRe0j3fx1I5dkqZnVgTAkv4F/Md7/r7THSXpEpoVAbB00RV88v03drsMSaqKH7ElqVIGgCRVygCQpEoZAJJUKQNAkiplAEhSpQwASaqUASBJlTIAJKlSBoAkVcoAkKRKGQCSVCkDQJIqZQBIUqUMAEmqlAEgSZWaMAAi4oqI+E5EPB8R34uIz5X2myLi2Yh4KSK+GhHzS/uCsn+gPH7jzP4KkqSpaOcM4ARwZ2beAtwKfCgi1gCfB76QmSuB14H7yvH3Aa9n5s8BXyjHSZJmmQkDIBveKrvzyi2BO4HHS/tWYF3ZvqfsUx6/KyKiYxVLkjqirWsAETEnIp4DRoCngb8BjmXmO+WQg8Bg2R4EXgEoj78BXDvGz9wQEbsiYtfRo0en91tIkiatrQDIzFOZeSuwArgdeO9Yh5X7sT7t5wUNmVsycygzhxYvXtxuvZKkDpnULKDMPAZ8C1gDDETE3PLQCuBQ2T4I3ABQHr8K+HEnipUkdU47s4AWR8RA2e4DfgXYB/wV8PFy2HrgibL9ZNmnPP5MZl5wBiBJ6q65Ex/CMmBrRMyhERhfy8ynIuL7wGMR8fvAHuCRcvwjwJcj4gCNT/73zkDdkqRpmjAAMvMFYPUY7T+kcT3g/Pa3gU90pDpJ0ozxm8CSVCkDQJIqZQBIUqXauQgsSZeV7XuG2bxjP4eOjbJ8oI+Na1exbvXgxE+8zBgAkqqyfc8wm7btZfTkKQCGj42yadtegOpCwCEgSVXZvGP/mc6/afTkKTbv2N+lirrHAJBUlUPHRifVfjkzACRVZflA36TaL2cGgKSqbFy7ir55c85p65s3h41rV3Wpou7xIrCkqjQv9DoLyACQVKF1qwer7PDP5xCQJFXKAJCkShkAklQpA0CSKmUASFKlDABJqpQBIEmVMgAkqVIGgCRVygCQpEoZAJJUKQNAkiplAEhSpQwASaqUASBJlTIAJKlS/kGYadi+Z9i/KiSpZxkAU7R9zzCbtu1l9OQpAIaPjbJp214AQ0BST3AIaIo279h/pvNvGj15is079nepIkmaHANgig4dG51UuyTNNgbAFC0f6JtUuyTNNgbAFG1cu4q+eXPOaeubN4eNa1d1qSJJmhwvAk9R80Kvs4Ak9aoJAyAibgAeBa4HTgNbMvOhiLgG+CpwI/Ay8I8z8/WICOAh4MPAT4Hfyszvzkz53bVu9aAdvqSe1c4Q0DvAv8jM9wJrgPsj4ueBB4CdmbkS2Fn2Ae4GVpbbBuDhjlctSZq2CQMgMw83P8Fn5pvAPmAQuAfYWg7bCqwr2/cAj2bDt4GBiFjW8colSdMyqYvAEXEjsBp4FliamYehERLAknLYIPBKy9MOlrbzf9aGiNgVEbuOHj06+colSdPSdgBExLuA/wl8NjOPj3foGG15QUPmlswcysyhxYsXt1uGJKlD2gqAiJhHo/P/s8zcVpqPNId2yv1IaT8I3NDy9BXAoc6UK0nqlAkDoMzqeQTYl5n/ueWhJ4H1ZXs98ERL+6eiYQ3wRnOoSJI0e7TzPYA7gE8CeyPiudL2r4EHga9FxH3Aj4BPlMe+QWMK6AEa00A/3dGKJUkdMWEAZOZfM/a4PsBdYxyfwP3TrEuSNMNcCkKSKmUASFKlDABJqpQBIEmVMgAkqVIGgCRVygCQpEoZAJJUKQNAkiplAEhSpQwASaqUASBJlTIAJKlSBoAkVcoAkKRKGQCSVCkDQJIqZQBIUqUMAEmqlAEgSZUyACSpUgaAJFXKAJCkShkAklQpA0CSKmUASFKlDABJqpQBIEmVMgAkqVIGgCRVygCQpEoZAJJUKQNAkiplAEhSpQwASarUhAEQEf8jIkYi4sWWtmsi4umIeKncX13aIyK+FBEHIuKFiLhtJouXJE1dO2cAfwp86Ly2B4CdmbkS2Fn2Ae4GVpbbBuDhzpQpSeq0CQMgM/8X8OPzmu8BtpbtrcC6lvZHs+HbwEBELOtUsZKkzpk7xectzczDAJl5OCKWlPZB4JWW4w6WtsPn/4CI2EDjLAHgROsQUw+6Dni120VMg/V3Ty/XDtbfbaum8+SpBsDFxBhtOdaBmbkF2AIQEbsyc6jDtVwy1t9dvVx/L9cO1t9tEbFrOs+f6iygI82hnXI/UtoPAje0HLcCODT18iRJM2WqAfAksL5srweeaGn/VJkNtAZ4ozlUJEmaXSYcAoqIrwAfBK6LiIPAvwceBL4WEfcBPwI+UQ7/BvBh4ADwU+DTbdaxZXJlzzrW3129XH8v1w7W323Tqj8yxxyilyRd5vwmsCRVygCQpEpd8gCIiCsi4jsR8XxEfC8iPlfab4qIZ8vyEl+NiPmXurZ2RcSciNgTEU+V/V6q/eWI2BsRzzWnkF1saY/ZKCIGIuLxiPhBROyLiPf3Sv0Rsaq87s3b8Yj4bK/UDxAR/6z8u30xIr5S/j33xPs/Ij5T6v5eRHy2tM3q136ml+LpxhnACeDOzLwFuBX4UJkx9HngC2V5ideB+7pQW7s+A+xr2e+l2gF+OTNvbZn/fLGlPWajh4BvZubNwC00/j/0RP2Zub+87rcC76MxUeLP6ZH6I2IQ+KfAUGb+AjAHuJceeP9HxC8Avw3cTuN985GIWMnsf+3/lJlciiczu3YDrgS+C/wDGt/Gm1va3w/s6GZt49S8orzodwJP0fjyW0/UXup7GbjuvLb9wLKyvQzY3+06L1L7IuBvKZMXeq3+82r+NeD/9FL9nP2m/zU0ZhA+Bazthfc/jZmKf9yy/2+B3+uF1x64EXixZX/MmoE/An5jrOMuduvKNYAyhPIcjS+QPQ38DXAsM98phzSXkJiNvkjjjXO67F9L79QOjW9m/2VE7C7LccB5S3sASy767O56D3AU+JMyBPfHEbGQ3qm/1b3AV8p2T9SfmcPAf6Ix9fsw8Aawm954/78IfCAiro2IK2lMV7+BHnntz3Oxmi+2FM9FdSUAMvNUNk6DV9A4JXvvWIdd2qomFhEfAUYyc3dr8xiHzrraW9yRmbfROF28PyI+0O2CJmEucBvwcGauBn7C7Dtln1AZI/8o8PVu1zIZZaz5HuAmYDmwkMb76Hyz7v2fmftoDFU9DXwTeB54Z9wn9Z5J90VdnQWUmceAbwFraKwc2vxi2mxdQuIO4KMR8TLwGI1hoC/SG7UDkJmHyv0IjfHn27n40h6zzUHgYGY+W/YfpxEIvVJ/093AdzPzSNnvlfp/BfjbzDyamSeBbcA/pEfe/5n5SGbelpkfoLHC8Uv0zmvfqmNL8XRjFtDiiBgo23003lT7gL8CPl4Oa11eYtbIzE2ZuSIzb6RxCv9MZv4mPVA7QEQsjIj+5jaNcegXufjSHrNKZv4d8EpENFdAvAv4Pj1Sf4vf4OzwD/RO/T8C1kTElRERnH39e+X9v6Tcvxv4GI3/B73y2rfq3FI8Xbig8UvAHuAFGp3Pvyvt7wG+Q2MZia8DC7p98WWC3+ODwFO9VHup8/ly+x7wb0r7tTQubL9U7q/pdq3j/A63ArvK+2c7cHWP1X8l8BpwVUtbL9X/OeAH5d/ul4EFPfT+/980Aut54K5eeO1phNRh4CSNT/j3XaxmGkNA/5XGNdW9NGZrjfvzXQpCkirlN4ElqVIGgCRVygCQpEoZAJJUKQNAkiplAEhSpQwASarU/wdsjqnlqVs8kwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.xlim(30, 100)\n",
    "plt.ylim(100, 600)\n",
    "\n",
    "X = np.array([40, 50, 60, 70, 80, 90])\n",
    "Y = np.array([216, 399, 496, 507, 432, 271])\n",
    "plt.scatter(X, Y)\n",
    "\n",
    "x = np.linspace(30, 100)\n",
    "y = 1.1*x+315.33\n",
    "plt.plot(x, y)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从图中很容易看出，线性回归对数据的拟合度不是很高，很显然散点的分布不是直线，而是一个二次曲线。我们知道二次曲线通常可以写成这样：\n",
    "\n",
    "$$\n",
    "y = a_0 + a_1x + a_2x^2\n",
    "$$\n",
    "\n",
    "到现在为止，我们一直在处理线性回归问题，那么对于这种不是线性的情况，该如何求解呢？首先，看到上面这个式子，我们不由得联想到了我们之前遇到的二元线性回归：\n",
    "\n",
    "$$\n",
    "y = a_0 + a_1x_1 + a_2x_2\n",
    "$$\n",
    "\n",
    "二元线性回归问题处理的是两个自变量，但是我们这里只有一个自变量，怎么办呢？答案是：自己造一个出来。\n",
    "\n",
    "| x |$x^2$| y |\n",
    "|---|-----|---|\n",
    "|40 |1600 |216|\n",
    "|50 |2500 |399|\n",
    "|60 |3600 |496|\n",
    "|70 |4900 |507|\n",
    "|80 |6400 |432|\n",
    "|90 |8100 |271|\n",
    "\n",
    "我们引入一个新变量 $x^2$，令：$x_1 = x, x_2 = x^2$，这样我们就可以像处理线性回归问题一样处理非线性回归问题了。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = [[-1376.  ]\n",
      " [   57.  ]\n",
      " [   -0.43]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "X = np.matrix([\n",
    "    [1, 1, 1, 1, 1, 1], \n",
    "    [40, 50, 60, 70, 80, 90],\n",
    "    [1600, 2500, 3600, 4900, 6400, 8100]]).T\n",
    "Y = np.matrix([216, 399, 496, 507, 432, 271]).T\n",
    "\n",
    "a = (X.T * X).I * X.T * Y\n",
    "print(\"a = {0}\".format(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "得到回归问题的解：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y &= 57x_1 - 0.43x_2 - 1376 \\\\\n",
    "&= 57x - 0.43x^2 - 1376\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "画出二次曲线图，拟合度非常高："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VeXhx/HPkz0YSYAwMggbZENAljhwo0LdVgUVRds6qpaK7a/DakvVqlWrVpxIq3UUBSciIMgmyN4BAhmQhJGwkpDx/P64BxsxmJ1zb+73/XrllXvPPTd+E673e885z3mOsdYiIiL+J8DtACIi4g4VgIiIn1IBiIj4KRWAiIifUgGIiPgpFYCIiJ+qUgEYY6KMMR8YY7YYYzYbY4YaY2KMMXOMMdud79HOusYY85wxJtUYs84YM6B+fwUREamJqm4BPAt8Ya3tDvQFNgOTgbnW2i7AXOc+wCVAF+drIvBSnSYWEZE6YSo7EcwY0wxYC3S05VY2xmwFzrHW7jXGtAW+ttZ2M8a87Nx+59T16u23EBGRaguqwjodgVzgDWNMX2AVcB/Q+uSbulMCsc76cUB6uednOMu+VwDGmIl4thCIjIwc2L1799r8HiIifmfVqlX7rbWtavr8qhRAEDAAuMdau9wY8yz/291TEVPBsh9sZlhrpwJTAZKTk21KSkoVooiIyEnGmN21eX5VjgFkABnW2uXO/Q/wFEK2s+sH53tOufUTyj0/HsiqTUgREal7lRaAtXYfkG6M6eYsGgVsAmYB451l44GZzu1ZwDhnNNAQIF/7/0VEvE9VdgEB3AP82xgTAuwEbsVTHu8ZYyYAe4BrnHU/Ay4FUoHjzroiIuJlqlQA1to1QHIFD42qYF0L/KKWuUREpJ7pTGARET+lAhAR8VMqABERP6UCEBHxUyoAERE/pQIQEfFTKgARET+lAhAR8VMqABERP6UCEBHxUyoAERE/pQIQEfFTKgARET+lAhAR8VMqABERP6UCEBHxUyoAERE/pQIQEfFTKgARET+lAhAR8VMqABERP6UCEBHxUyoAERE/pQIQEfFTKgARET+lAhAR8VMqABERP6UCEBHxUyoAERE/pQIQEfFTKgARET8V5HYAEW91/EQJ2YeLyDlcSPYRz/eikjICAwyBxhAQYAg0eO4HBNCySQjx0RHERYfTPDzY7fgilVIBiN/LP17Musw81uzJY21GHjv3HyPncBFHi0pq/DObhgURFxVOfHQECTHh9EuIIjkphriocAA+Wp3Jk7O3kpVXQLuocCZd1I2x/ePq6lcSqRIVgPidHblHWbR9P2vT81iT7nnDP6lTq0i6t2nKyC6taN0sjNbNQolt+r/vocEBlFlLaZmlrAxKndslZWXkHC4iM6+AjEPHyTxUQMahAtIPHmdx6n7eWJwGQLvmYbRpHsb6jHyKyywAmXkFPDxjPYBKQBqUCkAaPWstm/ce4YuN+/hiw162ZR8FoFXTUPolRHHVwHj6JUTRK655rXbdtG0eTt+EqB8sLy2zbNl3mJS0Q6xMO8hn6/fivPd/p6C4lCe+2KICkAZlrLWVr2RMGnAEKAVKrLXJxpgY4F0gCUgDrrXWHjLGGOBZ4FLgOHCLtfbbH/v5ycnJNiUlpRa/hsj3WWtZm5HP5xv2MnvDPtIOHCfAwKCkGC7p1Ybzz2hNXFQ4npdrw0qa/OlpH7tlWBKX923HgMQoV7KJbzHGrLLWJtf0+dXZAjjXWru/3P3JwFxr7V+NMZOd+w8BlwBdnK8zgZec7yL1rrC4lJlrMnljcRpb9h0hKMAwrHNLJo7sxIU9W9OySajbEYmLCiczr+AHy8OCA3h7xR7eXJJGfHQ4P+kfx81D2hPbLMyFlOIParMLaAxwjnN7GvA1ngIYA7xlPZsWy4wxUcaYttbavbUJKnKq8gdSWzcLo3d8c1btPsTBYyfo3qYpU67szaW92tI8wrtG5Ey6qBsPz1hPQXHpd8vCgwOZcmVvzusRy5cbs5m1NosX5qfy8oKd/KR/HHeM7Ejn2CYuppbGqKq7gHYBhwALvGytnWqMybPWRpVb55C1NtoY8wnwV2vtImf5XOAha23KKT9zIjARIDExceDu3bvr7JeSxu+j1Zk/eBMF6NWuGb8Z3YOhHVt49S6UqowCStt/jFcX7eT9lAyKSso4v0csd57dieT20V79u0nDqe0uoKoWQDtrbZYxJhaYA9wDzDpNAXwKTDmlAH5trV11up+vYwBSXYP//BU5R4p+sDwuKpzFk89zIVH9OXC0iGlLdzN9aRqHjhfTPzGKu8/tzHndY1UEfq62BVClM4GttVnO9xzgQ2AwkG2MaeuEaAvkOKtnAAnlnh4PZNU0oEh5OYcLeXjG+grf/AGyKti37utaNAnlgQu6smTyKP40picHjp5gwrQUbnptOZv3HnY7nviwSgvAGBNpjGl68jZwIbABmAWMd1YbD8x0bs8CxhmPIUC+9v9LbR0tKuHpOds4+8mveT8lnciQwArXa+ecaNUYhYcEMm5oEnMfPJtHrujJxqzDjH7uG37z4Xr2H624EEV+TFUOArcGPnQ2NYOAt621XxhjVgLvGWMmAHuAa5z1P8MzBDQVzzDQW+s8tfiNsjLL2yv28PevtrH/6AlG92nLpAu7sSY9r8IDqZMu6uZi2oYRHBjA+GFJjOnXjmfnbmf60t18vCaLe0Z1ZvywJEKDKi5HkVNV6RhAfdMxAKnIztyj/PqDdaTsPsTgpBgevrQ7/ROjv3tc0yl4pOYc5S+fbWbelhzat4jgj1f05NxusW7HkgbQIAeB65sKQMorLbO8tmgnT325jdCgAH5/eU+uGhCnA56VWLAtl0c/2URqzlGuH5TAb0f3oGmYdw2BlbrVkCeCidS77dlH+NUH61ibnsf5PVrz55/0orVOhKqSs7u2Ysi9I3hmznamLtzBN9v387dr+jK0Uwu3o4mX0vUAxCsUl5bxwvxURj+3iD0HjvHs9f14ZdxAvflXU2hQIJMv6c77dw0lONBwwyvLeOTjjRSecr6ECGgLQLxA+sHj3P32t6zNyGd077Y8MqanV0zZ4MsGto/hs/vO4vHPt/DG4jQWbMvl6Wv70a+CyerEf2kLQFw1f0sOlz2/iJ37j/HCTwfwwo0D9OZfRyJCgnhkTC/+NeFMCk+UcuWLi3lhfirecNxPvIMKQFxRWmZ56sut3PrmStpFhfPJPSMY3aet27EapRFdWvLF/SO5rE87npy9lbv+tYojhcVuxxIvoAKQBnfgaBHjX1/B8/NSuTY5ng9/Poz2LSLdjtWoNQsL5tnr+/G7y87gq805jH1hMak5R92OJS5TAUiDWrX7EKOfW8TKtIM8cVUfnri6L2HBOnGpIRhjmDCiA/+acCZ5x4sZ+8Jivty4z+1Y4iIVgDSY6ct2c93LSwkJCmDGz4dx7aCEyp8kdW5opxZ8fM8IOraKZOL0VTz15VZKT71EmfgFFYDUu7Iyy5TPN/O7jzZwdtdWfHzPCHq2a+52LL/WLiqc9+4cyrXJ8Tw/L5UJ01ZyWMcF/I4KQOpVUUkp9727hpcX7OTmIe2ZOi65VtfdlboTFhzI41f14bGxvVi0fT/XvbyMnMOFbseSBqQCkHqTX1DMuNdW8PHaLCZf0p0/jelJYICmc/AmxhhuGtKe128ZxO4Dx7jypSXszNXBYX+hApB6kZlXwDX/XMK3ew7x7PX9uOvsTprLx4uN7NqKd+4YwvETpVz9z6WsTc9zO5I0ABWA1LlNWYe58sXF7M0vZNptgxnTz/9m6PRFfROi+OCuoUSEBHLDK8tYuC3X7UhSz1QAUqeW7NjPtS8vJcAYPrhrGMM6tXQ7klRDx1ZNmPEzz3kZt725kplrMt2OJPVIBSB1ZtH2/dz6xkraRYXx4c+H061NU7cjSQ3ENgvj3TuHMLB9NPf9Zw2vfrPT7UhST1QAUie+2Z7LhGkr6dAyknfuGEKb5prF05c1Cwtm2m2DuaRXGx77dDMvfp3qdiSpByoAqbUF23KZMC2FDi0jefuOIbTQZG6NQlhwIP/46QDG9GvHE19sZerCHW5Hkjqm6aClVr7emsPE6avo1KoJb99+JtGRIW5HkjoUGGB46pq+lJRZ/vLZFgIDApgwooPbsaSOqACkxuZvzeHO6avo3KoJ/9abf6MVFBjA36/rR1mZ5dFPNhEcaBg3NMntWFIHVABSI/O3eN78u7T2vPlHRejNvzELDgzg2ev7U/L2t/x+5kbWZ+SzZMcBsvIKaBcVzqSLujG2v4b7+hodA5BqW7gtlzunr6JrG735+5OQoAD+8dP+nNG2Ge+vyiAzrwCL56S/h2es56PVGjLqa1QAUi1r0/O461+r6Ngqkn9PGKI3fz8TGhRI3vETP1heUFzKk7O3upBIakMFIFW2M/cot765kpjIEN66bTDNIzSpmz/am1/xhHFZeQUNnERqSwUgVZJzuJBxr6/AANMnnElsM43z91ftosKrtVy8lwpAKnW4sJhxr6/g4LETvHHrIDq01OUb/dmki7oRXsFV3G4ZltTwYaRWVADyowqLS7ljWgqpOUf5500D6RMf5XYkcdnY/nFMubI3cVHhGKB101AiQwOZvmw3uUeK3I4n1aBhoHJapWWW+99dw/JdB3n2+n6M7NrK7UjiJcb2j/vesM/Vew5xwyvLmDBtJe/cMYTIUL21+AJtAUiFrLX8YdYGPt+wj/8b3UNTOsuP6p8YzT9uGMCGzHzufvtbSkrL3I4kVaACkAq9tmgX/1q2hztHduT2szq6HUd8wPlntObRsb2YvzWX//toA9bqQvPeTttp8gPzt+bwl882c0mvNjx0cXe344gPufHM9uzNK+Qf81NpFxXOvaO6uB1JfoQKQL4nNeco9769mm5tmvHUtX0J0DV8pZoevLArWXkFPD1nG11bN+XiXm3cjiSnoV1A8p3848VMfCuFkKAAXhk3kIgQfT6Q6jPG8Jcre9M3IYoH31vDtuwjbkeS01ABCAAlpWXc/c63pB86zj9vHkh8dITbkcSHhQUH8vJNAwkPCeKOt1LIP17sdiSpQJULwBgTaIxZbYz5xLnfwRiz3Biz3RjzrjEmxFke6txPdR5Pqp/oUpemfL6Fb7bv57GxvRiUFON2HGkE2jQP4+WbB5CVV8A9/1lNaZkOCnub6mwB3AdsLnf/ceAZa20X4BAwwVk+AThkre0MPOOsJ17oo9WZDP/rPJImf8pri3YxsktLrhuU6HYsaUQGto/hT2N6sXBbLk/M3uJ2HDlFlQrAGBMPjAZede4b4DzgA2eVacBY5/YY5z7O46Oc9cWLfLQ6k4dnrCez3AReK3Yd1JS+UuduGJzIjWcm8vKCncxam+V2HCmnqlsAfwd+DZw8u6MFkGetLXHuZwAnzxSKA9IBnMfznfW/xxgz0RiTYoxJyc3NrWF8qaknZ2+loLj0e8sKS8o0pa/Uiz9c3pNBSdH8+oO1bMjMdzuOOCotAGPMZUCOtXZV+cUVrGqr8Nj/Flg71VqbbK1NbtVKUww0tMzTTN2rKX2lPoQEBfDijQOJCg/hzumrOHBUcwZ5g6psAQwHrjDGpAH/wbPr5+9AlDHm5DjBeODktl0GkADgPN4cOFiHmaUOND3NXC2a0lfqS6umoUwdN5Dco0Xc9581lOmgsOsqLQBr7cPW2nhrbRJwPTDPWnsjMB+42lltPDDTuT3LuY/z+Dyrc8K9ysJtuRwtKiHwlEMz4cGBTLqom0upxB/0iY/ikSt6sih1Py8t2OF2HL9Xm/MAHgIeMMak4tnH/5qz/DWghbP8AWBy7SJKXco+XMj9766ha+um35vSNy4qnClX9taFvaXeXT8ogcv7tuPpOdtISdPOATcZb/hwnpycbFNSUtyO0eiVlJbx01eWsyErn1l3j6BzbBO3I4mfOlJYzGXPL6K4pIxP7z2L6EhdW7omjDGrrLXJNX2+zgT2I898tY0VaQf580966c1fXNU0LJjnb+hP7tEiJn2wTjOHukQF4Ce+3prDC/N3cP2gBH7SP97tOCL0iY/i4Ut68NXmbN5YnOZ2HL+kAvADe/MLuP/dNXRv05Q/XtHT7Tgi37l1eBLn92jNlM83sy4jz+04fkcF0MiVllnue2cNJ0rKeOHGAYRVcDFvEbcYY/jbNX1o1SSUe95ZzZFCTRrXkFQAjdzUhTtZkXaQR8f2olMr7fcX7xMVEcJzN/Qn41ABv/lQVxJrSCqARmxjVj5Pz9nK6N5t+YmGd4oXS06K4YELuvLx2izeT8lwO47fUAE0UoXFpdz/7hqiI0J4bGwvNB+feLufnd2JoR1b8MjHG0k/eNztOH5BBdBI/W32VrZlH+WJq/tojLX4hIAAw5PX9MEYw6/eX6upIhqACqARWrJjP68t3sXNQ9pzTrdYt+OIVFl8dAS/v+wMlu86yBtL0tyO0+ipABqZw4XF/Oq9tSS1iOThS7u7HUek2q5JjmdU91ie+GILqTlH3Y7TqKkAGpk/ztxI9pEinr62ry7qLj7JGMOUq3oTERLIg++toaS0rPInSY2oABqRz9bvZcbqTO4+tzP9E6PdjiNSY7FNw3hsbG/WZuTz0teaNbS+qAAaiZzDhfzmw/X0iW/O3ed1djuOSK2N7tOWy/u249m523UVsXqiAmgErLVMnrGewuJSnrmuH8GB+meVxuHRMT2JiQzhwffWUlRSWvkTpFr0TtEIzFqbxbwtOfzqwm4621calaiIEB6/qg9bs4/w9JxtbsdpdFQAPu7gsRM88vEm+iZEcevwDm7HEalz53aP5YbBCUxduJNVuw+5HadRUQH4uEc/2cThgmIev6o3gQE621cap9+OPoO2zcJ4eMY6TpRoVFBdUQH4sK+35vDh6kx+fk4nurdp5nYckXrTJDSIR8f2Ylv2UV7WtYTrjArARx0tKuG3H26gc2wTfqFRP+IHRvVozeg+bXl+Xio7cnWCWF1QAfiov83eSlZ+AY9f1ZvQIM3xL/7hD5efQVhwAL+ZsV5zBdUBFYAPWrX7ENOWpjFuSHsGto9xO45Ig4ltGsZvLu3B8l0HeX9VuttxfJ4KwMcUlZTy0H/X0bZZGJMu1lw/4n+uTU5gcIcY/vzpZnKOFLodx6epAHzMC/N3kJpzlD9f2ZsmoZrrR/xPQIBhypW9KSwu408fb3I7jk9TAfiQrfuO8NLXqYzt145zNc2z+LFOrZpw93md+WTdXuZtyXY7js9SAfgIay3/99F6moQG8fvLe7odR8R1d53diS6xTfjdRxs5VlTidhyfpALwER+uzmRl2iEmX9KdGF3hS4SQoAD+elVvsvILeOpLTRNREyoAH5BfUMxfPttMv4QorhmY4HYcEa8xsH0MN53ZnjeX7NKMoTWgAvABz8zZxoFjJ3hsbC8CNN2DyPdMurgbMZEh/H7mBp0bUE0qAC+3Keswby1N46Yz29MrrrnbcUS8TrOwYB66uDvf7snjw9WZbsfxKSoAL1ZWZvn9zA1ERYTwqwu7uR1HxGtdNSCe/olRTPl8C4cLi92O4zNUAF5sxupMUnZ7Dvw2jwh2O46I1woIMDxyRU8OHCvi2a+2ux3HZ6gAvFR+QTFTPttM/8Qorh4Q73YcEa/XJz6K6wcl8uaSNLZlH3E7jk9QAXipZ+Zs49DxEzw6Rgd+Rapq0kXdaBIaxB9nbcRaHRCujArAC23Myvcc+B2iA78i1RETGcKvLuzKkh0H+Gz9PrfjeD0VgJfxHPjdSHRECA9eoAO/ItX10zPbc0bbZjz26SaOn9AZwj+m0gIwxoQZY1YYY9YaYzYaYx5xlncwxiw3xmw3xrxrjAlxloc691Odx5Pq91doXD5cnckqHfgVqbHAAMOfxvRkb34hL87X1cN+TFW2AIqA86y1fYF+wMXGmCHA48Az1touwCFggrP+BOCQtbYz8IyznlTB8RMlPDl7K33jm3OVDvyK1FhyUgxX9o9j6sKdpO0/5nYcr1VpAViPk9dfC3a+LHAe8IGzfBow1rk9xrmP8/goY4yOYlbBKwt3se9wIf932Rk68CtSS5Mv6U5IUAB/+kRTRp9OlY4BGGMCjTFrgBxgDrADyLPWntzBlgHEObfjgHQA5/F8oEUFP3OiMSbFGJOSm5tbu9+iEdiXX8g/F+xgdO+2DErSVb5Eaiu2WRj3jurMvC05fLNd7zEVqVIBWGtLrbX9gHhgMNCjotWc7xV9dP3BeCxr7VRrbbK1NrlVq1ZVzdto/e3LrZSWWR7SVb5E6sz4YUkkxITz5083U6p5gn6gWqOArLV5wNfAECDKGHPyklTxQJZzOwNIAHAebw4crIuwjdWGzHz++20Gtw5PIrFFhNtxRBqN0KBAHrq4O1v2HeG/32a4HcfrVGUUUCtjTJRzOxw4H9gMzAeudlYbD8x0bs9y7uM8Ps/qjIzTstby2KebiI4I4efndnY7jkijM7p3W/olRPG32Vs1LPQUVdkCaAvMN8asA1YCc6y1nwAPAQ8YY1Lx7ON/zVn/NaCFs/wBYHLdx2485mzKZtnOg9x/fheah2vYp0hdM8bwu8t6kHOkiFcW7nI7jlep9Kri1tp1QP8Klu/Eczzg1OWFwDV1kq6RO1FSxpTPt9A5tgk3DE50O45IozWwfQyX9m7Dywt3cMPgBGKbhbkdySvoTGAXTV+2m137j/Hb0T0ICtQ/hUh9+vVF3SkuLePpObp85El613FJ3vETPDd3O2d1ack5XTUKSqS+JbWM5OYhSbyXks6WfYfdjuMVVAAueXbudo4UFvPb0T3QeXIiDePeUZ1pEhrElM+2uB3FK1R6DEDqXtr+Y0xfupvrBiXSvU0zt+OI+I2oiBDuHdWFxz7dzMBH53Dw2AnaRYUz6aJujO0fV/kPaGRUAC54as42ggMDuP+CLm5HEfE7UeHBGODAsRMAZOYV8PCM9QB+VwLaBdTANmTm8/HaLCaM6EBsU41EEGloz3y1/QdTExQUl/Lk7K2u5HGTCqCBPTl7K1ERwUw8u6PbUUT8UlZeQbWWN2YqgAa0bOcBFmzL5Wdnd6JZmE76EnFDu6jwai1vzFQADcRayxNfbKF1s1DGD0tyO46I35p0UTfCgwO/tywsKIBJF/nfFfhUAA1k7uYcvt2Txy/P70rYKS8+EWk4Y/vHMeXK3sSV+8R/TrdYvzsADBoF1CBKyyxPzt5Kh5aRXDNQV/oScdvY/nHfveHfNX0Vi1L3c/DYCWIiQ1xO1rC0BdAAZq7JZGv2ER68sKumfBDxMg9e2JVjJ0r45wL/u36w3o3q2YkSz9wjveKacWmvtm7HEZFTdGndlJ/0j2PakjT25Re6HadBqQDq2Tsr9pBxqIBJF3XXdX5FvNT953elzFqen7fd7SgNSgVQj44VlfD8vO0M6RjDyC4t3Y4jIqeREBPB9YMSeXdlOnsOHHc7ToNRAdSjNxbvYv/RE/z64u6a8E3Ey91zXmeCAg1//8p/potWAdST/OPFvLxgJxec0ZoBidFuxxGRSsQ2C2P8sCQ+XJPJtuwjbsdpECqAevLaop0cKSrhgQu6uh1FRKrorpGdaBISxFNf+se8QCqAepB3/ARvLE7j0t5t6NFW0z2L+IroyBDuGNmR2RuzWZue53aceqcCqAevL9rFkaIS7h2l6Z5FfM1tIzoQExnCU35w6UgVQB0r/+lfF3sR8T1NQoOYOLIjC7flsnrPIbfj1CsVQB3Tp38R33fzkPZERwTz3NzGfV6ACqAO6dO/SOMQGRrE7Wd1ZP7W3EZ9LEAFUIf06V+k8Rg/LImoiOBGfXawCqCO6NO/SOPSJDSICcM78NXmHDZk5rsdp16oAOqIPv2LND7jhyfRLCyo0R4LUAHUAX36F2mcmoUFc9uIDny5KZtNWYfdjlPnVAC18NHqTIb/dR79/jSHI0Ul9GrX3O1IIlLHbh3egaahjXMrQAVQQx+tzuThGevJzCv4btnz81L5aHWmi6lEpK41Dw/m1uFJfLFxH1v2Na6tABVADT05eysFxaXfW1ZQXMqTs/1jDhERf3LbiA40CQ3i+bmpbkepUyqAGsoq98m/KstFxHdFRYRwy7AkPtuwt1HNFKoCqKF2UeHVWi4ivm3CiA5EBAfy/LzGsxWgAqih+yoY7hkeHMiki7q5kEZE6lt0ZAjjhiXxybosUnMax1aACqCGCks8+/9bNgnBAHFR4Uy5sjdj+8e5G0xE6s3tIzoQFhTIS1/vdDtKnQiqbAVjTALwFtAGKAOmWmufNcbEAO8CSUAacK219pDxXPvwWeBS4Dhwi7X22/qJ747i0jJeXrCTge2j+eCuobrco4ifaNEklBsGJ/LW0jQeuLArcT6+y7cqWwAlwIPW2h7AEOAXxpgzgMnAXGttF2Cucx/gEqCL8zUReKnOU7ts1posMvMK+MW5nfTmL+Jnbj+rAwCvfuP7WwGVFoC1du/JT/DW2iPAZiAOGANMc1abBox1bo8B3rIey4AoY0zbOk/ukrIyy0sLdtC9TVPO7RbrdhwRaWDtosIZ0y+O/6xI5+CxE27HqZVqHQMwxiQB/YHlQGtr7V7wlARw8t0wDkgv97QMZ9mpP2uiMSbFGJOSm5tb/eQumbM5m9Sco/zsHH36F/FXd53dkYLiUqYtSXM7Sq1UuQCMMU2A/wK/tNb+2OlwFb0r2h8ssHaqtTbZWpvcqlWrqsZwlbWWF+en0r5FBKN7N5qNGhGppi6tm3LBGa2ZtjSNY0UlbsepsSoVgDEmGM+b/7+ttTOcxdknd+0433Oc5RlAQrmnxwNZdRPXXUt2HGBtRj53juxEUKAGUIn4s5+d04m848X8Z2V65St7qUrfxZxRPa8Bm621T5d7aBYw3rk9HphZbvk44zEEyD+5q8jXvfh1KrFNQ7lqoIZ6ivi7AYnRnNkhhle/2cmJkjK349RIVT7GDgduBs4zxqxxvi4F/gpcYIzZDlzg3Af4DNgJpAKvAD+v+9gNb016HotTD3DHWR0JDQp0O46IeIGfndOJvfmFzFzjm5NAVnoegLV2ERXv1wcYVcH6FvhFLXN5nRfnp9I8PJgbzkx0O4qIeImzu7aiR9tm/HPBDq4aEE9AgG8NDNGO7CrYnn2ELzfkr5gYAAAKTElEQVRlM35YEk1CK+1MEfETxhh+dk4nduQeY87mbLfjVJsKoApeWrCDiJBAbh2W5HYUEfEyl/ZqQ2JMBC9+vQPPDhDfoQKoxL78QmatyeLa5ASiI0PcjiMiXiYoMICJIzuyNj2PZTsPuh2nWlQAlXhzSRpl1jJhRAe3o4iIl7p6YDwtm4Ty0oIdbkepFhXAjzhWVMLby3dzSa+2JMREuB1HRLxUWHAgt41IYuG2XDZk5rsdp8pUAD/ivZR0DheWfDf5k4jI6dx4ZnsiQwJ5bdEut6NUmQrgNEpKy3h98S6S20fTPzHa7Tgi4uWahwdz7aAEPl6bxb78QrfjVIkK4DRmb8wm/WABt5/V0e0oIuIjbh3WgTJrmbY0ze0oVaICqIC1lle+2Un7FhFccEZrt+OIiI9IbBHBRT3b8O9lu31ikjgVQAVW7T7EmvQ8JozoQKCPndknIu66/awOHC4s4b/fZrgdpVIqgAq88s1OmocHc/XAeLejiIiPGZAYTb+EKF5ftIvSMu8+MUwFcIq0/cf4clM2Nw9pT0SIpn0QkeoxxnD7WR1IO3CcuV4+PYQK4BSvLdpFcEAA44a1dzuKiPioi3u2IS4qnFe9fEioCqCcQ8dO8P6qdMb0a0ds0zC344iIjwoKDODW4Ums2HWQdRl5bsc5LRVAOf9evpvC4jIN/RSRWrt2UAJNQoO8+sQwFYCjqKSUaUt3M7JrK7q1aep2HBHxcc3CgrluUAKfrttLVl6B23EqpAJwzFyTRe6RIu7QtA8iUkduGZbk1SeGqQDwnPj1xuI0urVuyojOLd2OIyKNREJMBJf0asvby/d45YlhKgBgZdohNu89zC3DkzBGJ36JSN2ZcFYHjhSW8H5KuttRfkAFALy5ZBfNw4MZ2y/O7Sgi0sgMSIxmQGIUry9O87oTw/y+ALLyCpi9MZvrBycQHhLodhwRaYQmjOjInoPHmbclx+0o3+P3BfCvZbux1nLzEJ34JSL148KerWnTLIy3lqa5HeV7/LoACotLeWfFHi44ozXx0bril4jUj+DAAG4aksg32/eTmnPU7Tjf8esCmLUmi0PHi7llmIZ+ikj9un5wIiGBAUxfmuZ2lO/4bQFYa3lziWfo55COMW7HEZFGrmWTUC7r05YPVmVwpLDY7TiAHxfAyrRDbNLQTxFpQOOGJXHsRCn/XeUd1wrw2wLQ0E8RaWj9EqLomxDFW0t3U+YFQ0L9sgC+G/o5SEM/RaRhjR/anp37j7Eodb/bUfyzAE4O/bxJQz9FpIGN7tOWlk1CvGJIqN8VQPmhnwkxGvopIg0rNCiQGwYnMndLDnsOHHc1i98VwKy1nqGf44cluR1FRPzUT89MJMAY/rV8t6s5/KoArLW86cz6ObRjC7fjiIifats8nIt7tuHdlekUnCh1LYdfFUDKbg39FBHvMG5oe/ILipm5JtO1DH5VANOX7qZpWBBj+rVzO4qI+LnBHWLo3qYpby5Jw1p3hoT6TQHsP1rE5xv2ctWAeCJCgtyOIyJ+zhjDLcOS2LLvCCt2HXQlQ6UFYIx53RiTY4zZUG5ZjDFmjjFmu/M92llujDHPGWNSjTHrjDED6jN8dbyfkkFxqeWmIYluRxERAWBMvziahwfz1lJ3DgZXZQvgTeDiU5ZNBuZaa7sAc537AJcAXZyvicBLdROzdsrKLG+v2M2QjjF0jtUF30XEO4SHBHLdoAS+2LiP7MOFDf7fr7QArLULgVO3T8YA05zb04Cx5Za/ZT2WAVHGmLZ1FbamFm7PJf1ggU78EhGv89PBiZSWWd5b2fCXjKzpzvDW1tq9ANbavcaYWGd5HFD+t8hwlu099QcYYybi2UoAKCq/i6m+XP54vf3oloD753XXnPK7x5ezg/LXmXsfh3ur/7Rutflv1vXR0IrGVlZ4eNtaOxWYCmCMSbHWJtdxlgaj/O7y5fy+nB2U323GmJTaPL+mo4CyT+7acb6fvNBlBpBQbr14IKvm8UREpL7UtABmAeOd2+OBmeWWj3NGAw0B8k/uKhIREe9S6S4gY8w7wDlAS2NMBvAH4K/Ae8aYCcAe4Bpn9c+AS4FU4DhwaxVzTK1ebK+j/O7y5fy+nB2U3221ym/cOgNNRETc5TdnAouIyPepAERE/FSDF4AxJswYs8IYs9YYs9EY84izvIMxZrkzvcS7xpiQhs5WVcaYQGPMamPMJ859X8qeZoxZb4xZc3II2emm9vBGxpgoY8wHxpgtxpjNxpihvpLfGNPN+buf/DpsjPmlr+QHMMbc7/x/u8EY847z/7NPvP6NMfc5uTcaY37pLPPqv319T8XjxhZAEXCetbYv0A+42Bkx9DjwjDO9xCFgggvZquo+YHO5+76UHeBca22/cuOfTze1hzd6FvjCWtsd6Ivn38En8ltrtzp/937AQDwDJT7ER/IbY+LwnKuUbK3tBQQC1+MDr39jTC/gDmAwntfNZcaYLnj/3/5N6nMqHmuta19ABPAtcCaes/GCnOVDgdluZvuRzPHOH/084BM8J7/5RHYnXxrQ8pRlW4G2zu22wFa3c54mezNgF87gBV/Lf0rmC4HFvpSf/53pH4NnBOEnwEW+8PrHM1Lx1XL3fwf82hf+9kASsKHc/QozAy8DN1S03um+XDkG4OxCWYPnBLI5wA4gz1pb4qxycgoJb/R3PC+cMud+C3wnO3jOzP7SGLPKmY4DTpnaA4g97bPd1RHIBd5wdsG9aoyJxHfyl3c98I5z2yfyW2szgb/hGfq9F8gHVuEbr/8NwEhjTAtjTASe4eoJ+Mjf/hSny3y6qXhOy5UCsNaWWs9mcDyeTbIeFa3WsKkqZ4y5DMix1q4qv7iCVb0ueznDrbUD8Gwu/sIYM9LtQNUQBAwAXrLW9geO4X2b7JVy9pFfAbzvdpbqcPY1jwE6AO2ASDyvo1N53evfWrsZz66qOcAXwFqg5Eef5Huq/V7k6igga20e8DUwBM/MoSdPTPPWKSSGA1cYY9KA/+DZDfR3fCM7ANbaLOd7Dp79z4M5/dQe3iYDyLDWLnfuf4CnEHwl/0mXAN9aa7Od+76S/3xgl7U211pbDMwAhuEjr39r7WvW2gHW2pF4Zjjeju/87curs6l43BgF1MoYE+XcDsfzotoMzAeudlYrP72E17DWPmytjbfWJuHZhJ9nrb0RH8gOYIyJNMY0PXkbz37oDZx+ag+vYq3dB6QbY07OgDgK2ISP5C/nBv63+wd8J/8eYIgxJsIYY/jf399XXv+xzvdE4Eo8/wa+8rcvr+6m4nHhgEYfYDWwDs+bz++d5R2BFXimkXgfCHX74Eslv8c5wCe+lN3Judb52gj81lneAs+B7e3O9xi3s/7I79APSHFePx8B0T6WPwI4ADQvt8yX8j8CbHH+350OhPrQ6/8bPIW1FhjlC397PCW1FyjG8wl/wuky49kF9AKeY6rr8YzW+tGfr6kgRET8lM4EFhHxUyoAERE/pQIQEfFTKgARET+lAhAR8VMqABERP6UCEBHxU/8Pso6WSZ3porEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.xlim(30, 100)\n",
    "plt.ylim(100, 600)\n",
    "\n",
    "X = np.array([40, 50, 60, 70, 80, 90])\n",
    "Y = np.array([216, 399, 496, 507, 432, 271])\n",
    "plt.scatter(X, Y)\n",
    "\n",
    "x = np.linspace(30, 100)\n",
    "y = 57*x - 0.43*x*x - 1376\n",
    "plt.plot(x, y)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一元二次回归的求解可以转换为二元线性回归的求解。依葫芦画瓢，当我们遇到一元三次，一元四次，或 一元 N 次曲线时，也可以将自变量 $x$ 扩充为 $x^3$，$x^4$ 或 $x^N$，从而转换为 N 元线性回归的求解。如果自变量不止一个，譬如二元三次曲线，也可以使用同样的方法，只不过要注意自变量之间的组合。这被称为线性回归的扩展。\n",
    "\n",
    "我们用数学语言来描述的更清晰一点。对于一元线性模型，我们记为：\n",
    "\n",
    "$$\n",
    "y = ax + b\n",
    "$$\n",
    "\n",
    "其中，a、b、x 都是标量，如果把 x 变成关于 x 的函数向量 $\\phi(x)$，同时把参数 a 变成关于函数向量 $\\phi(x)$ 的参数向量 $\\bf{\\theta}$，那么可以得到线性模型的一个扩展：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y &= \\theta_1\\phi_1(x) + \\theta_2\\phi_2(x) + ... + \\theta_k\\phi_k(x) \\\\\n",
    "&= \\sum_{i=1}^k\\theta_i\\phi_i(x) \\\\\n",
    "&= \\bf{\\theta}^T\\phi(x)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "这里我们把 $\\phi(x)$ 叫做 **基函数**（basis function），当我们把 $\\phi(x)$ 定义为：\n",
    "\n",
    "$$\n",
    "\\phi(x) = (1, x, x^2)^T\n",
    "$$\n",
    "\n",
    "我们就得到了上面的一元二次回归模型，更一般的，当我们把 $\\phi(x)$ 定义为多项式形式：\n",
    "\n",
    "$$\n",
    "\\phi(x) = (1, x, x_2, ..., x^{k-1})^T\n",
    "$$\n",
    "\n",
    "我们就可以得到更通用的 **多项式回归** 模型，当我们把 $\\phi(x)$ 定义为三角多项式形式：\n",
    "\n",
    "$$\n",
    "\\phi(x) = (1, \\sin x, \\cos x, \\sin 2x, \\cos 2x, ..., \\sin mx, \\cos mx)^T\n",
    "$$\n",
    "\n",
    "我们就可以得到了 **三角多项式回归** 模型。像这样把本来是一维的模型扩展成多维模型，叫做 **升维**。\n",
    "\n",
    "> **基函数**是构成**函数空间**的基础，就好像**基向量**是构成**向量空间**（线性空间）的基础一样，函数空间中的任何一个函数都可以表示成基函数的线性组合。在数值分析和近似理论中，基函数也被称为**混合函数**（blending function），常常用于**插值**（interpolation）。上面的多项式回归使用的基函数叫做**多项式基**，三角多项式回归使用的基函数叫做**傅里叶基**。此外，基函数还是**核方法**（Kernel method）和**字典训练**的基础。\n",
    "\n",
    "但是这种扩展只能针对输入变量是一维的情况，对于多维的输入向量 $\\bf{x}$，基函数该如何选取呢？对于这种情况，一般有两种扩展方法：一种是通过基函数相乘来得到扩展模型，这种模型称为 **乘法模型**，另一种是通过基函数相加得到扩展模型，这种模型称为 **加法模型**。假设 $\\bf{x}$ 是一个 d 维向量 $(x_1, ... x_d)^T$，n 表示各维的参数个数，那么乘法模型可以记为：\n",
    "\n",
    "$$\n",
    "y = \\sum_{j_1=1}^{n} \\dots \\sum_{j_d=1}^{n} \\theta_{j_1 \\dots j_d} \\phi_{j_1}(x_1) \\dots \\phi_{j_d}(x_d)\n",
    "$$\n",
    "\n",
    "加法模型可以记为：\n",
    "\n",
    "$$\n",
    "y = \\sum_{i=1}^d \\sum_{j=1}^n \\theta_{i,j} \\phi_{j}(x_i)\n",
    "$$\n",
    "\n",
    "上面的公式看起来可能很复杂，不过我们可以通过一个简单的例子把上式展开来看看，假设维数 d = 2，也就是 $\\bf{x} = (x_1, x_2)^T$，基函数为多项式形式并且 n = 3，也就是 $\\phi(x) = (1, x, x^2, x^3)$，很容易得出乘法模型：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y &= \\sum_{i=1}^3 \\sum_{j=1}^3 \\theta_{ij} \\phi_i(x_1) \\phi_j(x_2) \\\\\n",
    "&= \\sum_{i=1}^3 (\\theta_{i1} \\phi_i(x_1) \\phi_1(x_2) + \\theta_{i2} \\phi_i(x_1) \\phi_2(x_2) + \\theta_{i3} \\phi_i(x_1) \\phi_3(x_2)) \\\\\n",
    "&= \\theta_{11} \\phi_1(x_1) \\phi_1(x_2) + \\theta_{12} \\phi_1(x_1) \\phi_2(x_2) + \\theta_{13} \\phi_1(x_1) \\phi_3(x_2) \\\\\n",
    "&+ \\theta_{21} \\phi_2(x_1) \\phi_1(x_2) + \\theta_{22} \\phi_2(x_1) \\phi_2(x_2) + \\theta_{23} \\phi_2(x_1) \\phi_3(x_2) \\\\\n",
    "&+ \\theta_{31} \\phi_3(x_1) \\phi_1(x_2) + \\theta_{32} \\phi_3(x_1) \\phi_2(x_2) + \\theta_{33} \\phi_3(x_1) \\phi_3(x_2) \\\\\n",
    "&= \\theta_{11} x_1 x_2 + \\theta_{12} x_1 x_2^2 + \\theta_{13} x_1 x_2^3 \\\\\n",
    "&+ \\theta_{21} x_1^2 x_2 + \\theta_{22} x_1^2 x_2^2 + \\theta_{23} x_1^2 x_2^3 \\\\\n",
    "&+ \\theta_{31} x_1^3 x_2 + \\theta_{32} x_1^3 x_2^2 + \\theta_{33} x_1^3 x_2^3\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "加法模型：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y &= \\sum_{i=1}^2 \\sum_{j=1}^3 \\theta_{ij} \\phi_j(x_i) \\\\\n",
    "&= \\sum_{i=1}^2 (\\theta_{i1} \\phi_1(x_i) + \\theta_{i2} \\phi_2(x_i) + \\theta_{i3} \\phi_3(x_i)) \\\\\n",
    "&= \\theta_{11} \\phi_1(x_1) + \\theta_{12} \\phi_2(x_1) + \\theta_{13} \\phi_3(x_1) \\\\\n",
    "&+ \\theta_{21} \\phi_1(x_2) + \\theta_{22} \\phi_2(x_2) + \\theta_{33} \\phi_3(x_3) \\\\\n",
    "&= \\theta_{11} x_1 + \\theta_{12} x_1^2 + \\theta_{13} x_1^3 \\\\\n",
    "&+ \\theta_{21} x_2 + \\theta_{22} x_2^2 + \\theta_{33} x_3^3\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "乘法模型是将基函数两两相乘，要回归的参数个数为 $b^d$，而加法模型只是将所有的基函数相加，所以要回归的参数个数为 $b \\times d$，很显然，乘法模型的表现力要丰富的多，但是乘法模型的参数个数是呈指数级増长的，当输入参数的维度增加到 100 时，要回归的参数个数就增加到 $2^100$ 个，这个计算量是非常大的。这种随着维度的增加，计算量呈指数级増长的现象，通常称为 **维数灾难** 或 **维数诅咒**（curse of dimensionality）。\n",
    "\n",
    "无论是上面介绍的一维情况下的多项式回归模型，还是多维情况下的乘法模型或加法模型，本质上都还是线性回归。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们假设基函数 $\\phi(x) = (\\phi_1(x), \\phi_2(x), \\dots, \\phi_b(x))$，于是有下面的 $n \\times b$ 阶矩阵，被称为**设计矩阵**：\n",
    "\n",
    "$$\n",
    "\\Phi = \n",
    "\\left ( \n",
    "\\begin{array}{ccc}\n",
    "\\phi_1(x_1) & \\phi_2(x_1) & \\dots & \\phi_b(x_1) \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "\\phi_1(x_n) & \\phi_2(x_n) & \\dots & \\phi_b(x_n)\n",
    "\\end{array} \n",
    "\\right )\n",
    "$$\n",
    "\n",
    "基于基函数的线性模型可以表示成：\n",
    "\n",
    "$$\n",
    "f_{\\theta}(x) = \\Phi \\theta\n",
    "$$\n",
    "\n",
    "损失函数为：\n",
    "\n",
    "$$\n",
    "J(\\theta) = \\frac{1}{2} \\| \\Phi \\theta - y \\|^2\n",
    "$$\n",
    "\n",
    "和前面多元线性回归的求解一样，我们得到其最小二乘解为：\n",
    "\n",
    "$$\n",
    "\\hat{\\theta} = (\\Phi^T\\Phi)^{-1}\\Phi^T y = \\Phi^{\\dagger}y\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用这种方法我们可以对线性模型进行无限的扩展，譬如我们把设计矩阵 $\\Phi$ 替换为下面的 **核矩阵** $K$：\n",
    "\n",
    "$$\n",
    "K = \n",
    "\\left ( \n",
    "\\begin{array}{ccc}\n",
    "K(x_1, x_1) & K(x_1, x_2) & \\dots & K(x_1, x_n) \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "K(x_n, x_1) & K(x_n, x_2) & \\dots & K(x_n, x_n)\n",
    "\\end{array} \n",
    "\\right )\n",
    "$$\n",
    "\n",
    "我们就得到了 **核模型**：\n",
    "\n",
    "$$\n",
    "f_{\\theta}(x) = \\sum_{i=1}^n \\theta_i K(x, x_i)\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}