{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 最小二乘\n", "\n", "## 1. 最小二乘的基本原理\n", "\n", "最小二乘法(Least Squares)是一种数学优化技术,它通过最小化误差的平方和找到一组数据的最佳函数匹配, 最小二乘法通常用于曲线拟合、求解模型。很多其他的优化问题也可通过最小化能量或最大化熵用最小二乘形式表达。\n", "\n", "![ls_theory](images/least_squares.png)\n", "\n", "最小二乘原理的一般形式为:\n", "$$\n", "L = \\sum (V_{obv} - V_{target}(\\theta))^2\n", "$$\n", "其中\n", "* $V_{obv}$是观测的多组样本值\n", "* $V_{target}$是假设拟合函数的输出值\n", "* $\\theta$为构造模型的参数\n", "* $L$是目标函数\n", "\n", "如果通过调整模型参数$\\theta$,使得$L$下降到最小则表明,拟合函数与观测最为接近,也就是找到了最优的模型。\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 示例\n", "\n", "假设我们有下面的一些观测数据,希望找到它们内在的规律。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU8ElEQVR4nO3df6yeZX3H8ffXw1EPuu2ANIS21DYZaUdwUj0hOLJNClicRBrnnEZNsxD5RzdGTLXsR8wSF+u6oSZzLI2i3SSCg64QZVYCdWZmOk8tDvnRwBBKD2jrpM5pM0r57o/zHCinzzk9z3Oe+9dzv19J0/Pc58d9nbvp57nu6/pe1x2ZiSSpPV5SdQMkSeUy+CWpZQx+SWoZg1+SWsbgl6SWMfglqWVOKfoEETECTAJTmXlFRKwCbgZeBewB3puZz8z3M84444xcuXJl0U2VpKGyZ8+eH2fmktnHCw9+4BrgQeCXO68/DnwiM2+OiL8HrgJumO8HrFy5ksnJyWJbKUlDJiIe73a80KGeiFgOvAX4TOd1AOuAWztfsh3YUGQbJEkvVvQY/yeBDwHPdV6/Cjicmc92Xh8AlhXcBknScQoL/oi4AjiYmXv6/P6rI2IyIiYPHTo04NZJUnsV2eO/CHhrRDzG9GTuOuBTwHhEzMwtLAemun1zZm7LzInMnFiy5IS5CUlSnwoL/sy8LjOXZ+ZK4J3APZn5bmA38PbOl20Ebi+qDZKkE5VR1TPbh4GbI+KjwF7gsxW0QZJqa+feKbbu2seTh4+wdHyMTetXs2Ht4KZDSwn+zPw68PXOx48CF5RxXklqmp17p7hux30cOXoMgKnDR7hux30AAwt/V+5KUo1s3bXv+dCfceToMbbu2jewcxj8klQjTx4+0tPxfhj8klQjS8fHejreD4Nfkmpk0/rVjI2OvOjY2OgIm9avHtg5qqjqkSTNYWYCt/FVPZKkhduwdtlAg342h3okqWUMfklqGYNfklrG4JeklnFyV9JQK3rfm7qcsxcGv6ShVca+N3U4Z68c6pE0tMrY96YO5+yVwS9paJWx700dztkrh3okDa2l42NMdQncbvveDGpcvpdzVsUev6ShtdB9b2bG5acOHyF5YVx+596uT4YdyDmrZPBLGlob1i7jY297DcvGxwhg2fgYH3vba07oyQ9yXH6h56ySQz2ShtpC9r0Z9Lh80XvtLJY9fkmtV8Ye+HVi8EtqvYvXLCFmHavbuPwgGfySWm3n3ilu2zNFHncsgN99fb2HaxbD4JfUat0mdhPY/dChahpUAoNfUqs1YcHVoFnVI6nV5ltwVffN1vplj19Sq8214OriNUsGtqirbgx+Sa0214Kr3Q8dqv1ma/1yqEdS63VbcHXtLfd2/dphGPu3xy9JXQzzoi6DX5K6aMJma/1yqEeSupgZ+hnGqh6DX9LQGVQZZt03W+uXwS9pqDThmbfdlLlmwDF+SUOlCc+8nW2QD4JZCINf0lBp4hYMZb9ZGfyShkoTyzDLfrMy+CUNlSaWYZb9ZmXwSxoqTXjm7Wxlv1lZ1SNp6DStDLPsNQMGvyTVQJlvVoUN9UTEyyPiPyLiexFxf0T8Ref4qoj4dkQ8EhG3RMRLi2qDJOlERY7x/x+wLjNfC5wPXB4RFwIfBz6Rmb8KPA1cVWAbJEmzFBb8Oe1/Oy9HO38SWAfc2jm+HdhQVBskSScqtKonIkYi4l7gIHAX8F/A4cx8tvMlB4Cug1oRcXVETEbE5KFDw/vQY0kqW6GTu5l5DDg/IsaBfwbW9PC924BtABMTE1lIAyU1xrA+/7YKpVT1ZObhiNgNvAEYj4hTOr3+5UDzH2ApqVBN3XitrgoL/ohYAhzthP4YcBnTE7u7gbcDNwMbgduLaoPUJPZoT7wGF69Zwu6HDjHVZeuCmb1s2naNBqHIHv9ZwPaIGGF6LuFLmfnliHgAuDkiPgrsBT5bYBukRrBH2/0afOFb++f9njpvvFZnhQV/Zv4nsLbL8UeBC4o6r9RE8+3O2Jbg73YNTma+vWy8g5qbK3elGmjiVsKD1uvvOt9eNt5Bzc9N2qQaaOJWwoPWy+96so3XmvgwljIZ/FINNHEr4UHrdg1mG31JcNqpozx5+Ahbd+2b8wlV3kHNz+CXaqCJWwkPWrdr8J4LVzz/enxsFAKe/sXRkz6e0Duo+UVm/ddGTUxM5OTkZNXNkFSBmUnabiWdMP0G8c3N6074nuPH+GH6Dqptb6YRsSczJ2Yfd3JXUm11C/DZug3flL2/fdMY/JJqayElnkvHx+Ys3TTouzP4JdXWySZjx0ZHuHjNEks3e+TkrqTamm8ydmYCfPdDhyzd7JHBL6m25ipz/eTvn883N69jw9pllm72weCXVFsLKXO1dLN3jvFLqrWTTdJuWr+6a+lmmxa/9crgl9Rolm72zuCXdFJ13+nS0s3eGPyS5uVOl8PH4Jc0rzo+K6DudyB1Z/BLBRqGgKpbuWRd7kCa/G9rOadUkD/beR/X3nIvU4ePnHQ3yTqrW7lkHfban3nzaeq/rcEvFWDn3ilu+tZ+Zu9928QVpXV7VkAd7kDq8OazGA71SB2DvHXfumvfCaE/o2krSutWLrl0fKzrFs1l3oHU4c1nMQx+icGPG88XAE1cUVqncsk6LNiqw5vPYjjUIzH4W/e5AiDAFaWLVIenldVt+KtX9vglBn/r3q1XGsC7L1xRm55zk1V9B1K34a9eGfwSg791b3ow6OSqfvNZDINfYrqHvunW73H02AtTsqMjsahb97mCoddJ5CbXi6ueDH5pxuwynLnKchah10nkuixW0nBxcldiekjm6HMvTvqjz+XA67J7nURuer246snglyivLrvX8zS9Xlz15FCPRHl12b2ep2714s43DAd7/BLl1WX3ep461Ys3fX8avcDglyhvUVCv56nDYqUZzjcMj8gsoHRhwCYmJnJycrLqZkittmrzV7oWOgXwgy1vKbs5WoCI2JOZE7OP2+OXtCB1255Z/TP4JS1IneYbtDhW9UgnYSXLtEFvQ+F1rY7BL82jypWzdQzGQe1P44rkajnUI82jqkqWYS+dtEKoWga/NI+qVs4OezC6IrlahQV/RJwdEbsj4oGIuD8irukcPz0i7oqIhzt/n1ZUG6TFqqqSZdiD0QqhahXZ438W+GBmngtcCLw/Is4FNgN3Z+Y5wN2d11ItVVXJMuzBaIVQtQoL/sx8KjO/2/n4Z8CDwDLgSmB758u2AxuKaoO0WFWtnB32YKzTiuQ2KmXlbkSsBL4BnAfsz8zxzvEAnp55Pet7rgauBlixYsXrH3/88cLbKc2oQ0XNYtpQh/arenOt3C08+CPilcC/An+ZmTsi4vDxQR8RT2fmvOP8btmgMs0uNYTp3nZTeqRNb78Gp5ItGyJiFLgNuCkzd3QO/ygizup8/izgYJFtkHrV9IqaprdfxSuyqieAzwIPZub1x33qDmBj5+ONwO1FtUHqR9MraprefhWvyB7/RcB7gXURcW/nz+8AW4DLIuJh4NLOa6k2ml5R0/T2q3hFVvX8W2ZGZv56Zp7f+XNnZv53Zl6Smedk5qWZ+ZOi2iD1o+kVNU1vv4rnXj3SLIPejKxsTW+/iueDWCRpSPkgFkkSME/wR8SdnYVXkqQhMt8Y/+eAr0XEduCvMvNoSW3SEHIlaf+8dhq0OYM/M/8pIv4F+HNgMiL+EXjuuM9fP9f3SsfzoRv989qpCCcb438G+DnwMuCXZv2RFsSVpP3z2qkIc/b4I+Jy4HqmV9q+LjN/UVqrNFRcSdo/r52KMF+P/0+B38vMzYa+FsOVpP3z2qkIcwZ/Zv5mZt5fZmM0nFxJ2j+vnYrgyl0VzpWk/Zvv2lnto365cldqIPfc10K4clcaIlb7aDEMfqmBrPbRYhj8UgNZ7aPFMPilBrLaR4thVY/UQFZKaTEMfqmhNqxdZtCrLwZ/zVmrLWnQDP4ac2dGSUVwcrfGrNWWVAR7/DVmrXb7OLSnMtjjrzFrtdtlZmhv6vARkheG9nbunaq6aRoyBn+NzVWrffGaJVy05R5Wbf4KF225x2AYEg7tqSwO9dRYt1rti9cs4bY9U074DiGH9lQWg7/mZtdqX7Tlnjl7hQZ/sy0dH2OqS8g7tKdBc6inBDv3Tg1saMZe4fByGwaVxR5/wQZdi9+0XqFVKgvnNgwqi8FfsPkm7Pr5D71p/equD+CoY6/QBWi9cxsGlcGhnoINemhmw9plfOxtr2HZ+BgBLBsfq+1Tl6xSkerJHn/BihiaaUqv0PkIqZ7s8ReszRN2C12ANsjJb0knZ4+/YG2bsDt+Mnf81FFGXxIcfS6f//zsNz3nAaTyGfwlaMrQzGLNDvGnf3GU0ZFgfGyUnx452vVNb9CT35JObuiD33LC8nQL8aPHkle87BTu/cibun6P8wBS+YZ6jN9Nr8rVT4i7EZ1UvqEN/p17p/jgl75nOWGJ+gnxNk9+S1UZyuCf6ekfy+z6eYcRitFPiDdpXYI0LAob44+IG4ErgIOZeV7n2OnALcBK4DHgHZn59KDP3W2s+XgOIxSj3wqmtkx+S3VR5OTu54G/Bf7huGObgbszc0tEbO68/vCgTzxfj95hhGIZ4lL9FRb8mfmNiFg56/CVwBs7H28Hvk4BwT/XatmRiFYPI1jhJAnKH+M/MzOf6nz8Q+DMIk4y11jz37zjta0NOiucJM2orI4/MzMius++AhFxNXA1wIoVK3r62W1bLbsQRSyU8g5Caqayg/9HEXFWZj4VEWcBB+f6wszcBmwDmJiYmPMNYi6ONb/YoBdKzbXVwuTjP2H3Q4d8M5BqrOyhnjuAjZ2PNwK3l3z+1hr0Qqm57iBu+tZ+h5Okmiss+CPii8C/A6sj4kBEXAVsAS6LiIeBSzuvVYJBL5Sa605h9q2ZC+ak+imyquddc3zqkqLOqbkNet5jrsqpblwwJ9XL0G/SphcMct6j2yMggxN7/OCCOaluhnLLBhWv21YL775whfvuSA1gj79Ew1b+2O0OYuLVpw/V7ygNI4N/ABYS6G150pRltFL9OdSzSAtdETvfAipJKpPBv0gLDXSfNCWpLgz+RVpooPukKUl1YfAv0kID3SdNSaoLg3+RFhroPmlKUl1Y1bNIvayIteJFUh0Y/ANgoEtqEod6JKllDH5JahmDX5JaxuCXpJYx+CWpZQx+SWoZg1+SWsbgl6SWcQFXCYbtASySms3gL1hbHsAiqTkc6imYD2CRVDcGf8F8AIukujH4C+YDWCTVjcFfMB/AIqlunNwtWC/79UtSGQz+Erhfv6Q6cahHklrG4JekljH4JallDH5JahmDX5JaxuCXpJYx+CWpZQx+SWoZg1+SWsbgl6SWMfglqWUMfklqmUqCPyIuj4h9EfFIRGyuog2S1FalB39EjACfBt4MnAu8KyLOLbsdktRWVfT4LwAeycxHM/MZ4GbgygraIUmtVEXwLwOeOO71gc6xF4mIqyNiMiImDx06VFrjJGnY1XZyNzO3ZeZEZk4sWbKk6uZI0tCoIvingLOPe728c0ySVIIqgv87wDkRsSoiXgq8E7ijgnZIUiuV/szdzHw2Ij4A7AJGgBsz8/6y2yFJbVXJw9Yz807gzirOLUltV9vJXUlSMQx+SWoZg1+SWsbgl6SWMfglqWUMfklqmUrKOeto594ptu7ax5OHj7B0fIxN61ezYe0JWwhJUuMZ/EyH/nU77uPI0WMATB0+wnU77gMw/CUNHYd6gK279j0f+jOOHD3G1l37KmqRJBXH4AeePHykp+OS1GQGP7B0fKyn45LUZAY/sGn9asZGR150bGx0hE3rV1fUIkkqjpO7vDCBa1WPpDYw+Ds2rF1m0EtqBYd6JKllDH5JahmDX5JaxuCXpJYx+CWpZSIzq27DSUXEIeDnwI+rbkvFzqDd16Dtvz94Ddr++0Nv1+DVmblk9sFGBD9ARExm5kTV7ahS269B239/8Bq0/feHwVwDh3okqWUMfklqmSYF/7aqG1ADbb8Gbf/9wWvQ9t8fBnANGjPGL0kajCb1+CVJA9CI4I+IyyNiX0Q8EhGbq25PmSLi7IjYHREPRMT9EXFN1W2qSkSMRMTeiPhy1W0pW0SMR8StEfFQRDwYEW+ouk1li4hrO/8Hvh8RX4yIl1fdpiJFxI0RcTAivn/csdMj4q6IeLjz92n9/OzaB39EjACfBt4MnAu8KyLOrbZVpXoW+GBmngtcCLy/Zb//8a4BHqy6ERX5FPDVzFwDvJaWXYeIWAb8ETCRmecBI8A7q21V4T4PXD7r2Gbg7sw8B7i787pntQ9+4ALgkcx8NDOfAW4Grqy4TaXJzKcy87udj3/G9H/41u0fHRHLgbcAn6m6LWWLiF8Bfgv4LEBmPpOZhyttVDVOAcYi4hTgVODJittTqMz8BvCTWYevBLZ3Pt4ObOjnZzch+JcBTxz3+gAtDD6AiFgJrAW+XXFTqvBJ4EPAcxW3owqrgEPA5zpDXZ+JiFdU3agyZeYU8NfAfuAp4KeZ+bVqW1WJMzPzqc7HPwTO7OeHNCH4BUTEK4HbgD/OzP+puj1liogrgIOZuafqtlTkFOB1wA2ZuZbp7UvaNtd1GtO93VXAUuAVEfGealtVrZwuyeyrLLMJwT8FnH3c6+WdY60REaNMh/5Nmbmj6vZU4CLgrRHxGNNDfesi4gvVNqlUB4ADmTlzp3cr028EbXIp8IPMPJSZR4EdwG9U3KYq/CgizgLo/H2wnx/ShOD/DnBORKyKiJcyPaFzR8VtKk1EBNNjuw9m5vVVt6cKmXldZi7PzJVM//vfk5mt6e1l5g+BJyJidefQJcADFTapCvuBCyPi1M7/iUto2QR3xx3Axs7HG4Hb+/khtX/mbmY+GxEfAHYxPZN/Y2beX3GzynQR8F7gvoi4t3PsTzLzzuqapAr8IXBTp/PzKPAHFbenVJn57Yi4Ffgu05VuexnyVbwR8UXgjcAZEXEA+AiwBfhSRFwFPA68o6+f7cpdSWqXJgz1SJIGyOCXpJYx+CWpZQx+SWoZg1+SWsbgl3rU2TH1BxFxeuf1aZ3XKytumrQgBr/Uo8x8AriB6ZpqOn9vy8zHKmuU1APr+KU+dLbR2APcCLwPOL+zlYBUe7VfuSvVUWYejYhNwFeBNxn6ahKHeqT+vZnpLYLPq7ohUi8MfqkPEXE+cBnTT0W7dmbHRKkJDH6pR53dIW9g+tkI+4GtTD8kRGoEg1/q3fuA/Zl5V+f13wG/FhG/XWGbpAWzqkeSWsYevyS1jMEvSS1j8EtSyxj8ktQyBr8ktYzBL0ktY/BLUssY/JLUMv8PXGSx1qZ3pgkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# 生成数据\n", "data_num = 50\n", "X = np.random.rand(data_num, 1)*10\n", "Y = X * 3 + 4 + 4*np.random.randn(data_num,1)\n", "\n", "# 画出数据的分布\n", "plt.scatter(X, Y)\n", "plt.xlabel(\"X\")\n", "plt.ylabel(\"Y\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 数学原理\n", "有$N$个观测数据为:\n", "$$\n", "\\mathbf{X} = \\{x_1, x_2, ..., x_N \\} \\\\\n", "\\mathbf{Y} = \\{y_1, y_2, ..., y_N \\}\n", "$$\n", "其中$\\mathbf{X}$为自变量,$\\mathbf{Y}$为因变量。\n", "\n", "希望找到一个模型能够解释这些数据,假设使用最简单的线性模型来拟合数据:\n", "$$\n", "y = ax + b\n", "$$\n", "那么问题就变成求解参数$a$, $b$能够使得模型输出尽可能和观测数据有比较小的误差。\n", "\n", "如何构建函数来评估模型输出与观测数据之间的误差是一个关键问题,这里我们使用观测数据与模型输出的平方和来作为评估函数(也被称为损失函数Loss function):\n", "$$\n", "L = \\sum_{i=1}^{N} \\{y_i - (a x_i + b)\\}^2 \\\\\n", "L = \\sum_{i=1}^{N} (y_i - a x_i - b)^2 \n", "$$\n", "\n", "使误差函数最小,那么我们就可以求出模型的参数:\n", "$$\n", "\\frac{\\partial L}{\\partial a} = -2 \\sum_{i=1}^{N} (y_i - a x_i - b) x_i \\\\\n", "\\frac{\\partial L}{\\partial b} = -2 \\sum_{i=1}^{N} (y_i - a x_i - b)\n", "$$\n", "\n", "即当偏微分为0时,误差函数为最小,因此我们可以得到:\n", "$$\n", "-2 \\sum_{i=1}^{N} (y_i - a x_i - b) x_i = 0 \\\\\n", "-2 \\sum_{i=1}^{N} (y_i - a x_i - b) = 0 \\\\\n", "$$\n", "\n", "将上式调整一下顺序可以得到:\n", "$$\n", "a \\sum x_i^2 + b \\sum x_i = \\sum y_i x_i \\\\\n", "a \\sum x_i + b N = \\sum y_i\n", "$$\n", "\n", "上式中$\\sum x_i^2$, $\\sum x_i$, $\\sum y_i$, $\\sum y_i x_i$都是已知的数据,而参数$a$, $b$是我们想要求得未知参数。通过求解二元一次方程组,我们即可求出模型的最优参数。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.3 求解程序" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a = 3.074602, b = 3.761519\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAm80lEQVR4nO3deXiU1dnH8e8xRFlcWLXspC6ALBpARRZlUUFBWapWqrKIpaLVumHBt0pbtWJZRCuiVBFaqYKAQEVFRZSlKoIgO4gaMEEFEZAlSAjn/eMkLCHL7M/zTH6f6/JiZpjMczKYO2fuc5/7GGstIiISPCd4PQAREYmMAriISEApgIuIBJQCuIhIQCmAi4gElAK4iEhAlQn1icaYFGAJkGWt7WqMSQNeBaoAS4GbrbUHinuNqlWr2nr16kUxXBGR0mfp0qU/WGurFXw85AAO/AFYC5yad/8J4Elr7avGmOeA/sDY4l6gXr16LFmyJIxLioiIMWZTYY+HlEIxxtQCugAv5N03QAdgat5TJgLdox6liIiELNQc+GjgAeBQ3v0qwE5r7cG8+5lAzdgOTUREilNiADfGdAW2WmuXRnIBY8wAY8wSY8ySbdu2RfISIiJSiFBy4K2Ba4wxVwFlcTnwp4CKxpgyebPwWkBWYV9srR0HjANo0aLFcY1XcnJyyMzMZP/+/RF+CxILZcuWpVatWqSmpno9FBEJUYkB3Fo7BBgCYIxpB9xvrb3RGPMacC2uEqUPMDOSAWRmZnLKKadQr149XGpdEs1ay/bt28nMzCQtLc3r4YhIiKKpA/8jcK8xZiMuJ/5iJC+yf/9+qlSpouDtIWMMVapU0acgkRibsSyL1sPeJ23wbFoPe58ZywpNVEQsnDJCrLUfAB/k3f4KuDAWg1Dw9p7+DURia8ayLIZMX0l2Ti4AWTuzGTJ9JQDd02NT86GdmCIicTB8zvrDwTtfdk4uw+esj9k1FMBjrF69evzwww9RP0dEgm3LzuywHo+EAriISBzUqFgurMcjoQAOZGRk0KBBA/r27cs555zDjTfeyHvvvUfr1q05++yzWbx4MT/++CPdu3enadOmtGzZkhUrVgCwfft2rrjiCho1asStt97K0UfUvfzyy1x44YWcf/75/O53vyM3N7eoIYhIkhnUqT7lUlOOeaxcagqDOtWP2TXCWsSMu7vvhuXLY/ua558Po0eX+LSNGzfy2muvMX78eC644AL+85//sHDhQmbNmsXf/vY3ateuTXp6OjNmzOD999+nd+/eLF++nL/85S+0adOGhx9+mNmzZ/Pii64YZ+3atUyePJlFixaRmprK7bffzqRJk+jdu3dsvz8R8aX8hcrhc9azZWc2NSqWY1Cn+jFbwAS/BXAPpaWl0aRJEwAaNWpEx44dMcbQpEkTMjIy2LRpE9OmTQOgQ4cObN++nZ9++on58+czffp0ALp06UKlSpUAmDt3LkuXLuWCCy4AIDs7m9NPP92D70xEvNI9vWZMA3ZB/grgIcyU4+Wkk046fPuEE044fP+EE07g4MGDYe9QtNbSp08fHn/88ZiOU0Qkn3LgIWrbti2TJk0C4IMPPqBq1aqceuqpXHLJJfznP/8B4K233mLHjh0AdOzYkalTp7J161YAfvzxRzZtKrQjpIhIRPw1A/exP//5z9xyyy00bdqU8uXLM3HiRACGDh1Kr169aNSoEa1ataJOnToAnHvuuTz66KNcccUVHDp0iNTUVMaMGUPdunW9/DZEJImYo6sm4q1Fixa24IEOa9eupWHDhgkbgxRN/xYi/mSMWWqtbVHwcc3ARSQQZizLimtFh1+uGQ4FcBHxvUT0FfHDNcOlRUwR8b1E9BXxwzXDpQAuIr6XiL4ifrhmuJRCERHfq1GxHFmFBM7C+orEKm8dzjW9ohm4iPheqH1F8vPWWTuzsRzJW0dykEIieplESwE8DFdddRU7d+4s9jkPP/ww7733XkSv/8EHH9C1a9cSn9euXTsKlmMWNHr0aPbt2xfROET8pnt6TR7v2YSaFcthgJoVy/F4zybHzaxjmbcO9ZpeUgolBNZarLW8+eabJT73r3/9awJGVLLRo0dz0003Ub58ea+HIhITofQViXXeOt69TKIVuBl4PM6YGzVqFI0bN6Zx48aMzuvHkpGRQf369enduzeNGzfmm2++OeYghkceeYT69evTpk0bevXqxYgRIwDo27cvU6dOBdzBDUOHDqVZs2Y0adKEdevWAbB48WIuvvhi0tPTadWqFevXFz87yM7O5oYbbqBhw4b06NGD7Owj/zMOHDiQFi1a0KhRI4YOHQrA008/zZYtW2jfvj3t27cv8nkiySYRPbj9JFAz8HjUZS5dupSXXnqJTz75BGstF110EZdeeimVKlXiiy++YOLEibRs2fKYr/n000+ZNm0an3/+OTk5OTRr1ozmzZsX+vpVq1bls88+49lnn2XEiBG88MILNGjQgAULFlCmTBnee+89HnzwwcOdDgszduxYypcvz9q1a1mxYgXNmjU7/HePPfYYlStXJjc3l44dO7JixQruuusuRo0axbx586hatWqRz2vatGlE75mIX7VvUI1JH2/m6P3lfstbx1KgZuDxqMtcuHAhPXr0oEKFCpx88sn07NmTBQsWAFC3bt3jgjfAokWL6NatG2XLluWUU07h6quvLvL1e/bsCUDz5s3JyMgAYNeuXVx33XU0btyYe+65h9WrVxc7xvnz53PTTTcB0LRp02MC75QpU2jWrBnp6emsXr2aNWvWFPoaoT5PJKhmLMti2tKsY4K3AX7V3N9pkGgEKoAnui6zQoUKUb9GflvalJQUDh48CMBDDz1E+/btWbVqFf/973/Zv39/RK/99ddfM2LECObOncuKFSvo0qVLoa8V6vNEgqywCZ4F5q3b5s2A8h06BAsXxuWlAxXA45Hfatu2LTNmzGDfvn3s3buX119/nbZt2xb7Na1btz4cePfs2cMbb7wR1jV37dpFzZpuRjBhwoQSn390y9pVq1YdPs7tp59+okKFCpx22ml8//33vPXWW4e/5pRTTmH37t0lPk8kWfhu483+/fDPf8K550LbtvDZZzG/RKBy4IM61T8mBw7R57eaNWtG3759ufDCCwG49dZbSU9PP5zuKMwFF1zANddcQ9OmTTnjjDNo0qQJp512WsjXfOCBB+jTpw+PPvooXbp0KfH5AwcOpF+/fjRs2JCGDRsezrefd955pKen06BBA2rXrk3r1q0Pf82AAQPo3LkzNWrUYN68eUU+TyRZFLfxJqFNqX74AZ59FsaMga1boVkzmDQJ8k78iqXAtZP1S3ewPXv2cPLJJ7Nv3z4uueQSxo0bd8ziYhCpnawEWcEiB3ATvF81r8m0pVnHPR7zmu4NG+DJJ2HiRMjOhi5d4L77oF07MCaql06adrJ+qcscMGAAa9asYf/+/fTp0yfwwVsk6Io6RLi44oeoY4m1sGgRjBwJM2dCair07g333ONSJ3EWuADuF/k5aRHxj8ImePdMXl7oc6PKjefmwuuvw4gR8MknULky/N//wR13wC9+EfnrhskXAdxai4nyI4ZEJ5GpNJFEimlTqj174KWXXKrk66/hzDNdrrtPH4hB1Vq4PK9CKVu2LNu3b1cA8ZC1lu3bt1O2bFmvhyISczFpSvXtt/Dgg1CnDtx1F1SvDtOnw/r1cPvtngRv8MEMvFatWmRmZrJtm8e1mqVc2bJlqVWrltfDEIm5onLjIeW/V62CUaPg5Zfh4EHo2dMtTF58cZxHHRrPq1BERIriSdWZtTB3rluYfPttKF8e+vVzC5NnnhnfaxchaapQRKR0SPiZlDk5MHmyW5j8/HM44wx49FG47TaoUiWscSfql47nOXARkcIk7EzKXbtg+HBIS4Obb3aB/MUXISPDVZaEGbxjdaBEKBTARcSX4r41ftMmuPdeqFULHngA6teHN9+ElSvhllsggkX9RB+ErBSKiPhS3M6kXLrU5benTHH3b7jBLUymp0f3uiS+H4tm4CLiSzE9k/LQIXjjDWjfHlq0cLfvvhu++spVmMQgeEPiD5RQABcRX4rJmZT798MLL0CjRnD11fDll26R8ptv3J916sR0zIk+CFkpFBHxrYh7H/3wA4wdC8884zoCpqe7joDXXef6lcRJVDXnEVAAF5Hk8cUXbpv7hAmuI+BVV8H998ekI2CoEtlwr8QAbowpC8wHTsp7/lRr7VBjTBrwKlAFWArcbK09EM/Biogcx1r43/9cSiS/I+DNN7sKkwR0BPRSKDnwn4EO1trzgPOBzsaYlsATwJPW2rOAHUD/uI1SRKSg3FyYOhVatYI2bWD+fNevZNMml/dO8uANIQRw6+zJu5ua958FOgBT8x6fCHSPxwBFRI6xd6/LbZ9zjstpb93q7m/e7HZOJrCdq9dCyoEbY1JwaZKzgDHAl8BOa+3BvKdkAoUmfYwxA4ABAHVivOIrIqXIt9+6QD12LOzY4RpKDR8O3bpBSkrJX5+EQgrg1tpc4HxjTEXgdaBBqBew1o4DxoFrZhXBGEUkiYTdK2T1arfxZtIkt829Rw+38aZVq8QN2qfCqkKx1u40xswDLgYqGmPK5M3CawHx2ewvIkkj5AZV1sL777vA/dZbUK4c/Pa3bvPNWWd5MHJ/CqUKpRqQkxe8ywGX4xYw5wHX4ipR+gAz4zlQkaDwy8HbXir4HrRvUI1567YVujX+mPMpc3LcFvcRI2D5ctcR8JFHYODAsJpKlRYl9gM3xjTFLVKm4BY9p1hr/2qM+SUueFcGlgE3WWt/Lu611A9ckl1RJ6PH/AR0HyvsPSjJqT/vZUXNzfDUU5CZCQ0bujTJjTdG1FQq2UTcD9xauwI4rlGAtfYr4MLYDE8kOcT1BPSAKOw9KEqNn7bSb8ksfrPiHfh5n+tV8vzz0LkznOCK5PSJpmjaiSkSQ4nuRudHoXyvjb7byIDFr9Nl3QIAvu18DRUefQiaNTvmeQk/1CFgFMBFYihuLVADpKj3wNhDtPtqKb9d/DqtNq9g94nlmNLmV1QZfD+drir8w7w+0RRPAVwkhgZ1ql9oDjxe3ej8qOB7cNLBA3RfPY9bP53B2du/YcspVRnW4RZmX9SVzEMnUmPFHrKrZxUakPWJpngK4CIxlOhudH6U/70+P30xl304nX7L3qDynp1sqH4Wd3e9j4XpHdiVCzm5roCiuLSIPtEUTwFcJMYS2Y3OlzZupPsLT9L9pZdcR8Arr4T772dNxXP49J0N/FBSKeFR9ImmeArgIhIb+R0BZ8xwHQFvusl1BGzUKKTSwsLSIvpEUzwFcBGJXG6uC9gjR8JHH0GlSq4j4O9/f0xTqVBKC2tULFdkyaACduEUwEUkfHv3wksvucMTvvoKfvlL+Mc/oF8/qFDhuKeXtOhYLjWF9g2qqWQwTDoTU0RC99138Kc/Qe3acOedcPrprif3hg1u1l1I8IbiFx3zz7qct25bkSWDUjgFcBEp2erV0L8/1K0Lf/ubO6Js0SKXNvnVr0ps51rUYb+jf30+iwZ3oHt6TZUMRkApFBEpnLUwb55bmMzvCHjrra4j4Nlnh/VSoSxGqmQwfArgInKsgh0BTz/ddQS87TaoWjXily1pMVIlg+FTABcRZ9cud5bk6NGuI2CDBvDPf7pywAR0BFTJYPgUwEVKkULL9Koecm1cx42D3btdfvu559wGnBMSu0ymksHwKICLlBIFN9NUXLeSlJcf4dC6Ba6a4frrXQ/u5s09HaeETgFcpJQYPmc9+w/k0O6rpQw4qiPg5It70GvSSPDg0HH1+o6OArhICAIfaPbvp82HM+n/6QzO2b6Zb0+uwmPtbuHV8zux56QK9PIoePth406Q/20VwEVK8KcZK5n08WbyDx8M1A7B7dth7Fh45hme+P57Vp/+S+7ueh+zG7QhJyUVcBtpvOCHXt9++SUSKQVwkWLMWJZ1TPDO5/tDBb780m1zHz/edQTs3JmF3fry28xTyT546PDTvCzT88PGHT/8EomGArgknVh+JB4+Z/1xwTufL3cIfvSRq99+/XXXEfDGG11HwMaNaQM87qN0gR827vjhl0g0FMAlqcT6I3FxP8i+2SGYmwszZ7rAnd8RcMgQ15ukevVjnuqnMj0/bNzxwy+RaKgXiiSV4j4SR6KoH2QD3u8Q3LsXxoyB+vVdP5LvvoOnn4bNm+Gxx44L3n7TPb0mj/dsQs2K5TAcaWqVyF8wRfVo8fzfNkSagUtSifVH4sJmiQa4sWUd72ay330HzzzjFid//BEuugiGDYMePUpsKuU3Xn8iCPruTwVwSSqx/kjsqx/wNWtg1Cj4979dv5Lu3d3Gm1atwJjEjydJeP1LJBoK4JJUBnWqz6Cpnx8+MBcgNcVE9ZG4qB/wcBdLI1pctRY++MDlt9980/Uk6d8f7rkn7I6AknwUwCX5FCwbKaqMJArhLpaGvbiakwOvveYC97JlUK0a/PWvMHBgVB0BJbloEVOSyvA568k5dGzEzjlkY36qS7iLpSE//6ef3PmSZ57pSgCzs11HwM2b4aGHFLzlGJqBS1JJVF1vuNcp8fFvvnEVJOPGuSDerh08+yxcdVXCOwJKcCiAS1JJVF1vuNcp6vmX7st0/bYnT3b57uuucwuTLVrEdLwFBbn/hxyhX+2SVBJV1xvudY55vrW0+3IJr0z+Pyb84za3CefOO93291deSUjwHjJ9JVk7s7EcycfPWJYV1+tK7GkGLkklUWV/4V6ne3pNTjjwM6uGj+XaD6ZwzvbNZJ/+C/j73+G3v4WKFWM6vuIEvf+HHGGsjcMSfRFatGhhlyxZkrDrifjC9u3uhJtnnnGbcM47D+6/3x2gcOKJCR9O2uDZhRbmGODrYV0SPRwJgTFmqbX2uI9mmoGLxMuXX7rzJcePh337oHNnl9/u2NHTjTdB7/8hRygHLhJrH30E117rNto8/7ybaa9YAW+9BZdd5vmuyaD3/5AjNAOXUiOulRe5uTBrltt487//uY6Agwe7joA1asTmGjES63UCVbR4RwFcSoW4nbyybx9MmOAOT9i4EdLSXD13v35w8smHr+23ABer/h9BP9Em6JRCkVIh1m1m+f57tzOyTh244w6oXBmmTIENG1xJ4FHBO5lL9mL+vkpYFMClVIjZDs01a+DWW6FuXddzu21bWLAAPv7YbcIpc+yH2mQPcEE/0SboSgzgxpjaxph5xpg1xpjVxpg/5D1e2RjzrjHmi7w/K8V/uCKRKW6HZInyOwJ27QqNGsGkSS5Fsm6dO7qsTZsiFyaTPcBF9b5K1EKZgR8E7rPWngu0BO4wxpwLDAbmWmvPBubm3RfxpYgqL3JyjuyMbN8eFi+Gv/zFNZYaOxbOOafE6yZ7gFNFi7dKDODW2m+ttZ/l3d4NrAVqAt2AiXlPmwh0j9MYRaIW1vFdP/3kDk4480z4zW/c0WXjxsGmTfDww661a4iSPcD54Vi00iysnZjGmHrAfKAxsNlaWzHvcQPsyL9f4GsGAAMA6tSp03zTpk1RD1okVGFVgBTsCHjppW7HZJQdAaOpQvFjBYskXlE7MUMO4MaYk4EPgcestdONMTuPDtjGmB3W2mLz4NpKL4lUsMQN3Oz3uBni8uWuB/errya0I2BJQh6/JL2iAnhI0wpjTCowDZhkrZ2e9/D3xpjqeX9fHdgaq8GKxEKxFSDWwttvu52R6ekwY4bbdLNxY0I6AoYi2StYJHqhVKEY4EVgrbV21FF/NQvok3e7DzAz9sMTiVxhlR4nHsyh1fxZ0KQJXHklrF0LTzzh0idPPgn16iV+oEVI9goWiV4oOzFbAzcDK40xy/MeexAYBkwxxvQHNgHXx2WEIhE6umnTadm7uXH5W/Rd+l9O37sDmjaFf/0Lfv1rTzoChkJNp6QkJQZwa+1CXKfJwnSM7XBEYmdQp/qMefEdbvxoOtevfJfyOT+z4JfN+WLQ/bT+3a89bypVkkGd6heaA0+WChaJnnqhSHL6+GO6jxxJt+nTOWhOYGbDS5nV4df07HtVYBYAE3U4hQSXDnSQ5JHfEXDkSFi0yJ1yc9ttrjeJzzoCioRDBzpI8tq3DyZOdJtvNm50C5FPPQW33HK4qZRIMlIAl+D6/nsYMwaefdYdW3bBBa4jYI8exzWVEklG+r9cwuKLnYFr17rZ9r//DQcOwDXXuI03xTSV8gNfvHeSVBTAJWSeNu+3Fj780J14M3s2lC3rOgLec09ITaW8poMPJB7UD1xC5snOwIMH3Rb3Cy440hHwz38OqyOgH2hXpcSDZuASsoTuDNy9G154wZ3qvnmzC9TPPw833wzlgreRRbsqJR40A5eQJaS3dWYmPPAA1KoF997rKkpmzXJ57wEDAhm8Ifn7gos3FMAlZHHtbb18uZtdp6W5Ou4rr3Tpkg8/hKuvjqqdqx8ke19w8YZSKBKymO8MtBbmzHEB+733oEIF1xHwD3/wVVOpWCjuvVN1ikRKOzEl8X7+2bVsHTkSVq1yuyTvusulSCqVrqNV1fNbQhFVP3CRmPjxR3j8cTe77tfP1WxPnAhffw1//GOpC96g6hSJjlIoEn9ffeWqSV580W17v+IK18r1sst8vfEmEVSdItFQAJf4+eQTt/Fm+nRISXEHBN97r+vFLYB6fkt0lEKR2MrNdceTtW0LLVvCu+/CoEEuTTJhgoJ3AapOkWhoBi6xkd8R8Mkn4YsvXJ579GjXEfCUU7wenW+p57dEQwFcorN1q+sIOGbMkY6AkydDz57qCBii7uk1FbAlIvoJS5Ckq/Vdt851BPzXv1xZ4DXXwP33+74joEgyUQBPgKTpRGctzJ/vFibfeMN1BOzb13UErK+crUiiaREzAQJf65vfEfDCC6FdO/j44yMdAZ97TsFbxCOagSdAYGt9d+92tdujR8OmTa4j4HPPQe/egW0qlShJlzITX1IAT4DA1fpmZcHTT7v2rbt2uZLAp5+Grl0D31QqEZImZSa+p5/GBCiq1rd9g2q0HvY+aYNn03rY+8xYluXRCPN8/rmbXder5/LcnTq5zTjz57tFSgXvkAQ+ZSaBoRl4AhRW69u+QTWmLc3yfpZmLbzzjgvY+R0B77jDdQRMS0vcOJJIYFNmEjgK4AlSsNa39bD3i5ylJSSAF+wIWL06DBtWKjsCxlrgUmYSWPpMHIYZy7JilvLwbJa2Y4cL1GlpRzoCTpgAGRmltiNgrGl7vCSKZuAhivXCVMJnaV9/faQj4N69cPnlLnBffnlIG29UVRE6bY+XRFEAD1FxC1OR/GAO6lS/0Eb+MZ+lffKJS5NMm+Y6Avbq5ToCnndeyC+hqorwaXu8JIJSKCGKdcqje3pNHu/ZhJoVy2GAmhXLxe4UlkOHYObMIx0B33nnSEfAiRPDCt6gqgoRv9IMPETxSHnEfJaWnX2kI+CGDVC3rrvdv39UHQFVVSHiT5qBh8jXC1Nbt7qt7XXqwMCBcOqpbuv7xo1w991Rt3Mt6pdUwcdjucgrIiXTDDxEvlyYWr/edQScONGVBV59tesI2LZt1B0Bj160rFg+ldQTDDmHjhyAXfCXl/LkIomnAB4GXyxMWQsLFriNN//9L5x0EvTp4zoCNmgQk0sUDMY79uWQmmKoWC6VXdk5hf7yivUir4iULDABvNSXsR086CpJRoyAJUugalUYOhRuvx1OPz2mlyosGOfkWiqcVIblQ68o9GuUJxdJvEAE8FL98bxgR8Czz4axY13PkvLl43LJSIKxdh+KJJ7vFzFnLMvivimfl74ytqwsGDwYatd26ZE6ddxhwevWwW23xS14Q+iLlkfz9SKvSJLydQDPn3nnWlvo3yflx/MVK1xOOy0Nhg+HK65wByjMnw/duiWkI2AkwTiude0iUqgSUyjGmPFAV2CrtbZx3mOVgclAPSADuN5auyPWgyssF3u0pPl4bi28+67Lb7/7rusIOHCgKwH0oCNgpBU3vljkFSlFQsmBTwCeAf511GODgbnW2mHGmMF59/8Y68EVN8NOio/nBw4c6Qi4cqXrCPj44/C733neVErBWMT/Sgzg1tr5xph6BR7uBrTLuz0R+IA4BPCiFsZSjAn2x/MdO9xpN//4B2zZAo0bu8ZSN9zgygJLUOorckQEiDwHfoa19tu8298BZ8RoPMcoKhc78vrzghmwMjJcWqR2bRgyBBo1grffPpL3DjF4D5m+kqyd2ViOVORo16NI6RN1GaG11hpjCl9lBIwxA4ABAHXq1AnrtX25+zESixe7NMnUqW4RslcvuO++sJtKQXw2zGhGLxJMkQbw740x1a213xpjqgNbi3qitXYcMA6gRYsWRQb6ogQ2F3voELzxhluYXLAATjvNbXO/806oVSvil431hpmiauyXbPqReeu2KaiL+FikKZRZQJ+8232AmbEZThLIznb57YYNXdnfpk2uI+A338ATT0QVvCGyGu3iFDWjn/TxZqVpRHyuxABujHkF+Aiob4zJNMb0B4YBlxtjvgAuy7tfum3bdqQj4G23uQ6Ar7wCX34Zk46A+WK9YaaomXvBj0pJv3FKJIBCqULpVcRfdYzxWIJp/Xo3w544Efbvh65dXarkkkui7ghYmFivCxRV6VOYpNw4JRJggeiF4jv5HQFHjoRZs1z1SO/e7qiyGHUELE4s1wUKO9rNcPwMHJJo45RIkvD1VnrfOXgQpkyBiy6CSy+FRYvg4Ydh82YYNy4hwTvWCtsCf2PLOuprIhIAmoGHYvduGD/edQTMyGBPnTTGdr+b8b9sQ+VylRiUlUP32HZ0TajCZvQt6lZWaaGIzymAH6VgPfRD6afSed5UV1Wycye0acMnd/6JW7b/gr15GYdkbW0b2PJNkVJEKZQ8R+9wPGdbBvdM+hsdrroYO3w4XH656wi4YAH3Hkg7HLzzqUJDRLygGXie4W+vo/mGJQxYPJ1LMpaxL/UkJqVfyZsdrue1v990+Hk6eUZE/EIB/MABePVVXnhyKA23ZbC1QiX+fklvJp1/JbvKnULBQkCdPCMiflF6A/jOnS63/fTTsGULJ51Rj/uvuptZDS/lQJnUw08rGJgLK7tThYaIeKH0BfCMDFdN8uKLsGcPXHYZjB/PimqNmP36Kg6UEJiTpsGWiARe6Qngn37qNt689tqRjoD33gvnnw9AdwBjQgrMqtAQET9I7gB+6BDMnu06As6fD6ee6tq43nVXoU2lFJhFJEiSM4BnZ8O//w2jRrleJXXquNv9+7sgLiKSBJIrgG/bBs8+C2PGuNvNm7uOgNdeC2WS61sVEUmOqLZhg5thJ6gjoIiIHwQ3gFsLCxce6Qh44olw881uYbJhQ69HJyISd8EL4AcPwvTpLnAvXgxVqsCf/gR33AFnxOVsZRERXwpOAN+zx3UEfPJJV8t91lku392nD5Qv7/XoREQSLhgB/JFHXI57505o3drdvuYaSEkp8UtFRJJVMAL4tm1ux+R990HLll6PRkTEF4IRwJ96yhfVJAX7hWsLvYh4KRgB3CfB++gmVsl6kIOIBIcOdAjR8Dnrj+lACDrIQUS8pQAeIh3kICJ+owAeoqIObNBBDiLiFQXwEA3qVJ9yqceWLeogBxHxUjAWMX1ABzmIiN8ogIdB/cJFxE+UQhERCSgFcBGRgFIAFxEJKAVwEZGAUgAXEQkoBXARkYBSABcRCSgFcBGRgFIAFxEJKAVwEZGAUgAXEQkoBXARkYCKKoAbYzobY9YbYzYaYwbHalAiIlKyiAO4MSYFGANcCZwL9DLGnBurgYmISPGimYFfCGy01n5lrT0AvAp0i82wRESkJNEE8JrAN0fdz8x77BjGmAHGmCXGmCXbtm2L4nIiInK0uC9iWmvHWWtbWGtbVKtWLd6XExEpNaIJ4FlA7aPu18p7TEREEiCaAP4pcLYxJs0YcyJwAzArNsMSEZGSRHwmprX2oDHm98AcIAUYb61dHbORiYhIsaI61Nha+ybwZozGIiIiYdBOTBGRgFIAFxEJKAVwEZGAUgAXEQkoBXARkYBSABcRCaioygj9aMayLIbPWc+WndnUqFiOQZ3q0z39uBYtIiKBl1QBfMayLIZMX0l2Ti4AWTuzGTJ9JYCCuIgknaRKoQyfs/5w8M6XnZPL8DnrPRqRiEj8JFUA37IzO6zHRUSCLKkCeI2K5cJ6XEQkyJIqgA/qVJ9yqSnHPFYuNYVBnep7NCIRkfhJqkXM/IVKVaGISGmQVAEcXBBXwBaR0iCpUigiIqWJAriISEApgIuIBJQCuIhIQCmAi4gElLHWJu5ixmwD9gI/JOyi/lSV0v0elPbvH/QelPbvH8J7D+paa6sVfDChARzAGLPEWtsioRf1mdL+HpT27x/0HpT27x9i8x4ohSIiElAK4CIiAeVFAB/nwTX9prS/B6X9+we9B6X9+4cYvAcJz4GLiEhsKIUiIhJQCQ3gxpjOxpj1xpiNxpjBiby214wxtY0x84wxa4wxq40xf/B6TF4xxqQYY5YZY97weiyJZoypaIyZaoxZZ4xZa4y52OsxJZox5p68n4FVxphXjDFlvR5TPBljxhtjthpjVh31WGVjzLvGmC/y/qwUyWsnLIAbY1KAMcCVwLlAL2PMuYm6vg8cBO6z1p4LtATuKGXf/9H+AKz1ehAeeQp421rbADiPUvY+GGNqAncBLay1jYEU4AZvRxV3E4DOBR4bDMy11p4NzM27H7ZEzsAvBDZaa7+y1h4AXgW6JfD6nrLWfmut/Szv9m7cD26p63trjKkFdAFe8HosiWaMOQ24BHgRwFp7wFq709NBeaMMUM4YUwYoD2zxeDxxZa2dD/xY4OFuwMS82xOB7pG8diIDeE3gm6PuZ1IKAxiAMaYekA584vFQvDAaeAA45PE4vJAGbANeykshvWCMqeD1oBLJWpsFjAA2A98Cu6y173g7Kk+cYa39Nu/2d8AZkbyIFjETzBhzMjANuNta+5PX40kkY0xXYKu1dqnXY/FIGaAZMNZam45rK1Ha1oIq4WafaUANoIIx5iZvR+Ut60oBIyoHTGQAzwJqH3W/Vt5jpYYxJhUXvCdZa6d7PR4PtAauMcZk4FJoHYwxL3s7pITKBDKttfmfvKbiAnppchnwtbV2m7U2B5gOtPJ4TF743hhTHSDvz62RvEgiA/inwNnGmDRjzIm4hYtZCby+p4wxBpf7XGutHeX1eLxgrR1ira1lra2H+/d/31pbamZf1trvgG+MMfmnbHcE1ng4JC9sBloaY8rn/Ux0pJQt5OaZBfTJu90HmBnJiyTsTExr7UFjzO+BObiV5/HW2tWJur4PtAZuBlYaY5bnPfagtfZN74YkHrgTmJQ3ifkK6OfxeBLKWvuJMWYq8BmuMmsZSb4r0xjzCtAOqGqMyQSGAsOAKcaY/sAm4PqIXls7MUVEgkmLmCIiAaUALiISUArgIiIBpQAuIhJQCuAiIgGlAC4iElAK4CIiAaUALiISUP8PT1ZrL1HOVhoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "N = X.shape[0]\n", "\n", "S_X2 = np.sum(X*X)\n", "S_X = np.sum(X)\n", "S_XY = np.sum(X*Y)\n", "S_Y = np.sum(Y)\n", "\n", "A1 = np.array([[S_X2, S_X], \n", " [S_X, N]])\n", "B1 = np.array([S_XY, S_Y])\n", "\n", "# numpy.linalg模块包含线性代数的函数。\n", "# 使用这个模块,可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等。\n", "coeff = np.linalg.inv(A1).dot(B1)\n", "\n", "print('a = %f, b = %f' % (coeff[0], coeff[1]))\n", "\n", "x_min = np.min(X)\n", "x_max = np.max(X)\n", "y_min = coeff[0] * x_min + coeff[1]\n", "y_max = coeff[0] * x_max + coeff[1]\n", "\n", "plt.scatter(X, Y, label='original data')\n", "plt.plot([x_min, x_max], [y_min, y_max], 'r', label='model')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. 如何使用迭代的方法求出模型参数\n", "\n", "当数据比较多的时候,或者模型比较复杂,无法直接使用解析的方式求出模型参数。因此更为常用的方式是,通过迭代的方式逐步逼近模型的参数。\n", "\n", "### 2.1 梯度下降法\n", "在机器学习算法中,对于很多监督学习模型,需要对原始的模型构建损失函数,接下来便是通过优化算法对损失函数进行优化,以便寻找到最优的参数。在求解机器学习参数的优化算法中,使用较多的是基于梯度下降的优化算法(Gradient Descent, GD)。\n", "\n", "梯度下降法有很多优点,其中最主要的优点是,**在梯度下降法的求解过程中只需求解损失函数的一阶导数,计算的代价比较小,这使得梯度下降法能在很多大规模数据集上得到应用。**\n", "\n", "梯度下降法的含义是通过当前点的梯度方向寻找到新的迭代点。梯度下降法的基本思想可以类比为一个下山的过程。假设这样一个场景:\n", "* 一个人被困在山上,需要从山上下来,找到山的最低点,也就是山谷;\n", "* 但此时山上的浓雾很大,导致可视度很低。因此,下山的路径就无法全部确定,他必须利用自己周围的信息去找到下山的路径。\n", "* 以他当前的所处的位置为基准,寻找这个位置最陡峭的地方,然后朝着山的高度下降的地方走\n", "* 每走一段距离,都反复采用同一个方法,最后就能成功的抵达山谷。\n", "\n", "\n", "一般情况下,这座山最陡峭的地方是无法通过肉眼立马观察出来的,而是需要一个工具来测量;同时,这个人此时正好拥有测量出最陡峭方向的能力。所以,此人每走一段距离,都需要一段时间来测量所在位置最陡峭的方向,这是比较耗时的。那么为了在太阳下山之前到达山底,就要尽可能的减少测量方向的次数。这是一个两难的选择,如果测量的频繁,可以保证下山的方向是绝对正确的,但又非常耗时;如果测量的过少,又有偏离轨道的风险。所以需要找到一个合适的测量方向的频率,来确保下山的方向不错误,同时又不至于耗时太多!\n", "\n", "\n", "![gradient_descent](images/gradient_descent.png)\n", "\n", "如上图所示,得到了最优解。$x$,$y$表示的是$\\theta_0$和$\\theta_1$,$z$方向表示的是花费函数,很明显出发点不同,最后到达的收敛点可能不一样。当然如果是碗状的,那么收敛点就应该是一样的。\n", "\n", "对于最小二乘的损失函数\n", "$$\n", "L = \\sum_{i=1}^{N} (y_i - a x_i - b)^2\n", "$$\n", "\n", "更新的策略是:\n", "$$\n", "\\theta^1 = \\theta^0 - \\eta \\triangledown L(\\theta)\n", "$$\n", "其中$\\theta$代表了模型中的参数,例如$a$, $b$\n", "\n", "此公式的意义是:$L$是关于$\\theta$的一个函数,我们当前所处的位置为$\\theta_0$点,要从这个点走到L的最小值点,也就是山底。首先我们先确定前进的方向,也就是梯度的反向,然后走一段距离的步长,也就是$\\eta$,走完这个段步长,就到达了$\\theta_1$这个点!\n", "\n", "最终的更新方程是:\n", "\n", "$$\n", "a^1 = a^0 + 2 \\eta [ y - (ax+b)]*x \\\\\n", "b^1 = b^0 + 2 \\eta [ y - (ax+b)] \n", "$$\n", "\n", "下面就这个公式的几个常见的疑问:\n", "\n", "* **$\\eta$是什么含义?**\n", "$\\eta$在梯度下降算法中被称作为学习率或者步长,意味着我们可以通过$\\eta$来控制每一步走的距离,以保证不要步子跨的太大,错过了最低点。同时也要保证不要走的太慢,导致太阳下山了,还没有走到山下。所以$\\eta$的选择在梯度下降法中往往是很重要的。\n", "![gd_stepsize](images/gd_stepsize.png)\n", "\n", "* **为什么要梯度要乘以一个负号?**\n", "梯度前加一个负号,就意味着朝着梯度相反的方向前进!梯度的方向实际就是函数在此点上升最快的方向,而我们需要朝着下降最快的方向走,自然就是负的梯度的方向,所以此处需要加上负号。\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 示例代码" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "epoch 0: loss = 964.462432, a = 3.605114, b = 1.371685\n", "epoch 100: loss = 856.680487, a = 3.044999, b = 3.374801\n", "epoch 200: loss = 900.182987, a = 3.261635, b = 3.670704\n", "epoch 300: loss = 841.107710, a = 3.082096, b = 3.703232\n", "epoch 400: loss = 848.034818, a = 3.023324, b = 3.703601\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD5CAYAAAA+0W6bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAl6ElEQVR4nO3deXxV1bn/8c8SUplUENCrDEIVwyzBQFEERVrRAhrRVrkytVqq1lsHxIptpVpbaRm1xYHigIoiAkbqcPmVwaL2KrPMKCJDgsqgAZSgAdbvj3VCIJwkZ957n/N9v16+yNmcnL3OwTxZe+3nWY+x1iIiIsFzgtcDEBGR2CiAi4gElAK4iEhAKYCLiASUAriISEApgIuIBFT1SJ9ojKkGLAEKrbV9jDHNgWlAfWApMNBa+11lr9GgQQPbrFmzOIYrIpJ5li5dusta27D88YgDOHA7sA44OfT4L8B4a+00Y8wTwI3A45W9QLNmzViyZEkUpxQREWPMlnDHI1pCMcY0BnoDk0OPDXApMCP0lClAXtyjFBGRiEW6Bj4BuAc4HHpcHyiy1h4MPS4AGiV2aCIiUpkqA7gxpg+ww1q7NJYTGGOGGmOWGGOW7Ny5M5aXEBGRMCJZA+8KXGmM+TFQA7cG/ghQ1xhTPTQLbwwUhvtma+0kYBJAbm7ucRuvlJSUUFBQwIEDB2J8C5IINWrUoHHjxmRlZXk9FBGJUJUB3Fo7AhgBYIy5BLjbWnuDMeYV4FpcJspg4LVYBlBQUMBJJ51Es2bNcEvrkmrWWnbv3k1BQQHNmzf3ejgiEqF48sB/A9xljNmIWxN/KpYXOXDgAPXr11fw9pAxhvr16+sqSCTB8pcX0nXUfJrf+wZdR80nf3nYhYqYRZNGiLX2beDt0NebgM6JGISCt/f0byCSWPnLCxkxaxXFJYcAKCwqZsSsVQDk5SQm50OVmCIiSTB6zoYjwbtUcckhRs/ZkLBzKIAnWLNmzdi1a1fczxGRYNteVBzV8VgogIuIJMGZdWtGdTwWCuDA5s2badmyJUOGDOHcc8/lhhtuYO7cuXTt2pUWLVqwaNEivvzyS/Ly8mjfvj1dunRh5cqVAOzevZvLLruMNm3acNNNN3F0i7oXXniBzp0706FDB375y19y6NChioYgImlmeK9samZVO+ZYzaxqDO+VnbBzRHUTM+nuuANWrEjsa3boABMmVPm0jRs38sorr/D000/TqVMnXnzxRd59911mz57Nn//8Z5o0aUJOTg75+fnMnz+fQYMGsWLFCh544AEuuugi7r//ft544w2eesol46xbt46XX36Z9957j6ysLG699VamTp3KoEGDEvv+RMSXSm9Ujp6zge1FxZxZtybDe2Un7AYm+C2Ae6h58+a0a9cOgDZt2tCzZ0+MMbRr147NmzezZcsWZs6cCcCll17K7t272bt3LwsXLmTWrFkA9O7dm3r16gEwb948li5dSqdOnQAoLi7mtNNO8+CdiYhX8nIaJTRgl+evAB7BTDlZTjzxxCNfn3DCCUcen3DCCRw8eDDqCkVrLYMHD+bhhx9O6DhFREppDTxC3bp1Y+rUqQC8/fbbNGjQgJNPPpnu3bvz4osvAvDWW2/x1VdfAdCzZ09mzJjBjh07APjyyy/ZsiXsjpAiIjHx1wzcx/7whz/w85//nPbt21OrVi2mTJkCwMiRI+nfvz9t2rThwgsvpGnTpgC0bt2ahx56iMsuu4zDhw+TlZXFxIkTOeuss7x8GyKSRszRWRPJlpuba8s3dFi3bh2tWrVK2RikYvq3EPEnY8xSa21u+eOagYtIIOQvL0xqRodfzhkNBXAR8b1U7Cvih3NGSzcxRcT3UrGviB/OGS0FcBHxvVTsK+KHc0ZLSygi4ntn1q1JYZjAGW5fkUStW0dzTq9oBi4ivhfpviKl69aFRcVYytatY2mkkIq9TOKlAB6FH//4xxQVFVX6nPvvv5+5c+fG9Ppvv/02ffr0qfJ5l1xyCeXTMcubMGEC+/fvj2kcIn6Tl9OIh/u1o1HdmhigUd2aPNyv3XEz60SuW0d6Ti9pCSUC1lqstbz55ptVPvfBBx9MwYiqNmHCBAYMGECtWrW8HopIQkSyr0ii162TvZdJvAI3A09Gj7lx48bRtm1b2rZty4TQfiybN28mOzubQYMG0bZtW7Zt23ZMI4Y//vGPZGdnc9FFF9G/f3/GjBkDwJAhQ5gxYwbgGjeMHDmSjh070q5dO9avXw/AokWLuOCCC8jJyeHCCy9kw4bKZwfFxcVcf/31tGrViquvvpri4rL/GW+55RZyc3Np06YNI0eOBODRRx9l+/bt9OjRgx49elT4PJF0k4o9uP0kUDPwZORlLl26lGeeeYYPPvgAay0/+MEPuPjii6lXrx4ff/wxU6ZMoUuXLsd8z+LFi5k5cyYffvghJSUldOzYkfPPPz/s6zdo0IBly5bx2GOPMWbMGCZPnkzLli155513qF69OnPnzuW+++47stNhOI8//ji1atVi3bp1rFy5ko4dOx75uz/96U+ceuqpHDp0iJ49e7Jy5Up+/etfM27cOBYsWECDBg0qfF779u1j+sxE/KpHy4ZMfX8rR9eX+23dOpECNQNPRl7mu+++y9VXX03t2rWpU6cO/fr145133gHgrLPOOi54A7z33ntcddVV1KhRg5NOOom+fftW+Pr9+vUD4Pzzz2fz5s0A7Nmzh5/85Ce0bduWO++8kzVr1lQ6xoULFzJgwAAA2rdvf0zgnT59Oh07diQnJ4c1a9awdu3asK8R6fNEgip/eSEzlxYeE7wNcM35/l4GiUegAniq8zJr164d92uUbktbrVo1Dh48CMDvf/97evTowerVq/nnP//JgQMHYnrtTz/9lDFjxjBv3jxWrlxJ7969w75WpM8TCbJwEzwLLFi/05sBpUCgAngy1re6detGfn4++/fv55tvvuHVV1+lW7dulX5P165djwTer7/+mtdffz2qc+7Zs4dGjdyM4Nlnn63y+UdvWbt69eoj7dz27t1L7dq1OeWUU/jiiy946623jnzPSSedxL59+6p8nki68G3hzdq1MHw4hCZwiRSoNfDhvbKPWQOH+Ne3OnbsyJAhQ+jcuTMAN910Ezk5OUeWO8Lp1KkTV155Je3bt+f000+nXbt2nHLKKRGf85577mHw4ME89NBD9O7du8rn33LLLfzsZz+jVatWtGrV6sh6+3nnnUdOTg4tW7akSZMmdO3a9cj3DB06lMsvv5wzzzyTBQsWVPg8kXRRWeFNyjelshYWLICxY+HNN6FGDbj+eqjgXlmsAredrF92B/v666+pU6cO+/fvp3v37kyaNOmYm4tBpO1kJcjKJzmAm+Bdc34jZi4tPO54UnK6S0pg+nQYM8b19z3tNLjtNrjlFgglFMQibbaT9Ute5tChQ1m7di0HDhxg8ODBgQ/eIkFXURPhypIfEhZL9uyBSZPg0UehoABatYJ//AMGDHCz7yQJXAD3i9I1aRHxj3ATvDtfXhH2uQlZG9+8GR55BCZPhq+/hh494Mkn4fLL4YTk32L0RQC31mKM8XoYGS2VS2kiqZSUTakWL3br2zNmgDFw3XVw112Q4itxz7NQatSowe7duxVAPGStZffu3dRI4qWeiFcStinV4cPw2mvQvTt07gxvveWC9qZN8MILKQ/e4IMZeOPGjSkoKGDnzvTN1QyCGjVq0LhxY6+HIZJwFa2NR7z+vX8/PPccjB8PH30ETZvCuHFw441w8slJHHnVPA/gWVlZNG/e3OthiIgPJSrrLKbkhy++gIkT4bHHYPduyM2FadPgmmuguuehE/BBABcRCceznpTr1rkZ9vPPw3ffQd++MGwYdOvm1rsjGHeqUp09XwMXEQknpT0prYX586F3b2jd2q1pDxnignnpuneEwTtRDSUioQAuIr6UktL4khKYOtVVSPbs6bJLHngAtm6FJ56A7OhudKa6EbKWUETEl5LakzJJhTep3o9FM3AR8aWk9KTcssWl/jVpAvfcAy1awOuvw+rVcNNNcVdNprqhhAK4iPhSQntSLl7sNpM6+2w36+7bF5YuLVv3TlDVZKobIWsJRUR8K669jw4fdrPrsWNh4UKXs33nnfDrX7sZeBLEnXMeJQVwEUkvxcWu8GbcOE8Kb1K54V6VAdwYUwNYCJwYev4Ma+1IY0xzYBpQH1gKDLTWfpfMwYqIVGjHjrLCm127XOHNSy/Btdf6pvAm0SJZ+PkWuNRaex7QAbjcGNMF+Asw3lp7DvAVcGPSRikiUpF16+AXv3Az7QcfhAsugH//GxYtcuveaRq8IYIAbp2vQw+zQv9Z4FJgRuj4FCAvGQMUETlOacebPn2OLbxZvx5mz4648CboIvrVZIyphlsmOQeYCHwCFFlrS5u8FQBhF32MMUOBoQBNmzaNd7wikslKSuCVV9yNyWXLoGFDV3hzyy3u6wwTUQC31h4COhhj6gKvAi0jPYG1dhIwCVxLtRjGKCJpJKa9QvbscYU2jzziCm9atnSFOAMGQM3k5FgHQVSLQ9baImPMAuACoK4xpnpoFt4YSE6xv4ikjag3qNqypazjzb59cMklrsT9iitS0vHG7yLJQmkIlISCd03gR7gbmAuAa3GZKIOB15I5UJGg8EvjbS+V/wx6tGzIgvU7w5bGh+1PuWSJWyZ55RX3+Lrr3I6A6j17jEhm4GcAU0Lr4CcA0621rxtj1gLTjDEPAcuBp5I4TpFA8GwLVB8J9xm88P7WSr9ne1GxK7x54w3X0T1FhTdBV2UAt9auBHLCHN8EdE7GoESCKiUd0H0u3GdQmRNLvuXGTe9A62GwYYNLBxw71u1NcvLJbjY/dX5GX9FUJH0TJEU8kOrd6Pwo0vda/5siBi17g4Er3uDU/Xvdlq4vvugKb7KyAF3RVEUBXCSBkroFakBU9BmUOnv3Nm5cnM81q+dz4qESPuv+I3jwt2Fzt3VFUzndxhVJoFTvRudH4T4DrKXL1pU8NeMB5k2+hX5rFvDPnB/R86YnuPaKe8k/+ZywhTe6oqmcZuAiCZTq3ej86OjPYMfufdywbRG/WPwqjTatY1etU3j8koE8e94VfFEjtLFUJcsiuqKpnAK4SIKlcjc6v8r7fh3yqi+Dlx+FbdsgO5vlvx3FXSe259NvDh/3/IqWRYb3yj5mDRwy74qmMgrgIpI4W7e6wpt//KOs8Oaxx8g/4zxG5K+h+JuKs1PCLYvoiqZyCuAiEr9whTd33eUyS4DRo+ZXmVp4Zt2aFRZBKWCHpwAuIrEpLbwZO9Zt33rSSXDHHa7wptzGdVXddKyZVY0eLRsqZTBKykIRkegUF8OTT7ptXK+8EjZtckG8oMBVUYbZdbSym46lvS4XrN9ZYcqghKcALiKR2bED/vAHF6Bvvhnq1HGFN5984pZLKmlXVlF65YTrOvDevZeSl9NIKYMx0BKKiFRu/XrXU/K55+Dbb10ThbvvjqppQiQ3I5UyGD0FcBE5nrVuXXvsWNfZvUYNGDzYbS7VMuJ2AMeo6makUgajpwAuImXKd7xp0MAtm9x6a9I73ihlMHoK4CIZpMK9yvfscU0THnnkSOENTz4JAwemtOONUgajowAukiHC7ez36DPzaDvuPc557SVXeHPxxfDYY/DjH6vjTQAogItkiKN39mv7+UZ+sehVeq9/BwxlHW9ChTepou5F8VEAF4lAOgSaz776hks/WcLQRbPosm01+75Xk6dzr2JKbl/ee+xnKR+PX/b6DvK/rQK4SBV+l7+Kqe9vxYYeB65CsLgYnn+eBc88xFk7t1F4UkMe6vFzXj6vF/tOrE0jj9L0/LDXt19+icRKAVykEvnLC48J3qUC0VRg5063nj1xIuzcySkt2zGs63/z2jkXcLCa+9H3Mk3PD4U7fvglEg8FcEk7ibwkHj1nw3HBu5RvKwQ3bCgrvDlwwBXeDBtG3YsvptuK7bzvk+UCPxTu+OGXSDwUwCWtJPqSuLIfZF9VCFrrOrmPGeMKb048EQYNciXuRxXe+ClNzw+FO374JRIP5QlJWqnskjgWFf0gG/BHhWBJCbz0EnTq5Pbefv99GDnS7cs9aVLMVZOpkJfTiIf7taNR3ZoYyja1SuUvmKC3wNMMXNJKoi+Jw80SDXBDl6bezmT37nVNE0oLb84915PCm3h5fUUQ9OpPBXBJK4m+JPbdD/i2bWUdb/budYU3EydC794qvImR179E4qEALmlleK9shs/4kJJDZbces6qZuC6JK/oBj/ZmaVw3V5cudfuTTJ/uHv/kJ67wJjc3lrckaUIBXNJP+bSRitJI4hDtzdKYbq4ePgxvvukC99tvu443t9/uOt6cdVbi35QEjq65JK2MnrOBksPHRuySwzbhXV2ivVka1fOLi90NyDZtoG9f2LgRRo92yydjxyp4yxGagUtaSVVeb7Tnieh4ucIbOnaEqVPdcklWVtxjlvSjAC5pJVV5vdGep9LnV1B4w8UXR9zxJlpB3v9DymgJRdJKqvJ6oz3Pcc+3lu7b1zDjrVEuV3vKFJcCuHYt/POfLqc7icF7xKxVFBYVYylbj89fXpiU80nyaAYuaSVVaX/Rnqf0+Lg319Dhg7ncuuw1WhZ+5Dre3H8//OpXcNppCR1jRYK+/4eUUQCXtJOqvN6ozrN3L3kLXiZv0iOuSvLcc+GJJ1y5e4oLb4K+/4eUUQAXSaZt2+DRR11WSWnhzd//7mnhTdD3/5AyWgMXSYZly+CGG+D734fx412LssWLXT53376eVk0Gff8PKaMZuGSMpGdehCu8+Z//ccU3PsrdTvR9AmW0eMdYm4QytQrk5ubaJUuWpOx8IqXKV0KCm3UmZPe7Awfg+eddKuD69dC4sQvav/gFnHJKWge4pH6ucoQxZqm19rh9E7SEIhkh0dvMAq7Y5sEHoWlTGDrU3YycOhU2bYK77z4SvNM5ZS8pn6tETAFcMkJCMy82bICbb3aBe+RI6NwZ5s93G079938fUzWZ7gFOGS3eqjKAG2OaGGMWGGPWGmPWGGNuDx0/1RjzL2PMx6E/6yV/uCKxqaxCMiLWwr//DVde6Qpvnn0WBgxwhTevvw49eoQtvEn3ABf35ypxiWQGfhAYZq1tDXQBfmWMaQ3cC8yz1rYA5oUei/hSzJkXBw/CtGluln3JJfB//+cKb7ZscXtyt2pV6bene4BTRou3qgzg1trPrLXLQl/vA9YBjYCrgCmhp00B8pI0RpG4Rd2+a+9el/539tnQv797/MQTrgjngQfg9NMjOm+6Bzg/tEXLZFFloRhjmgELgbbAVmtt3dBxA3xV+rjc9wwFhgI0bdr0/C1btsQ9aJFIRZ0BUr7wpnt3t7FUnz4x527Hk4WSzhksErmKslAiDuDGmDrAv4E/WWtnGWOKjg7YxpivrLWVroMrjVBSKaoUt2XLyjreWAvXXusCd6dOKR51GaXoSam40giNMVnATGCqtXZW6PAXxpgzQn9/BrAjUYMVSYQqM0AOH4Y33oBLL4Xzz4fZs13hzSefuHVvD4M3pH8Gi8QvkiwUAzwFrLPWjjvqr2YDg0NfDwZeS/zwRGJXUabHrl17YPJk1/GmTx/4+GPX8aagwBXj+KRqMt0zWCR+kZTSdwUGAquMMStCx+4DRgHTjTE3AluAnyZlhCIxKr9pU739exi4/E2GrHgDvi6CnBx44QX46U992fFGm05JVaoM4Nbad4GKdpbvmdjhiCTO8F7ZjJi1ijM+38KNS/K5ZvV8ahz8js8vuhQe+G2Fudt+UTr+8mvg6ZLBIvHTZlaSnqwlb98n/OCd8Zy+8F+UnFCdOR1/RO3f3E3Pa3p4PbqIpKo5hQSXNrOS9HLwIMycCWPGwJIlUL++63Zz660R526L+E1FWSiagUt62LfP3Zh85BFXJdmiBTz+uOt4U6uW16MTSQoFcAm28oU33bq5x3EU3ogEhQK4RMU3lYHLl7vCm5df9k3hTVV889lJ2lAAl4iVrwws3dsaSE0gOnwY/vd/3fr2ggVQpw7cdptrntCsWfLPHwfPPztJS7rGlIh5Vhl44IBb327b1jUD/ugj+Otf3fLJ+PG+D96gqkpJDs3AJWIprwzctQseewwmToQdO6BDB18X3lRGVZWSDArgErGUVQZ+9JGbWU+ZAsXFrqP7sGG+L7ypjKoqJRm0hCIRS+re1tbCO+/AVVe5jjdPP+3ak61ZU7bhVECDN6T/vuDiDc3AJWJJqQwsLbwZOxYWL3aFN7/7nSu+SaPCm8o+O2WnSKxUiSne2LcPnnoKJkwoK7y5666MK7zRnt8Sibj2AxdJmIICuOceaNIE7rzTdXbPz4f1612n9wwK3qDsFImPllAkNY4uvDl8uKzwpnNnr0fmKWWnSDwUwCV5Sgtvxo6F+fMDVXiTKspOkXhoCUUS78ABt75dWnizYUPgCm9SRdkpEg/NwCVxdu1yOwD+/e9lhTfPP+8Kb773Pa9H50va81vioQAu8StfeHPFFXD33YEuvEmlvJxGCtgSEwXwFEm7XF9r4d133fr27NmutH3gQJcK2Lq116MTyQgK4CmQVjvRHTwIs2a5HQEXL4ZTT4Xf/tYV3vzXf3k9OpGMopuYKZAWub779rmim3POgeuug6Iit9HUtm3wxz8qeIt4QDPwFAh0rm9BQVnHmz174KKLXCDv2xeqVavy2zNV2i2ZiS8pgKdAIHN9V6xw69vTpqnwJkpptWQmvqYllBSoKNe3R8uGdB01n+b3vkHXUfPJX17o0QhDrIW33oKePSEnx5W433YbfPKJq6BU8I5IWiyZSSBoBp4C4XJ9e7RsyMylhf6YpR04AFOnwrhxsHYtNGoEf/kLDB0KdeumdixpINBLZhIoCuApUj7Xt+uo+RXO0lIWwMsX3px3Hjz3nLtJqcKbmAVyyUwCSQE8Com8MeXpLO3jj13hzbPPlhXeDBsW+KYJfjG8V3bYLWJVHi+JpgAeoUTfmEr5LM1aeO89d2PytdfKCm/uvBPatKny25VVETmVx0uqKIBHqLIbU7H8YKZsllZaeDN2LCxaFFPhjbIqoqfyeEkFZaFEKNFLHnk5jXi4Xzsa1a2JARrVrZnYLiylhTctWrg17S+/dN3dt26NuvBGWRUi/qQZeISSseSRlFlaYaErvHnyybLCm/Hj4yq8UVaFiD9pBh4h3+/b/OGHrp9ks2Zun5LLLoP333ed3vPy4qqarOiXVPnj+csL/ZXXLpLmNAOPkC9vTFlb1vFm3jyoXdutbd9+OzRvHtdLH33Tsm6tLLJOMJQcLmuAXf6Xl9bJRVJPATwKvrkxVb7w5swzYdQoV3hTr17cL18+GH+1v4Ssaoa6NbPYU1wS9pdXom/yikjVAhPAlcYG7N5dVnjzxRdJK7wJF4xLDllqn1idFSMvC/s9WicXSb1ABPCMvzwvX3hz+eWu402SCm9iCcaqPhRJPd/fxMxfXsiw6R9mXhpbacebq6+G7GzXJLh/f1i9umzDqSRVTUZ60/Jovr/JK5KGfB3AS2feh6wN+/dpeXl+8CC88gp06QLdusHChXDffbBliwviEVRNxiuWYJz0vHYROU6VSyjGmKeBPsAOa23b0LFTgZeBZsBm4KfW2q8SPbhwa7FHS6vL83374OmnXfHN5s2u883EiTB4sMsuSaFYM258c5NXJENEsgb+LPB34Lmjjt0LzLPWjjLG3Bt6/JtED66yGXbaXJ4XFsLf/gZPPOEKb7p2jbvwJhEUjEX8r8oAbq1daIxpVu7wVcAloa+nAG+ThABe0Y2xasYE//L8ww9d/vZLL7mON9dc43YE/MEPqvxWZeSICMS+Bn66tfaz0NefA6cnaDzHqGgtduxPzwtmwCotvPnhD6FDB7fJ1K23wsaNMH16xMF7xKxVFBYVYynLyFHVo0jmiTuN0FprjTHh7zICxpihwFCApk2bRvXavqx+jMW335YV3qxZE1fhTTIKZjSjFwmmWAP4F8aYM6y1nxljzgB2VPREa+0kYBJAbm5uhYG+IoFei929261t/+1vrvCmffu4C28SXTBTUY79ki1fsmD9TgV1ER+LdQllNjA49PVg4LXEDCdNbNzo9iRp0gR+9zvXIPhf/3Kd3gcOjKtqMpYc7cpUNKOf+v5WLdOI+FyVAdwY8xLwf0C2MabAGHMjMAr4kTHmY+CHoceZrbTjzdVXw7nnwuTJcP31sGqVK7z54Q8TUniT6IKZimbu5S+V0r5wSiSAIslC6V/BX/VM8FiC6eBBePVVl1HywQeu481998Ftt0XVNCFSib4vUFGmTzhpWTglEmCB2AvFl77+2hXejB/vCm/OPtttMjVkSNILbxJ5XyBcazfD8TNwSLPCKZE0oAAerdLCmyefhKIiV3gzbhxceaWnhTexCjej79GyITOXFqqruojPKYBHqrTwZto07KFDLGh9EX/rcyU72uQwvFk2eQEM3qXCzehzzzpVqYUiPqcAfpTj8qEvO5e8Hatd4J47F2rX5pNrB3LzqRfxcZ3T3Del6da2gU7fFMkQvt6NMJWOrnDMOljChQtn0+byi+CKK1zXm4cfhm3bGNT+hrLgHaIMDRHxgmbgIaPnbODEvV/x8+VvMXjZ65z2zVesa9iMP/7kN/z+hQeP5G6r84yI+IUCOMDGjfzylXFcu2outUq+5d/NO3Jn52G8d9Z5GGP4/VGFN+o8IyJ+kbkB3Fr4z3/c+nZ+PtefUJ381hczuVMeHzVsduRp5QNzuLQ7ZWiIiBcyL4CXL7ypVw9GjGDexdcw8p0dVQbmtNlgS0QCL3MCeGnhzYQJ8OmnxxXeXAF82zCyXfmUoSEifpD+AXz79rKON0VFcOGFbvYdpvBGgVlEgiR9A/jKlWUdbw4dgn79XMebLl28HpmISEKkVwC3FubMOabwhptvhjvugO9/3+vRiYgkVHoE8G+/hRdfdHuSrF4NZ5zhCm9++cuoO96IiARFsAP4l1+Wdbz5/HNo1w6mTHH7cMfRNEFEJAiCGcA/+cRt4/rMM7B/P/Tq5VqVJahpgohIEAQrgP/nPzBmDOTnQ/XqcMMNcNddbuYtIpJhghHAX30V/vpXeP/9I4U33HabW+sWEclQwQjgU6fCjh0p63gjIhIEwQjgTz4Jdet63vHmuP3CVUIvIh4KRgCvX9/rERzZL7x0r5TCNG3kICLBoYYOERo9Z8MxG12BGjmIiLcUwCOkRg4i4jcK4BGqqGGDGjmIiFcUwCM0vFc2NbOOvYmqRg4i4qVg3MT0ATVyEBG/UQCPgvYLFxE/0RKKiEhAKYCLiASUAriISEApgIuIBJQCuIhIQCmAi4gElAK4iEhAKYCLiASUAriISEApgIuIBJQCuIhIQCmAi4gEVFwB3BhzuTFmgzFmozHm3kQNSkREqhZzADfGVAMmAlcArYH+xpjWiRqYiIhULp4ZeGdgo7V2k7X2O2AacFVihiUiIlWJJ4A3ArYd9bggdOwYxpihxpglxpglO3fujON0IiJytKTfxLTWTrLW5lprcxs2bJjs04mIZIx4Angh0OSox41Dx0REJAXiCeCLgRbGmObGmO8B1wOzEzMsERGpSsw9Ma21B40xtwFzgGrA09baNQkbmYiIVCqupsbW2jeBNxM0FhERiYIqMUVEAkoBXEQkoBTARUQCSgFcRCSgFMBFRAJKAVxEJKDiSiP0o/zlhYyes4HtRcWcWbcmw3tlk5dz3BYtIiKBl1YBPH95ISNmraK45BAAhUXFjJi1CkBBXETSTlotoYyes+FI8C5VXHKI0XM2eDQiEZHkSasAvr2oOKrjIiJBllYB/My6NaM6LiISZGkVwIf3yqZmVrVjjtXMqsbwXtkejUhEJHnS6iZm6Y1KZaGISCZIqwAOLogrYItIJkirJRQRkUyiAC4iElAK4CIiAaUALiISUArgIiIBZay1qTuZMTuBb4BdKTupPzUgsz+DTH//oM8g098/RPcZnGWtbVj+YEoDOIAxZom1NjelJ/WZTP8MMv39gz6DTH//kJjPQEsoIiIBpQAuIhJQXgTwSR6c028y/TPI9PcP+gwy/f1DAj6DlK+Bi4hIYmgJRUQkoFIawI0xlxtjNhhjNhpj7k3lub1mjGlijFlgjFlrjFljjLnd6zF5xRhTzRiz3BjzutdjSTVjTF1jzAxjzHpjzDpjzAVejynVjDF3hn4GVhtjXjLG1PB6TMlkjHnaGLPDGLP6qGOnGmP+ZYz5OPRnvVheO2UB3BhTDZgIXAG0BvobY1qn6vw+cBAYZq1tDXQBfpVh7/9otwPrvB6ERx4B/tda2xI4jwz7HIwxjYBfA7nW2rZANeB6b0eVdM8Cl5c7di8wz1rbApgXehy1VM7AOwMbrbWbrLXfAdOAq1J4fk9Zaz+z1i4Lfb0P94ObcfveGmMaA72ByV6PJdWMMacA3YGnAKy131lrizwdlDeqAzWNMdWBWsB2j8eTVNbahcCX5Q5fBUwJfT0FyIvltVMZwBsB2456XEAGBjAAY0wzIAf4wOOheGECcA9w2ONxeKE5sBN4JrSENNkYU9vrQaWStbYQGANsBT4D9lhr/5+3o/LE6dbaz0Jffw6cHsuL6CZmihlj6gAzgTustXu9Hk8qGWP6ADustUu9HotHqgMdgcettTm4bSUy7V5QPdzsszlwJlDbGDPA21F5y7pUwJjSAVMZwAuBJkc9bhw6ljGMMVm44D3VWjvL6/F4oCtwpTFmM24J7VJjzAveDimlCoACa23pldcMXEDPJD8EPrXW7rTWlgCzgAs9HpMXvjDGnAEQ+nNHLC+SygC+GGhhjGlujPke7sbF7BSe31PGGINb+1xnrR3n9Xi8YK0dYa1tbK1thvv3n2+tzZjZl7X2c2CbMaa0y3ZPYK2HQ/LCVqCLMaZW6GeiJxl2IzdkNjA49PVg4LVYXiRlPTGttQeNMbcBc3B3np+21q5J1fl9oCswEFhljFkROnaftfZN74YkHvgfYGpoErMJ+JnH40kpa+0HxpgZwDJcZtZy0rwq0xjzEnAJ0MAYUwCMBEYB040xNwJbgJ/G9NqqxBQRCSbdxBQRCSgFcBGRgFIAFxEJKAVwEZGAUgAXEQkoBXARkYBSABcRCSgFcBGRgPr/CXhh3arQuCIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import random\n", "\n", "n_epoch = 500 # epoch size\n", "a, b = 1, 1 # initial parameters\n", "epsilon = 0.001 # learning rate\n", "\n", "for i in range(n_epoch):\n", " data_idx = list(range(N))\n", " random.shuffle(data_idx)\n", " \n", " for j in data_idx:\n", " a = a + epsilon*2*(Y[j] - a*X[j] - b)*X[j]\n", " b = b + epsilon*2*(Y[j] - a*X[j] - b)\n", "\n", " L = 0\n", " for j in range(N):\n", " L = L + (Y[j]-a*X[j]-b)**2\n", " \n", " if i % 100 == 0:\n", " print(\"epoch %4d: loss = %f, a = %f, b = %f\" % (i, L, a, b))\n", " \n", "x_min = np.min(X)\n", "x_max = np.max(X)\n", "y_min = a * x_min + b\n", "y_max = a * x_max + b\n", "\n", "plt.scatter(X, Y, label='original data')\n", "plt.plot([x_min, x_max], [y_min, y_max], 'r', label='model')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. 如何可视化迭代过程" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "/* global mpl */\n", "window.mpl = {};\n", "\n", "mpl.get_websocket_type = function () {\n", " if (typeof WebSocket !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof MozWebSocket !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert(\n", " 'Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.'\n", " );\n", " }\n", "};\n", "\n", "mpl.figure = function (figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = this.ws.binaryType !== undefined;\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById('mpl-warnings');\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent =\n", " 'This browser does not support binary websocket messages. ' +\n", " 'Performance may be slow.';\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = document.createElement('div');\n", " this.root.setAttribute('style', 'display: inline-block');\n", " this._root_extra_style(this.root);\n", "\n", " parent_element.appendChild(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message('supports_binary', { value: fig.supports_binary });\n", " fig.send_message('send_image_mode', {});\n", " if (fig.ratio !== 1) {\n", " fig.send_message('set_dpi_ratio', { dpi_ratio: fig.ratio });\n", " }\n", " fig.send_message('refresh', {});\n", " };\n", "\n", " this.imageObj.onload = function () {\n", " if (fig.image_mode === 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function () {\n", " fig.ws.close();\n", " };\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "};\n", "\n", "mpl.figure.prototype._init_header = function () {\n", " var titlebar = document.createElement('div');\n", " titlebar.classList =\n", " 'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n", " var titletext = document.createElement('div');\n", " titletext.classList = 'ui-dialog-title';\n", " titletext.setAttribute(\n", " 'style',\n", " 'width: 100%; text-align: center; padding: 3px;'\n", " );\n", " titlebar.appendChild(titletext);\n", " this.root.appendChild(titlebar);\n", " this.header = titletext;\n", "};\n", "\n", "mpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n", "\n", "mpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n", "\n", "mpl.figure.prototype._init_canvas = function () {\n", " var fig = this;\n", "\n", " var canvas_div = (this.canvas_div = document.createElement('div'));\n", " canvas_div.setAttribute(\n", " 'style',\n", " 'border: 1px solid #ddd;' +\n", " 'box-sizing: content-box;' +\n", " 'clear: both;' +\n", " 'min-height: 1px;' +\n", " 'min-width: 1px;' +\n", " 'outline: 0;' +\n", " 'overflow: hidden;' +\n", " 'position: relative;' +\n", " 'resize: both;'\n", " );\n", "\n", " function on_keyboard_event_closure(name) {\n", " return function (event) {\n", " return fig.key_event(event, name);\n", " };\n", " }\n", "\n", " canvas_div.addEventListener(\n", " 'keydown',\n", " on_keyboard_event_closure('key_press')\n", " );\n", " canvas_div.addEventListener(\n", " 'keyup',\n", " on_keyboard_event_closure('key_release')\n", " );\n", "\n", " this._canvas_extra_style(canvas_div);\n", " this.root.appendChild(canvas_div);\n", "\n", " var canvas = (this.canvas = document.createElement('canvas'));\n", " canvas.classList.add('mpl-canvas');\n", " canvas.setAttribute('style', 'box-sizing: content-box;');\n", "\n", " this.context = canvas.getContext('2d');\n", "\n", " var backingStore =\n", " this.context.backingStorePixelRatio ||\n", " this.context.webkitBackingStorePixelRatio ||\n", " this.context.mozBackingStorePixelRatio ||\n", " this.context.msBackingStorePixelRatio ||\n", " this.context.oBackingStorePixelRatio ||\n", " this.context.backingStorePixelRatio ||\n", " 1;\n", "\n", " this.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband_canvas = (this.rubberband_canvas = document.createElement(\n", " 'canvas'\n", " ));\n", " rubberband_canvas.setAttribute(\n", " 'style',\n", " 'box-sizing: content-box; position: absolute; left: 0; top: 0; z-index: 1;'\n", " );\n", "\n", " // Apply a ponyfill if ResizeObserver is not implemented by browser.\n", " if (this.ResizeObserver === undefined) {\n", " if (window.ResizeObserver !== undefined) {\n", " this.ResizeObserver = window.ResizeObserver;\n", " } else {\n", " var obs = _JSXTOOLS_RESIZE_OBSERVER({});\n", " this.ResizeObserver = obs.ResizeObserver;\n", " }\n", " }\n", "\n", " this.resizeObserverInstance = new this.ResizeObserver(function (entries) {\n", " var nentries = entries.length;\n", " for (var i = 0; i < nentries; i++) {\n", " var entry = entries[i];\n", " var width, height;\n", " if (entry.contentBoxSize) {\n", " if (entry.contentBoxSize instanceof Array) {\n", " // Chrome 84 implements new version of spec.\n", " width = entry.contentBoxSize[0].inlineSize;\n", " height = entry.contentBoxSize[0].blockSize;\n", " } else {\n", " // Firefox implements old version of spec.\n", " width = entry.contentBoxSize.inlineSize;\n", " height = entry.contentBoxSize.blockSize;\n", " }\n", " } else {\n", " // Chrome <84 implements even older version of spec.\n", " width = entry.contentRect.width;\n", " height = entry.contentRect.height;\n", " }\n", "\n", " // Keep the size of the canvas and rubber band canvas in sync with\n", " // the canvas container.\n", " if (entry.devicePixelContentBoxSize) {\n", " // Chrome 84 implements new version of spec.\n", " canvas.setAttribute(\n", " 'width',\n", " entry.devicePixelContentBoxSize[0].inlineSize\n", " );\n", " canvas.setAttribute(\n", " 'height',\n", " entry.devicePixelContentBoxSize[0].blockSize\n", " );\n", " } else {\n", " canvas.setAttribute('width', width * fig.ratio);\n", " canvas.setAttribute('height', height * fig.ratio);\n", " }\n", " canvas.setAttribute(\n", " 'style',\n", " 'width: ' + width + 'px; height: ' + height + 'px;'\n", " );\n", "\n", " rubberband_canvas.setAttribute('width', width);\n", " rubberband_canvas.setAttribute('height', height);\n", "\n", " // And update the size in Python. We ignore the initial 0/0 size\n", " // that occurs as the element is placed into the DOM, which should\n", " // otherwise not happen due to the minimum size styling.\n", " if (fig.ws.readyState == 1 && width != 0 && height != 0) {\n", " fig.request_resize(width, height);\n", " }\n", " }\n", " });\n", " this.resizeObserverInstance.observe(canvas_div);\n", "\n", " function on_mouse_event_closure(name) {\n", " return function (event) {\n", " return fig.mouse_event(event, name);\n", " };\n", " }\n", "\n", " rubberband_canvas.addEventListener(\n", " 'mousedown',\n", " on_mouse_event_closure('button_press')\n", " );\n", " rubberband_canvas.addEventListener(\n", " 'mouseup',\n", " on_mouse_event_closure('button_release')\n", " );\n", " rubberband_canvas.addEventListener(\n", " 'dblclick',\n", " on_mouse_event_closure('dblclick')\n", " );\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband_canvas.addEventListener(\n", " 'mousemove',\n", " on_mouse_event_closure('motion_notify')\n", " );\n", "\n", " rubberband_canvas.addEventListener(\n", " 'mouseenter',\n", " on_mouse_event_closure('figure_enter')\n", " );\n", " rubberband_canvas.addEventListener(\n", " 'mouseleave',\n", " on_mouse_event_closure('figure_leave')\n", " );\n", "\n", " canvas_div.addEventListener('wheel', function (event) {\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " on_mouse_event_closure('scroll')(event);\n", " });\n", "\n", " canvas_div.appendChild(canvas);\n", " canvas_div.appendChild(rubberband_canvas);\n", "\n", " this.rubberband_context = rubberband_canvas.getContext('2d');\n", " this.rubberband_context.strokeStyle = '#000000';\n", "\n", " this._resize_canvas = function (width, height, forward) {\n", " if (forward) {\n", " canvas_div.style.width = width + 'px';\n", " canvas_div.style.height = height + 'px';\n", " }\n", " };\n", "\n", " // Disable right mouse context menu.\n", " this.rubberband_canvas.addEventListener('contextmenu', function (_e) {\n", " event.preventDefault();\n", " return false;\n", " });\n", "\n", " function set_focus() {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "};\n", "\n", "mpl.figure.prototype._init_toolbar = function () {\n", " var fig = this;\n", "\n", " var toolbar = document.createElement('div');\n", " toolbar.classList = 'mpl-toolbar';\n", " this.root.appendChild(toolbar);\n", "\n", " function on_click_closure(name) {\n", " return function (_event) {\n", " return fig.toolbar_button_onclick(name);\n", " };\n", " }\n", "\n", " function on_mouseover_closure(tooltip) {\n", " return function (event) {\n", " if (!event.currentTarget.disabled) {\n", " return fig.toolbar_button_onmouseover(tooltip);\n", " }\n", " };\n", " }\n", "\n", " fig.buttons = {};\n", " var buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'mpl-button-group';\n", " for (var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " /* Instead of a spacer, we start a new button group. */\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", " buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'mpl-button-group';\n", " continue;\n", " }\n", "\n", " var button = (fig.buttons[name] = document.createElement('button'));\n", " button.classList = 'mpl-widget';\n", " button.setAttribute('role', 'button');\n", " button.setAttribute('aria-disabled', 'false');\n", " button.addEventListener('click', on_click_closure(method_name));\n", " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", "\n", " var icon_img = document.createElement('img');\n", " icon_img.src = '_images/' + image + '.png';\n", " icon_img.srcset = '_images/' + image + '_large.png 2x';\n", " icon_img.alt = tooltip;\n", " button.appendChild(icon_img);\n", "\n", " buttonGroup.appendChild(button);\n", " }\n", "\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", "\n", " var fmt_picker = document.createElement('select');\n", " fmt_picker.classList = 'mpl-widget';\n", " toolbar.appendChild(fmt_picker);\n", " this.format_dropdown = fmt_picker;\n", "\n", " for (var ind in mpl.extensions) {\n", " var fmt = mpl.extensions[ind];\n", " var option = document.createElement('option');\n", " option.selected = fmt === mpl.default_extension;\n", " option.innerHTML = fmt;\n", " fmt_picker.appendChild(option);\n", " }\n", "\n", " var status_bar = document.createElement('span');\n", " status_bar.classList = 'mpl-message';\n", " toolbar.appendChild(status_bar);\n", " this.message = status_bar;\n", "};\n", "\n", "mpl.figure.prototype.request_resize = function (x_pixels, y_pixels) {\n", " // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n", " // which will in turn request a refresh of the image.\n", " this.send_message('resize', { width: x_pixels, height: y_pixels });\n", "};\n", "\n", "mpl.figure.prototype.send_message = function (type, properties) {\n", " properties['type'] = type;\n", " properties['figure_id'] = this.id;\n", " this.ws.send(JSON.stringify(properties));\n", "};\n", "\n", "mpl.figure.prototype.send_draw_message = function () {\n", " if (!this.waiting) {\n", " this.waiting = true;\n", " this.ws.send(JSON.stringify({ type: 'draw', figure_id: this.id }));\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", " var format_dropdown = fig.format_dropdown;\n", " var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n", " fig.ondownload(fig, format);\n", "};\n", "\n", "mpl.figure.prototype.handle_resize = function (fig, msg) {\n", " var size = msg['size'];\n", " if (size[0] !== fig.canvas.width || size[1] !== fig.canvas.height) {\n", " fig._resize_canvas(size[0], size[1], msg['forward']);\n", " fig.send_message('refresh', {});\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_rubberband = function (fig, msg) {\n", " var x0 = msg['x0'] / fig.ratio;\n", " var y0 = (fig.canvas.height - msg['y0']) / fig.ratio;\n", " var x1 = msg['x1'] / fig.ratio;\n", " var y1 = (fig.canvas.height - msg['y1']) / fig.ratio;\n", " x0 = Math.floor(x0) + 0.5;\n", " y0 = Math.floor(y0) + 0.5;\n", " x1 = Math.floor(x1) + 0.5;\n", " y1 = Math.floor(y1) + 0.5;\n", " var min_x = Math.min(x0, x1);\n", " var min_y = Math.min(y0, y1);\n", " var width = Math.abs(x1 - x0);\n", " var height = Math.abs(y1 - y0);\n", "\n", " fig.rubberband_context.clearRect(\n", " 0,\n", " 0,\n", " fig.canvas.width / fig.ratio,\n", " fig.canvas.height / fig.ratio\n", " );\n", "\n", " fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n", "};\n", "\n", "mpl.figure.prototype.handle_figure_label = function (fig, msg) {\n", " // Updates the figure title.\n", " fig.header.textContent = msg['label'];\n", "};\n", "\n", "mpl.figure.prototype.handle_cursor = function (fig, msg) {\n", " var cursor = msg['cursor'];\n", " switch (cursor) {\n", " case 0:\n", " cursor = 'pointer';\n", " break;\n", " case 1:\n", " cursor = 'default';\n", " break;\n", " case 2:\n", " cursor = 'crosshair';\n", " break;\n", " case 3:\n", " cursor = 'move';\n", " break;\n", " }\n", " fig.rubberband_canvas.style.cursor = cursor;\n", "};\n", "\n", "mpl.figure.prototype.handle_message = function (fig, msg) {\n", " fig.message.textContent = msg['message'];\n", "};\n", "\n", "mpl.figure.prototype.handle_draw = function (fig, _msg) {\n", " // Request the server to send over a new figure.\n", " fig.send_draw_message();\n", "};\n", "\n", "mpl.figure.prototype.handle_image_mode = function (fig, msg) {\n", " fig.image_mode = msg['mode'];\n", "};\n", "\n", "mpl.figure.prototype.handle_history_buttons = function (fig, msg) {\n", " for (var key in msg) {\n", " if (!(key in fig.buttons)) {\n", " continue;\n", " }\n", " fig.buttons[key].disabled = !msg[key];\n", " fig.buttons[key].setAttribute('aria-disabled', !msg[key]);\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_navigate_mode = function (fig, msg) {\n", " if (msg['mode'] === 'PAN') {\n", " fig.buttons['Pan'].classList.add('active');\n", " fig.buttons['Zoom'].classList.remove('active');\n", " } else if (msg['mode'] === 'ZOOM') {\n", " fig.buttons['Pan'].classList.remove('active');\n", " fig.buttons['Zoom'].classList.add('active');\n", " } else {\n", " fig.buttons['Pan'].classList.remove('active');\n", " fig.buttons['Zoom'].classList.remove('active');\n", " }\n", "};\n", "\n", "mpl.figure.prototype.updated_canvas_event = function () {\n", " // Called whenever the canvas gets updated.\n", " this.send_message('ack', {});\n", "};\n", "\n", "// A function to construct a web socket function for onmessage handling.\n", "// Called in the figure constructor.\n", "mpl.figure.prototype._make_on_message_function = function (fig) {\n", " return function socket_on_message(evt) {\n", " if (evt.data instanceof Blob) {\n", " var img = evt.data;\n", " if (img.type !== 'image/png') {\n", " /* FIXME: We get \"Resource interpreted as Image but\n", " * transferred with MIME type text/plain:\" errors on\n", " * Chrome. But how to set the MIME type? It doesn't seem\n", " * to be part of the websocket stream */\n", " img.type = 'image/png';\n", " }\n", "\n", " /* Free the memory for the previous frames */\n", " if (fig.imageObj.src) {\n", " (window.URL || window.webkitURL).revokeObjectURL(\n", " fig.imageObj.src\n", " );\n", " }\n", "\n", " fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n", " img\n", " );\n", " fig.updated_canvas_event();\n", " fig.waiting = false;\n", " return;\n", " } else if (\n", " typeof evt.data === 'string' &&\n", " evt.data.slice(0, 21) === 'data:image/png;base64'\n", " ) {\n", " fig.imageObj.src = evt.data;\n", " fig.updated_canvas_event();\n", " fig.waiting = false;\n", " return;\n", " }\n", "\n", " var msg = JSON.parse(evt.data);\n", " var msg_type = msg['type'];\n", "\n", " // Call the \"handle_{type}\" callback, which takes\n", " // the figure and JSON message as its only arguments.\n", " try {\n", " var callback = fig['handle_' + msg_type];\n", " } catch (e) {\n", " console.log(\n", " \"No handler for the '\" + msg_type + \"' message type: \",\n", " msg\n", " );\n", " return;\n", " }\n", "\n", " if (callback) {\n", " try {\n", " // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n", " callback(fig, msg);\n", " } catch (e) {\n", " console.log(\n", " \"Exception inside the 'handler_\" + msg_type + \"' callback:\",\n", " e,\n", " e.stack,\n", " msg\n", " );\n", " }\n", " }\n", " };\n", "};\n", "\n", "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n", "mpl.findpos = function (e) {\n", " //this section is from http://www.quirksmode.org/js/events_properties.html\n", " var targ;\n", " if (!e) {\n", " e = window.event;\n", " }\n", " if (e.target) {\n", " targ = e.target;\n", " } else if (e.srcElement) {\n", " targ = e.srcElement;\n", " }\n", " if (targ.nodeType === 3) {\n", " // defeat Safari bug\n", " targ = targ.parentNode;\n", " }\n", "\n", " // pageX,Y are the mouse positions relative to the document\n", " var boundingRect = targ.getBoundingClientRect();\n", " var x = e.pageX - (boundingRect.left + document.body.scrollLeft);\n", " var y = e.pageY - (boundingRect.top + document.body.scrollTop);\n", "\n", " return { x: x, y: y };\n", "};\n", "\n", "/*\n", " * return a copy of an object with only non-object keys\n", " * we need this to avoid circular references\n", " * http://stackoverflow.com/a/24161582/3208463\n", " */\n", "function simpleKeys(original) {\n", " return Object.keys(original).reduce(function (obj, key) {\n", " if (typeof original[key] !== 'object') {\n", " obj[key] = original[key];\n", " }\n", " return obj;\n", " }, {});\n", "}\n", "\n", "mpl.figure.prototype.mouse_event = function (event, name) {\n", " var canvas_pos = mpl.findpos(event);\n", "\n", " if (name === 'button_press') {\n", " this.canvas.focus();\n", " this.canvas_div.focus();\n", " }\n", "\n", " var x = canvas_pos.x * this.ratio;\n", " var y = canvas_pos.y * this.ratio;\n", "\n", " this.send_message(name, {\n", " x: x,\n", " y: y,\n", " button: event.button,\n", " step: event.step,\n", " guiEvent: simpleKeys(event),\n", " });\n", "\n", " /* This prevents the web browser from automatically changing to\n", " * the text insertion cursor when the button is pressed. We want\n", " * to control all of the cursor setting manually through the\n", " * 'cursor' event from matplotlib */\n", " event.preventDefault();\n", " return false;\n", "};\n", "\n", "mpl.figure.prototype._key_event_extra = function (_event, _name) {\n", " // Handle any extra behaviour associated with a key event\n", "};\n", "\n", "mpl.figure.prototype.key_event = function (event, name) {\n", " // Prevent repeat events\n", " if (name === 'key_press') {\n", " if (event.key === this._key) {\n", " return;\n", " } else {\n", " this._key = event.key;\n", " }\n", " }\n", " if (name === 'key_release') {\n", " this._key = null;\n", " }\n", "\n", " var value = '';\n", " if (event.ctrlKey && event.key !== 'Control') {\n", " value += 'ctrl+';\n", " }\n", " else if (event.altKey && event.key !== 'Alt') {\n", " value += 'alt+';\n", " }\n", " else if (event.shiftKey && event.key !== 'Shift') {\n", " value += 'shift+';\n", " }\n", "\n", " value += 'k' + event.key;\n", "\n", " this._key_event_extra(event, name);\n", "\n", " this.send_message(name, { key: value, guiEvent: simpleKeys(event) });\n", " return false;\n", "};\n", "\n", "mpl.figure.prototype.toolbar_button_onclick = function (name) {\n", " if (name === 'download') {\n", " this.handle_save(this, null);\n", " } else {\n", " this.send_message('toolbar_button', { name: name });\n", " }\n", "};\n", "\n", "mpl.figure.prototype.toolbar_button_onmouseover = function (tooltip) {\n", " this.message.textContent = tooltip;\n", "};\n", "\n", "///////////////// REMAINING CONTENT GENERATED BY embed_js.py /////////////////\n", "// prettier-ignore\n", "var _JSXTOOLS_RESIZE_OBSERVER=function(A){var t,i=new WeakMap,n=new WeakMap,a=new WeakMap,r=new WeakMap,o=new Set;function s(e){if(!(this instanceof s))throw new TypeError(\"Constructor requires 'new' operator\");i.set(this,e)}function h(){throw new TypeError(\"Function is not a constructor\")}function c(e,t,i,n){e=0 in arguments?Number(arguments[0]):0,t=1 in arguments?Number(arguments[1]):0,i=2 in arguments?Number(arguments[2]):0,n=3 in arguments?Number(arguments[3]):0,this.right=(this.x=this.left=e)+(this.width=i),this.bottom=(this.y=this.top=t)+(this.height=n),Object.freeze(this)}function d(){t=requestAnimationFrame(d);var s=new WeakMap,p=new Set;o.forEach((function(t){r.get(t).forEach((function(i){var r=t instanceof window.SVGElement,o=a.get(t),d=r?0:parseFloat(o.paddingTop),f=r?0:parseFloat(o.paddingRight),l=r?0:parseFloat(o.paddingBottom),u=r?0:parseFloat(o.paddingLeft),g=r?0:parseFloat(o.borderTopWidth),m=r?0:parseFloat(o.borderRightWidth),w=r?0:parseFloat(o.borderBottomWidth),b=u+f,F=d+l,v=(r?0:parseFloat(o.borderLeftWidth))+m,W=g+w,y=r?0:t.offsetHeight-W-t.clientHeight,E=r?0:t.offsetWidth-v-t.clientWidth,R=b+v,z=F+W,M=r?t.width:parseFloat(o.width)-R-E,O=r?t.height:parseFloat(o.height)-z-y;if(n.has(t)){var k=n.get(t);if(k[0]===M&&k[1]===O)return}n.set(t,[M,O]);var S=Object.create(h.prototype);S.target=t,S.contentRect=new c(u,d,M,O),s.has(i)||(s.set(i,[]),p.add(i)),s.get(i).push(S)}))})),p.forEach((function(e){i.get(e).call(e,s.get(e),e)}))}return s.prototype.observe=function(i){if(i instanceof window.Element){r.has(i)||(r.set(i,new Set),o.add(i),a.set(i,window.getComputedStyle(i)));var n=r.get(i);n.has(this)||n.add(this),cancelAnimationFrame(t),t=requestAnimationFrame(d)}},s.prototype.unobserve=function(i){if(i instanceof window.Element&&r.has(i)){var n=r.get(i);n.has(this)&&(n.delete(this),n.size||(r.delete(i),o.delete(i))),n.size||r.delete(i),o.size||cancelAnimationFrame(t)}},A.DOMRectReadOnly=c,A.ResizeObserver=s,A.ResizeObserverEntry=h,A}; // eslint-disable-line\n", "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Left button pans, Right button zooms\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n", "\n", "mpl.extensions = [\"eps\", \"jpeg\", \"pgf\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n", "\n", "mpl.default_extension = \"png\";/* global mpl */\n", "\n", "var comm_websocket_adapter = function (comm) {\n", " // Create a \"websocket\"-like object which calls the given IPython comm\n", " // object with the appropriate methods. Currently this is a non binary\n", " // socket, so there is still some room for performance tuning.\n", " var ws = {};\n", "\n", " ws.binaryType = comm.kernel.ws.binaryType;\n", " ws.readyState = comm.kernel.ws.readyState;\n", " function updateReadyState(_event) {\n", " if (comm.kernel.ws) {\n", " ws.readyState = comm.kernel.ws.readyState;\n", " } else {\n", " ws.readyState = 3; // Closed state.\n", " }\n", " }\n", " comm.kernel.ws.addEventListener('open', updateReadyState);\n", " comm.kernel.ws.addEventListener('close', updateReadyState);\n", " comm.kernel.ws.addEventListener('error', updateReadyState);\n", "\n", " ws.close = function () {\n", " comm.close();\n", " };\n", " ws.send = function (m) {\n", " //console.log('sending', m);\n", " comm.send(m);\n", " };\n", " // Register the callback with on_msg.\n", " comm.on_msg(function (msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", " var data = msg['content']['data'];\n", " if (data['blob'] !== undefined) {\n", " data = {\n", " data: new Blob(msg['buffers'], { type: data['blob'] }),\n", " };\n", " }\n", " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(data);\n", " });\n", " return ws;\n", "};\n", "\n", "mpl.mpl_figure_comm = function (comm, msg) {\n", " // This is the function which gets called when the mpl process\n", " // starts-up an IPython Comm through the \"matplotlib\" channel.\n", "\n", " var id = msg.content.data.id;\n", " // Get hold of the div created by the display call when the Comm\n", " // socket was opened in Python.\n", " var element = document.getElementById(id);\n", " var ws_proxy = comm_websocket_adapter(comm);\n", "\n", " function ondownload(figure, _format) {\n", " window.open(figure.canvas.toDataURL());\n", " }\n", "\n", " var fig = new mpl.figure(id, ws_proxy, ondownload, element);\n", "\n", " // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n", " // web socket which is closed, not our websocket->open comm proxy.\n", " ws_proxy.onopen();\n", "\n", " fig.parent_element = element;\n", " fig.cell_info = mpl.find_output_cell(\"
\");\n", " if (!fig.cell_info) {\n", " console.error('Failed to find cell for figure', id, fig);\n", " return;\n", " }\n", " fig.cell_info[0].output_area.element.on(\n", " 'cleared',\n", " { fig: fig },\n", " fig._remove_fig_handler\n", " );\n", "};\n", "\n", "mpl.figure.prototype.handle_close = function (fig, msg) {\n", " var width = fig.canvas.width / fig.ratio;\n", " fig.cell_info[0].output_area.element.off(\n", " 'cleared',\n", " fig._remove_fig_handler\n", " );\n", " fig.resizeObserverInstance.unobserve(fig.canvas_div);\n", "\n", " // Update the output cell to use the data from the current canvas.\n", " fig.push_to_output();\n", " var dataURL = fig.canvas.toDataURL();\n", " // Re-enable the keyboard manager in IPython - without this line, in FF,\n", " // the notebook keyboard shortcuts fail.\n", " IPython.keyboard_manager.enable();\n", " fig.parent_element.innerHTML =\n", " '';\n", " fig.close_ws(fig, msg);\n", "};\n", "\n", "mpl.figure.prototype.close_ws = function (fig, msg) {\n", " fig.send_message('closing', msg);\n", " // fig.ws.close()\n", "};\n", "\n", "mpl.figure.prototype.push_to_output = function (_remove_interactive) {\n", " // Turn the data on the canvas into data in the output cell.\n", " var width = this.canvas.width / this.ratio;\n", " var dataURL = this.canvas.toDataURL();\n", " this.cell_info[1]['text/html'] =\n", " '';\n", "};\n", "\n", "mpl.figure.prototype.updated_canvas_event = function () {\n", " // Tell IPython that the notebook contents must change.\n", " IPython.notebook.set_dirty(true);\n", " this.send_message('ack', {});\n", " var fig = this;\n", " // Wait a second, then push the new image to the DOM so\n", " // that it is saved nicely (might be nice to debounce this).\n", " setTimeout(function () {\n", " fig.push_to_output();\n", " }, 1000);\n", "};\n", "\n", "mpl.figure.prototype._init_toolbar = function () {\n", " var fig = this;\n", "\n", " var toolbar = document.createElement('div');\n", " toolbar.classList = 'btn-toolbar';\n", " this.root.appendChild(toolbar);\n", "\n", " function on_click_closure(name) {\n", " return function (_event) {\n", " return fig.toolbar_button_onclick(name);\n", " };\n", " }\n", "\n", " function on_mouseover_closure(tooltip) {\n", " return function (event) {\n", " if (!event.currentTarget.disabled) {\n", " return fig.toolbar_button_onmouseover(tooltip);\n", " }\n", " };\n", " }\n", "\n", " fig.buttons = {};\n", " var buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'btn-group';\n", " var button;\n", " for (var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " /* Instead of a spacer, we start a new button group. */\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", " buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'btn-group';\n", " continue;\n", " }\n", "\n", " button = fig.buttons[name] = document.createElement('button');\n", " button.classList = 'btn btn-default';\n", " button.href = '#';\n", " button.title = name;\n", " button.innerHTML = '';\n", " button.addEventListener('click', on_click_closure(method_name));\n", " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", " buttonGroup.appendChild(button);\n", " }\n", "\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = document.createElement('span');\n", " status_bar.classList = 'mpl-message pull-right';\n", " toolbar.appendChild(status_bar);\n", " this.message = status_bar;\n", "\n", " // Add the close button to the window.\n", " var buttongrp = document.createElement('div');\n", " buttongrp.classList = 'btn-group inline pull-right';\n", " button = document.createElement('button');\n", " button.classList = 'btn btn-mini btn-primary';\n", " button.href = '#';\n", " button.title = 'Stop Interaction';\n", " button.innerHTML = '';\n", " button.addEventListener('click', function (_evt) {\n", " fig.handle_close(fig, {});\n", " });\n", " button.addEventListener(\n", " 'mouseover',\n", " on_mouseover_closure('Stop Interaction')\n", " );\n", " buttongrp.appendChild(button);\n", " var titlebar = this.root.querySelector('.ui-dialog-titlebar');\n", " titlebar.insertBefore(buttongrp, titlebar.firstChild);\n", "};\n", "\n", "mpl.figure.prototype._remove_fig_handler = function (event) {\n", " var fig = event.data.fig;\n", " if (event.target !== this) {\n", " // Ignore bubbled events from children.\n", " return;\n", " }\n", " fig.close_ws(fig, {});\n", "};\n", "\n", "mpl.figure.prototype._root_extra_style = function (el) {\n", " el.style.boxSizing = 'content-box'; // override notebook setting of border-box.\n", "};\n", "\n", "mpl.figure.prototype._canvas_extra_style = function (el) {\n", " // this is important to make the div 'focusable\n", " el.setAttribute('tabindex', 0);\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " } else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "};\n", "\n", "mpl.figure.prototype._key_event_extra = function (event, _name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager) {\n", " manager = IPython.keyboard_manager;\n", " }\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which === 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", " fig.ondownload(fig, null);\n", "};\n", "\n", "mpl.find_output_cell = function (html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i = 0; i < ncells; i++) {\n", " var cell = cells[i];\n", " if (cell.cell_type === 'code') {\n", " for (var j = 0; j < cell.output_area.outputs.length; j++) {\n", " var data = cell.output_area.outputs[j];\n", " if (data.data) {\n", " // IPython >= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] === html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "};\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel !== null) {\n", " IPython.notebook.kernel.comm_manager.register_target(\n", " 'matplotlib',\n", " mpl.mpl_figure_comm\n", " );\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib nbagg\n", "\n", "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "\n", "n_epoch = 300 # epoch size\n", "a, b = 1, 1 # initial parameters\n", "epsilon = 0.0001 # learning rate\n", "\n", "fig = plt.figure()\n", "imgs = []\n", "\n", "for i in range(n_epoch):\n", " data_idx = list(range(N))\n", " random.shuffle(data_idx)\n", " \n", " for j in data_idx[:10]:\n", " a = a + epsilon*2*(Y[j] - a*X[j] - b)*X[j]\n", " b = b + epsilon*2*(Y[j] - a*X[j] - b)\n", "\n", "\n", " if i<80 and i % 5 == 0:\n", " x_min = np.min(X)\n", " x_max = np.max(X)\n", " y_min = a * x_min + b\n", " y_max = a * x_max + b\n", "\n", " img = plt.scatter(X, Y, label='original data')\n", " img = plt.plot([x_min, x_max], [y_min, y_max], 'r', label='model')\n", " imgs.append(img)\n", " \n", "ani = animation.ArtistAnimation(fig, imgs)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. 如何使用批次更新的方法?\n", "\n", "如果有一些数据包含比较大的错误(异常数据),因此每次更新仅仅使用一个数据会导致不精确,同时每次仅仅使用一个数据来计算更新也导致计算效率比较低。\n", "\n", "\n", "* [梯度下降方法的几种形式](https://blog.csdn.net/u010402786/article/details/51188876)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. 如何拟合多项式函数?\n", "\n", "需要设计一个弹道导弹防御系统,通过观测导弹的飞行路径,预测未来导弹的飞行轨迹,从而完成摧毁的任务。按照物理学,可以得知模型为:\n", "$$\n", "y = at^2 + bt + c\n", "$$\n", "我们需要求解三个模型参数$a, b, c$。\n", "\n", "损失函数的定义为:\n", "$$\n", "L = \\sum_{i=1}^N (y_i - at_i^2 - bt_i - c)^2\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEGCAYAAACZ0MnKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoPElEQVR4nO3dd3xV9f3H8dfnZpLJSAiQICBEMGwIiIoTVOpCrbUuBKVS66jUts7+rG1/rl+tVmsdCAIqiruOunC0isoIIDMCEQgbwgohkP39/ZGDokK4kNyc3OT9fDzuI/eec+4976vAO2d9jznnEBERCUbA7wAiIhI+VBoiIhI0lYaIiARNpSEiIkFTaYiISNAi/Q4QSikpKa5jx45+xxARCStz5szZ4pxL3d+8Rl0aHTt2JCcnx+8YIiJhxczyDzRPu6dERCRoKg0REQmaSkNERIKm0hARkaCpNEREJGgqDRERCZpKQ0REgqbSkJBYv2MPz87Ip7Si0u8oIlKHGvXFfVL/Kqscz365ir++v5TiskqmLy/gn5f2IzJCv5+INAb6myx1ZunGIi584gvuemsJ/Tu2ZOzQTN5fvInfv7KAqird7EukMdCWRiO2fscePl1WQM+MZI5uk0QgYCFZT0l5JY9+nMcT//2GpGZRPHxxH87t3Q4zIyoiwF/fX0qz6AjuPq8HZqHJICL1Q6XRCFVWOSZ9sYq/fbCU3WXVxxRaxUdzXJcUBndpxeDMVNKbN6vVOrYVl/H1xp18vaGI52bks2JLMRf0S+cPZ2XRMj762+WuO6ULu0orePw/3xAfHcHtZx6t4hAJYyqNRmbRukJue20hC9cVcnLXVG467SiWb9rF53lbmJ63hbfmrwfgyNR4/ve8HhzXOSWoz52Tv50PFm8kd2MRX2/Yyeai0m/ndUqJ55mrBnLiUfsdFJObz+jK7tIKnvpsJfExkYwdelTtv6iI+EKl0UgUl1bw4LRlTPx8Ja0SYnj00r6c1bMtZkavjOb8tH8GzjmWb97FZ8u3MGVmPqOens1DP+/DWb3a1vjZz87I5643FxNhRpfWCQzOTOHoNkl0bZNIt7aJpCbE1Lj1YGb88ZzuFJdV8vcPl5MQE8kvTjiyrv8TiEg9MOca7wHK7Oxs19iGRt9dVkFBUel3j12lbN5Zyuvz1rFuxx4uPeYIbhnWjeRmUTV+TuHuckZPns2c1du565zujDyu44+Wqais4n//ncukL1ZxarfWPHxxHxJja/7cmlRUVvHrqfN4Z+FGWsRFkRgbRWJspPeIIik2ioGdWnBu73SaRUcc9npEpHbMbI5zLnu/81QaDdvmnSVM93YtfZG3lY07S360TMAgq10Sd53TneyOLYP+7JLySq5/fh4f5m7iulM687vTu367xVBUUs4NL8zjP0sLGD24E7efeTQRdXAgvayiiomfr2TN9t0UlVR4j3KKSirYVlzG5qJSkptFcfGA9lw+qAPtW8bVep0icmhUGmHmi7wtTMvdxOd5W1i2aRcALeKiOK5LCt3bJdE6MZbUxBhSE2JITYyhZXz0Yf+DXlFZxf+8sYgXZq3houwM7jm/JxsKS/jF5BzyCnbx5+HdueyYDnX59Q7IOcesldt45st83lu8kSrnGNItjVHHdeS4zq0oq6yiuLSC3WWV7C6rpLisguiIAFltQ3dmmEhTVFNp6JhGA7K5qIQ7/7WY9xZvJCYywMBOLbmgXwaDu6SE7B/GyIgA95zfk9SEGB75OI91O/bw9YYiyiurmHzlQAZnBnegvC6YGccc2YpjjmzFhsI9TJmxmhdmrebD3E2YwYF+v0lJiOHUbqkMOTqNEzJTiIvWH2uRUNGWRgPgnOPVuev4y9tL2FNeydihmVx1fCdio+p3v/6zM/K5841FHNEyjgkjB9CldUK9rn9/SsoreXfRBr7ZXExcTATx0ZHERUcQHxNJs+gIduwu46Pczfx3WQFFJRVERwY4vnMrhhydxmlZaaQlxfr9FUTCjnZP1YOqKkeVc4c8XMa6HXu47bWFfLqsgOwOLbj/wl50TvXvH+ulG4to2zyWpFoc8PZDeWUVs1duY1ruJj7K3czqbbsB6N2+OadnpXF6VhpdWifoGhGRIKg0QmxXaQW/mDyb/K27+dtFvYO69qGyyvH8zHzue/drHHDLsG6MGNRB++brwN5Ti6ct2cQHSzYxf80OADq2iuP07m0Yc+KRpCTE+BtSpAFTaYTQzpJyRj09i/lrC2mbHMu6HXsYc+KR/Pa0rkRH7n+rI2fVNu56azGL1u3khMwU7jm/p84SCqGNhSV8mFtdIF9+s4Wj0hJ56ZfHEh+jYx8i+6PSCJHC3eVc8fRMFq/fyaOX9uXEo1L5y9u5vDBrNT3Sk3j44r7f29W0sbCE+97N5V9fradNUiy3ndnt2zGapH588vVmRk+ezSldWzPuiuw6OY1YpLFRaYTA9uIyRjw9k2Ubd/HYZf0YmpX27bz3F2/k1lcXsKe8kjvP7s4F/dKZMH0l//wkj4oqx5gTjuTaUzrrLB+fPDsjn//51yJGHtuBu87tfsDSLquo4h8fLyd3QxGnZbXm9Kw2tNhnXC2RxkqlUce27irl8gmz+KZgF0+O6M8pXVv/aJlNO0v47UvzmZ63hcTYSIpKKjg9K40/nJXFEa20K8pvd/97CU99tpL/OTuL0YM7/Wh+/tZibnhhHgvWFtImKZaNO0uIDBjHd0nhrJ5tOb17Gs3jVCDSODXY0jCz3wC/ABywELgSaAtMBVoBc4ARzrkyM4sBngH6A1uBnzvnVtX0+aEojYKiUi4bP4P8rbsZPzKbEzL3P0gfVJ9R9fTnK5m2ZBPXn9qlxmWlflVVOa6dMpf3l2zkicv7c0b3Nt/Oe+Orddzx+iICBv93YW/O6J7G4vU7eXvBBt5ZuIHV23YTGTBOy0rjwYv6aMgTaXQaZGmYWTowHchyzu0xs5eAd4Azgdecc1PN7AlgvnPucTO7FujlnLvGzC4GznfO/bymddS2NLbsKmXZxiKWbipi2aYilm6sflQ5mDAqO+gRYqVh2lNWycVPzWDpxp28OOZYMtMSuOvNxbyUs5b+HVrw8MV9yGjx/a1C5xyL1u3kzfnrGD99JRcPOIJ7L+jp0zcQCY2GfEV4JNDMzMqBOGADcCpwqTd/MnAX8Dgw3HsO8ArwqJmZC0HrbSjcw9mPTGdrcdm305rHRdE1LZGf9s/gwv4Z9MpoXterlXrWLDqC8Vdkc/5jnzN6cg7JzSJZsaWY60/pwtihmfu95sbM6JmRTM+MZAIB48n/ruCko1IY1qPmkYJFGgvfSsM5t87MHgBWA3uAD6jeHbXDOVfhLbYWSPeepwNrvPdWmFkh1buwttR1ttSEGE7vnkaX1ol0TUvkqDYJBx3+W8JTamIMk64cwAWPfcHOEnhu9DEc3yW4LcjfntaVL7/Zyi2vLqRXRnPa1fLGViLhwLfSMLMWVG89dAJ2AC8Dw+rgc8cAYwCOOOKIw/qMyIgA917Qq7ZRJEx0aZ3IhzedRGx0xCFdCR8dGeCRi/ty1iOfMXbqV7wwZpBO4ZVG79DGvKhbQ4GVzrkC51w58BpwPNDczPaWWQawznu+DmgP4M1PpvqA+Pc458Y557Kdc9mpqTrwLMFpnXR4Q6d0TInnL+f1YNaqbTz6cV4Ikok0LH6WxmpgkJnFWfV+nyHAEuAT4EJvmZHAG97zN73XePM/DsXxDJFDdUG/DM7r046HP1pGzqptfscRCSnfSsM5N5PqA9pzqT7dNgCMA24BbjKzPKqPWUzw3jIBaOVNvwm4td5DixzAX87rQUaLOG6c+hWFe8r9jiMSMrq4T6SOzFu9nZ898SVndG/DPef3JKlZpE6ekLDUkE+5FWk0+h7RgptOP4r/e28p/164gagIo1V8DCmJ0dU/E2KIi44gKiJAVKQRHRGofh4RYHCXFHpmJPv9FUQOSqUhUoeuObEzR7VOJH/bbrbsKmVLUSlbdpWytbiM5ZuKKKmooryiirLK6sfeDf2HIgMHHJJGpCFRaYjUoUDAvjd45cFUVjm27irlykmzGfNMDo9e2u97Q5qINDR+nj0l0uRFBIzWSbE8f/UgurdL5topc3lr/nq/Y4kckEpDpAFIbhbFc784hv5HtODGqfN4Zc5avyOJ7JdKQ6SBSIiJZNJVAziucwq/e3k+U2bm+x1J5EdUGiINSFx0JONHZnNqt9bc8foiJn6+0u9IIt+j0hBpYGKjIrx7fKTxp7eW8KaOcUgDotIQaYCiIwM8cklfBnRswe9ens+cfA1PIg2DSkOkgYqJjGDciGzaJcdy9TNzyN9a7HckEZWGSEPWIj6aiVcOpMo5rpw0m8LdGtdK/KXSEGngOqXEM25ENmu37eGXz+VQVlHldyRpwlQaImFgYKeW3H9hT2as2MZtry2kMQ80Kg2bhhERCRPn980gf+tu/v7hcjq2iuOGIZl+R5ImSKUhEkZuHJJJ/tbd/G3aMrYWl3Hbmd2IiYzwO5Y0ISoNkTBiZtz/0160iIvm6c9XMm/1dh69tB/tW8b5HU2aCB3TEAkz0ZEB7jwniydH9GfFlmLOfOQz3lu00e9Y0kSoNETC1Bnd2/DOr0/gyJR4rnluDn96a7HOrJKQ0+4pkTDWvmUcL19zHPe+m8vEz1cxJ3875/VJp1NqPEemxJPevBmREfrdUOqO7hEu0ki8t2gjf/jXIrbsKv12WlSE0b5lHEemxPOrk7vQv0MLHxNKuKjpHuEqDZFGxDnH1uIyVm4p/u5RUExO/nbA8d7YE0lJiPE7pjRwNZWGdk+JNCJmRkpCDCkJMQzo2PLb6Us3FnHOo9O59dUFPHVFNmbmY0oJZ9rZKdIEdG2TyK3DuvFh7mZemLXG7zgSxlQaIk3EqOM6ckJmCn95ewkrCnb5HUfClK+lYWbNzewVM/vazHLN7Fgza2lm08xsufezhbesmdkjZpZnZgvMrJ+f2UXCTSBgPPCz3sREBfjNi19RXqnTc+XQ+b2l8TDwnnOuG9AbyAVuBT5yzmUCH3mvAX4CZHqPMcDj9R9XJLylJcVyz/k9mb+2kH98tNzvOBKGfCsNM0sGTgQmADjnypxzO4DhwGRvscnAed7z4cAzrtoMoLmZta3X0CKNwJk923Jh/wwe/SRPdwSUQ+bnlkYnoACYaGbzzGy8mcUDac65Dd4yG4E073k6sO8RvLXetO8xszFmlmNmOQUFBSGMLxK+/nhOFuktmjH2xa8oKtGNnSR4fpZGJNAPeNw51xco5rtdUQC46otIDulCEufcOOdctnMuOzU1tc7CijQmibFRPHRRH9Zt38P1z89jQ+EevyNJmPCzNNYCa51zM73Xr1BdIpv27nbyfm725q8D2u/z/gxvmogchuyOLfnTud2ZsWIrpz7wXx79eDkl5ZV+x5IGzrfScM5tBNaYWVdv0hBgCfAmMNKbNhJ4w3v+JnCFdxbVIKBwn91YInIYRhzbkQ9vOomTu6bywAfLOO2h//L+4o26M6AckK/DiJhZH2A8EA2sAK6kusheAo4A8oGLnHPbrPoS1keBYcBu4ErnXI1jhGgYEZHgfZG3hbveWsyyTbs4ITOFP56TRZfWiX7HEh9o7CkRCUpFZRXPzcjnwWnLKK90vHzNsfRIT/Y7ltSzmkrD7+s0RKQBiYwIMOr4Tky76SRaxEVx9TM5bN5Z4ncsaUBUGiLyI2lJsYwfOYDCPeVc/ewcHSCXb6k0RGS/stol8dDP+7Bg7Q5+/8oCHRwXQKUhIjU4o3sbbj6jG2/NX88/Ps7zO440ALqfhojU6JqTjmT55iIenLaMI1PjObtXO78jiY+0pSEiNTIz7r2gJ9kdWvDbl+Yzf80OvyOJj1QaInJQMZERPDGiPykJMVz9TA75W4v9jiQ+UWmISFBSEmJ4etQAyiqruOjJL8nbXOR3JPGBSkNEgta1TSIvjjmWKgcXPTmDxesL/Y4k9UylISKHpGubRF765bHERga4ZNwM5q7e7nckqUcqDRE5ZJ1S4nnpmmNpER/NiPEz+fKbrX5Hknqi0hCRw5LRIo6Xf3ks7Zo3Y9TEWXyydPPB3yRhTwMWikitbCsuY8SEmSzbVMQJmalktk6gS+sEMtMS6dI6gYQYXQ4WbmoasFD/N0WkVlrGR/P81YO4951cvlqzg+nLt1BWWfXt/HbJsdw8rBvn9f3R3ZklDKk0RKTWkptFcd9PewHVw6uv2b6H5ZuKWL55F2/NX88f31zMkKNbkxgb5XNSqS0d0xCROhUZEaBTSjynd2/Ddad04f8u7EXhnnImfb7K72hSB1QaIhJSvTKaM6Rba8ZPX8nOknK/40gtqTREJOTGDj1KWxuNhEpDREKuZ0YyQ49OY/xnK7S1EeZUGiJSL8YOzWRnSQUTp6/yO4rUgkpDROpFj/RkTstKY8L0FRTu0dZGuFJpiEi9uXGIt7Xx+Uq/o8hhUmmISL3pkZ7M6VlpTJi+UlsbYUqlISL16sahmRSVVPD0dG1thCPfS8PMIsxsnpm97b3uZGYzzSzPzF40s2hveoz3Os+b39HX4CJyWLq3S+aM7mk8/bm2NsKR76UB3Ajk7vP6fuAh51wXYDsw2ps+GtjuTX/IW05EwtDYoUdRVFLBA+8v1Sm4YcbX0jCzDOAsYLz32oBTgVe8RSYD53nPh3uv8eYP8ZYXkTBzdNskLuiXzrMz8un/l2lc8fQspszMZ/POEr+jyUEEVRpm9mww0w7D34Gbgb1DYrYCdjjnKrzXa4G9Q2OmA2sAvPmF3vI/zDXGzHLMLKegoKAOIopIKDxwYW9e/dWxXHl8J/K3FnPH64sYeM9HnP/Y50yZme93PDmAYEe57b7vCzOLAPrXZsVmdjaw2Tk3x8xOrs1n7cs5Nw4YB9X306irzxWRuhUIGP07tKR/h5bc9pNuLNu0iw8Wb+SdRRu54/VFtIqPZliPtn7HlB+ocUvDzG4zsyKgl5nt9B5FwGbgjVqu+3jgXDNbBUylerfUw0BzM9tbZhnAOu/5OqC9lysSSAZ0j0mRRsDM6NomkRuGZPLW9ceT1TaJO99YrOMdDVCNpeGcu9c5lwj81TmX5D0SnXOtnHO31WbFzrnbnHMZzrmOwMXAx865y4BPgAu9xUbyXTm96b3Gm/+xa8y3HRRpoiIjAtz3055s2VXK/e9+7Xcc+YGgjmk4524zs3QzO87MTtz7CFGmW4CbzCyP6mMWE7zpE4BW3vSbgFtDtH4R8VmvjOZceXwnpsxczexV2/yOI/sI6h7hZnYf1VsDS4BKb7Jzzp0bwmy1pnuEi4Sv4tIKTn/oU5pFR/DvXw8mJjLC70hNRk33CA/2lNvzga7OuTOdc+d4jwZdGCIS3uJjIvnf83qQt3kXT/xnhd9xxBNsaawAdHNfEalXp3RrzTm92/HPT/LI27zL7zjCQU65NbN/AA7YDXxlZh8BpXvnO+d+Hdp4ItLU3Xl2Fp8uK+D21xYydcwgAgFd0+ung12nsfeAwByqz14SEalXqYkx3HHm0dz86gJezFnDJQOP8DtSk1ZjaTjnJtc0X0SkPvwsO4PX5q3lnndyOSEzhYwWcX5HarKCHUZkoZkt+MHjMzN7yMx+NJSHiEhdMjPuvaAXAFdMmMXWXaUHeYeESrAHwt8F/g1c5j3eonrX1UZgUkiSiYjso1NKPE+PGsD6wj2MmjibIl0t7otgS2OodwX3Qu9xB3CSc+5+oGPo4omIfGdAx5Y8fll/cjfsZMwzcygprzz4m6ROBVsaEWY2cO8LMxsA7L3SpmL/bxERqXundGvNAz/rzZcrtvLrF+ZRUVl18DdJnQm2NH4BTDCzld4AgxOAq80sHrg3VOFERPbnvL7p3HVOFh8s2cRtry1Ew9DVn6CGRnfOzQZ6mlmy97pwn9kvhSKYiEhNRh3fie27y3n4o+U0j4vi9jOPRvdlC72DXdx3uXPuOTO76QfTAXDOPRjCbCIiNRo7NJPtu8t46rOVtEqI4ZqTOvsdqdE72JZGvPczMdRBREQOlZlx1znd2VZcxn3vfk1aUgzn983wO1ajdrCL+570fv6pfuKIiByaQMD420W92bqrjN+/vICUhBhOyEz1O1ajFezFfUeZ2Udmtsh73cvM/hDaaCIiwYmJjODJK/rTpXUC1zw7h0XrCg/+JjkswZ499RRwG1AO4JxbQPX9NUREGoSk2CgmXzWQ5nHRjJo4m9Vbd/sdqVEKtjTinHOzfjBN12eISIOSlhTL5KsGUF5ZxciJs9hWXOZ3pEYn2NLYYmadqR4mHTO7ENgQslQiIoepS+tEJozMZv2OPVw1aTa7y/T7bV0KtjSuA54EupnZOmAscE2oQomI1EZ2x5Y8cklfFqzdwdXP5GicqjoUbGmsAyYCdwNTgWnAyFCFEhGprTO6t+GBn/Vmxopt/PzJGWzeWeJ3pEYh2NJ4AziH6gPh64FdQHGoQomI1IUL+mUwYWQ2q7YWc8HjX/BNgW4ZW1sWzJgtZrbIOdejHvLUqezsbJeTk3PwBUWkUVuwdgdXTZpNRZVjwsgB9O/Qwu9IDZqZzXHOZe9vXrBbGl+YWc86zCQiUm96ZTTn1V8dR/NmUVz61AymLdnkd6SwVWNp7L1jHzAYmGtmS7279u2dftjMrL2ZfWJmS8xssZnd6E1vaWbTzGy597OFN93M7BEzy/My9KvN+kWkaenQKp5XfnUc3dok8stnc3h+5mq/I4Wlg409dXYI110B/NY5N9fMEoE5ZjYNGAV85Jy7z8xuBW4FbgF+AmR6j2OAx72fIiJBSUmI4YUxg7h2ylxuf30h0ZEBLuyvsaoOxcHGnsoP1YqdcxvwrvVwzhWZWS6QDgwHTvYWmwz8h+rSGA4846oPwswws+Zm1tb7HBGRoMRFR/LkiP5cNWk2t7y6gMTYSM7o3sbvWGEj2GMaIWVmHYG+wEwgbZ8i2Aikec/TgTX7vG2tN01E5JDEREYwbkQ2PdOTueH5eXyRt8XvSGHD99IwswTgVWCsc27nvvO8rYpDuiWXmY0xsxwzyykoKKjDpCLSmMTHRDLpygF0Sonn6mdy+GrNDr8jhQVfS8PMoqgujCnOude8yZvMrK03vy2w2Zu+Dmi/z9szvGnf45wb55zLds5lp6ZqeGQRObDmcdE8M3ogrRJiGDVxFss3FfkdqcHzrTSs+vZ/E4DcH9wB8E2+u9p8JNUXFu6dfoV3FtUgoFDHM0SkttKSYnlu9DFERwS4fMJM1mzT6Lg18XNL43hgBHCqmX3lPc4E7gNOM7PlwFDvNcA7wAogj+qh2q/1IbOINEJHtIrj2dHHUFJexYgJMyku1SCHBxLUFeHhSleEi8ihmLFiKxePm8Fvhh7FjUMz/Y7jm7q4IlxEpNEbdGQrhnVvw7hPv2HLrlK/4zRIKg0RkX38flhXSiqqePTjPL+jNEgqDRGRfXROTeDnA9ozZWY++Vs1mPcPqTRERH5g7JBMIgMBHvhgmd9RGhyVhojID7ROiuUXJ3TirfnrWbi20O84DYpKQ0RkP8aceCQt4qK4/72v/Y7SoKg0RET2IzE2ihtOzWR63hY+W64hifZSaYiIHMBlg44go0Uz7nv3a6qqGu81bYdCpSEicgAxkRH8/oyuLF6/k7cWrPc7ToOg0hARqcE5vdqR1TaJv76/lNKKSr/j+E6lISJSg0DAuO3Mbqzdvoc/v7WExjz0UjBUGiIiB3FCZirXnNSZKTNXM2H6Sr/j+Opg9wgXERHg5jO6snpbMXe/k0tGiziG9Wiat4jVloaISBACAePBi/rQO6M5Y1+cx/wmeqc/lYaISJBioyIYPzKblIQYRk/OYe32pnfDJpWGiMghSEmIYeKoAZRWVHLVpNnsLCn3O1K9UmmIiByizLREnri8PysKirluylzKK6v8jlRvVBoiIofh+C4p3HN+Tz5bvoU/v7XE7zj1RmdPiYgcposGtCevYBfjPl1Bz/RkLhrQ3u9IIactDRGRWrj5jK4M7pLCH/61iK+awBlVKg0RkVqIjAjwj0v60jophmuenUNBUeO+t7hKQ0SkllrER/PkiP7s2FPW6A+MqzREROpA93bJ3P/TXsxatY27/53rd5yQ0YFwEZE6MrxPOgvXFjJ++kq6t0viZ9mN78B42G1pmNkwM1tqZnlmdqvfeURE9nXrT7pxXOdW3PGvRY1yqJGwKg0ziwD+CfwEyAIuMbMsf1OJiHwnMiLAo5f2IzUhhltfW9johlIPq9IABgJ5zrkVzrkyYCow3OdMIiLf0zI+ml8P6ULuhp18uWKr33HqVLiVRjqwZp/Xa71p3zKzMWaWY2Y5BQW6GbyI+GN4n3RaxUfzdCO7/0a4lcZBOefGOeeynXPZqampfscRkSYqNiqCywZ14MPczawo2OV3nDoTbqWxDtj3dIQMb5qISIMzYlAHoiMCTPx8ld9R6ky4lcZsINPMOplZNHAx8KbPmURE9is1MYbhfdrxypy17Nhd5necOhFWpeGcqwCuB94HcoGXnHOL/U0lInJgo0/oxJ7ySp6ftdrvKHUirEoDwDn3jnPuKOdcZ+fc3X7nERGpSbc2SQzuksLkL1ZRVhH+w4uEXWmIiISb0YM7sWlnKe8s3OB3lFpTaYiIhNhJR6XSOTWeCdNXhv3FfioNEZEQCwSMqwZ3YuG6Qmav2u53nFpRaYiI1IML+mbQPC6K8Z+t8DtKrag0RETqQbPoCC4/pgPTcjeRv7XY7ziHTaUhIlJPrji2A5EBC+uL/VQaIiL1pHVSLOf2TmfKzHw+yt3kd5zDotIQEalHd56dxdFtk7jmuTlMWxJ+xaHSEBGpR8lxUTw7+hiy2iXzq+fm8N6ijX5HOiQqDRGRepbcLIpnRw+kZ0Yy1z8/l3fD6KI/lYaIiA+SYqN45qqB9G7fnOtfmMfbC9b7HSkoKg0REZ8kxkYx+aqB9DuiOTdO/Yo35zf84lBpiIj4KCEmkklXDqR/hxaMnTqPuasb9hXjKg0REZ/Fx0QyYWQ2bZJiueWVBZRWVPod6YBUGiIiDUBibBR3n9+T5Zt38dgn3/gd54BUGiIiDcQp3VozvE87HvtPHss2FfkdZ79UGiIiDcidZ2eREBPJza8soLKq4Q2jrtIQEWlAWiXE8MdzuvPVmh1M/mKV33F+RKUhItLADO/TjpO7pvLX95eyZttuv+N8j0pDRKSBMTPuPr8nAYPbX1/YoO72p9IQEWmA0ps34+Zh3fhs+RZem7vO7zjfUmmIiDRQIwZ1ILtDC/789hIKikr9jgOoNEREGqxAwLjvpz0pKiln4ucr/Y4DqDRERBq0Lq0TGXJ0Gi/OXtMgrhT3pTTM7K9m9rWZLTCz182s+T7zbjOzPDNbamZn7DN9mDctz8xu9SO3iIgfRgzqwNbiMt5d6P+9N/za0pgG9HDO9QKWAbcBmFkWcDHQHRgGPGZmEWYWAfwT+AmQBVziLSsi0ugN7pJCp5R4nvlyld9R/CkN59wHzrkK7+UMIMN7PhyY6pwrdc6tBPKAgd4jzzm3wjlXBkz1lhURafQCAePyQR2Yu3oHi9YV+pvF17VXuwp413ueDqzZZ95ab9qBpv+ImY0xsxwzyykoKAhBXBGR+ndhvwxiowI8NyPf1xwhKw0z+9DMFu3nMXyfZe4AKoApdbVe59w451y2cy47NTW1rj5WRMRXyXFRnNcnnX99tY7C3eW+5YgM1Qc754bWNN/MRgFnA0Pcd5c7rgPa77NYhjeNGqaLiDQJI47twNTZa3hl7lpGD+7kSwa/zp4aBtwMnOuc23dglTeBi80sxsw6AZnALGA2kGlmncwsmuqD5W/Wd24RET91b5dMvyOa89yMfKp8GgHXr2MajwKJwDQz+8rMngBwzi0GXgKWAO8B1znnKr2D5tcD7wO5wEvesiIiTcoVx3Zk5ZZipudt8WX91pAGwqpr2dnZLicnx+8YIiJ1prSikuPu/Zh+HVrw1BXZIVmHmc1xzu33wxvC2VMiIhKkmMgIfj6gPR/lbmLdjj31vn6VhohImLlsUAcApvhw+q1KQ0QkzKQ3b+bbeFQqDRGRMLR3PKopM1bX63pVGiIiYWhwlxT6e/faGD1pdr3dFlalISIShgIBY+qYQdx+Zje+XLGVoQ/+l4c/XE5JeWh3V6k0RETCVFREgDEnduaj357E0Kw0HvpwGWf8/VP+s3RzyNap0hARCXNtk5vxz0v78dzoY4gIGKMmzua65+eG5KrxkI09JSIi9WtwZgrv3ngCE6avZHdpJYGA1fk6VBoiIo1ITGQE157cJWSfr91TIiISNJWGiIgETaUhIiJBU2mIiEjQVBoiIhI0lYaIiARNpSEiIkFTaYiISNAa9e1ezawAqM1dSlIAf27E65+m9p2b2vcFfeemojbfuYNzLnV/Mxp1adSWmeUc6D65jVVT+85N7fuCvnNTEarvrN1TIiISNJWGiIgETaVRs3F+B/BBU/vOTe37gr5zUxGS76xjGiIiEjRtaYiISNBUGiIiEjSVxn6Y2TAzW2pmeWZ2q995Qs3M2pvZJ2a2xMwWm9mNfmeqL2YWYWbzzOxtv7PUBzNrbmavmNnXZpZrZsf6nSnUzOw33p/rRWb2gpnF+p2prpnZ02a22cwW7TOtpZlNM7Pl3s8WdbEulcYPmFkE8E/gJ0AWcImZZfmbKuQqgN8657KAQcB1TeA773UjkOt3iHr0MPCec64b0JtG/t3NLB34NZDtnOsBRAAX+5sqJCYBw34w7VbgI+dcJvCR97rWVBo/NhDIc86tcM6VAVOB4T5nCinn3Abn3FzveRHV/5Ck+5sq9MwsAzgLGO93lvpgZsnAicAEAOdcmXNuh6+h6kck0MzMIoE4YL3Peeqcc+5TYNsPJg8HJnvPJwPn1cW6VBo/lg6s2ef1WprAP6B7mVlHoC8w0+co9eHvwM1Alc856ksnoACY6O2SG29m8X6HCiXn3DrgAWA1sAEodM594G+qepPmnNvgPd8IpNXFh6o05FtmlgC8Cox1zu30O08omdnZwGbn3By/s9SjSKAf8Lhzri9QTB3tsmiovP34w6kuzHZAvJld7m+q+ueqr62ok+srVBo/tg5ov8/rDG9ao2ZmUVQXxhTn3Gt+56kHxwPnmtkqqndBnmpmz/kbKeTWAmudc3u3Il+hukQas6HASudcgXOuHHgNOM7nTPVlk5m1BfB+bq6LD1Vp/NhsINPMOplZNNUHzd70OVNImZlRvZ871zn3oN956oNz7jbnXIZzriPV/48/ds416t9AnXMbgTVm1tWbNARY4mOk+rAaGGRmcd6f8yE08oP/+3gTGOk9Hwm8URcfGlkXH9KYOOcqzOx64H2qz7R42jm32OdYoXY8MAJYaGZfedNud869418kCZEbgCneL0QrgCt9zhNSzrmZZvYKMJfqswTn0QiHFDGzF4CTgRQzWwv8EbgPeMnMRlN9i4iL6mRdGkZERESCpd1TIiISNJWGiIgETaUhIiJBU2mIiEjQVBoiIhI0lYZIHfJGkb3We97OO91TpNHQKbcidcgbu+ttb0RVkUZHF/eJ1K37gM7eRZLLgaOdcz3MbBTVo4zGA5lUD6IXTfVFlaXAmc65bWbWmeqh+VOB3cDVzrmv6/tLiByIdk+J1K1bgW+cc32A3/9gXg/gAmAAcDew2xs48EvgCm+ZccANzrn+wO+Ax+ojtEiwtKUhUn8+8e5XUmRmhcBb3vSFQC9vlOHjgJerh0kCIKb+Y4ocmEpDpP6U7vO8ap/XVVT/XQwAO7ytFJEGSbunROpWEZB4OG/07mGy0sx+BtWjD5tZ77oMJ1JbKg2ROuSc2wp8bmaLgL8exkdcBow2s/nAYhr5rYYl/OiUWxERCZq2NEREJGgqDRERCZpKQ0REgqbSEBGRoKk0REQkaCoNEREJmkpDRESC9v9m9VTr3O/UlgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "pa = -20\n", "pb = 90\n", "pc = 800\n", "\n", "t = np.linspace(0, 10) \n", "y = pa*t**2 + pb*t + pc + np.random.randn(np.size(t))*15\n", "\n", "\n", "plt.plot(t, y)\n", "plt.xlabel(\"time\")\n", "plt.ylabel(\"height\")\n", "plt.savefig(\"fig-res-missle_taj.pdf\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.1 如何得到更新项?\n", "\n", "$$\n", "L = \\sum_{i=1}^N (y_i - at_i^2 - bt_i - c)^2\n", "$$\n", "\n", "\\begin{eqnarray}\n", "\\frac{\\partial L}{\\partial a} & = & - 2\\sum_{i=1}^N (y_i - at_i^2 - bt_i -c) t_i^2 \\\\\n", "\\frac{\\partial L}{\\partial b} & = & - 2\\sum_{i=1}^N (y_i - at_i^2 - bt_i -c) t_i \\\\\n", "\\frac{\\partial L}{\\partial c} & = & - 2\\sum_{i=1}^N (y_i - at_i^2 - bt_i -c)\n", "\\end{eqnarray}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.2 程序" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "epoch 0: loss = 2.47402e+07, a = -3.30002, b = 6.91887, c = 4.30695\n", "epoch 500: loss = 1.09272e+06, a = -29.2101, b = 226.526, c = 431.351\n", "epoch 1000: loss = 366835, a = -24.1553, b = 158.837, c = 599.169\n", "epoch 1500: loss = 135985, a = -21.2759, b = 120.282, c = 694.733\n", "epoch 2000: loss = 63728, a = -19.6362, b = 98.3262, c = 749.154\n", "epoch 2500: loss = 41779.5, a = -18.7024, b = 85.8231, c = 780.145\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4OklEQVR4nO3dd3xN9//A8dcnsoiV2FuMWrGjVUFdtIISNRJbdejQlrY/oUOpamtUldaoraW4YlM7V61vCapqFZEaoQQRSciSz++Pe90mxMy4Se77+XjcR+4953PPeZ+E9/3czzpKa40QQgj74mDrAIQQQmQ9Sf5CCGGHJPkLIYQdkuQvhBB2SJK/EELYIUdbB/AoihYtqitWrGjrMIQQIkfZv3//Fa11sbT25YjkX7FiRfbt22frMIQQIkdRSp253z5p9hFCCDskyV8IIeyQJH8hhLBDOaLNXwiR+RITEzl//jxxcXG2DkU8JldXV8qWLYuTk9Mjv0eSvxACgPPnz1OgQAEqVqyIUsrW4YhHpLXm6tWrnD9/Hk9Pz0d+nzT7CCEAiIuLo0iRIpL4cxilFEWKFHnsb2yS/IUQVpL4c6Yn+btJ8hcPdu4cTJsG8fG2jkQIkYEk+Yu03b4N338PNWvC229Djx6QlGTrqEQulydPHurVq4eXlxcdOnTg+vXrT3ScefPm8c477zy0XMWKFbly5coDy3z11VdPFEN2J8lf3OvwYWjaFN57D3x8YORIWLEC+veH5GRbRydysbx583Lw4EEOHz6Mh4cHU6ZMsXVIkvyFDZ07B7NmwcGDmZt84+Lg00+hfn04dQoWLoT162HECPjyS1iwwPwtQO7+JrLAs88+S3h4OAChoaH4+vrSsGFDmjVrxvHjxwFYs2YNzzzzDPXr16d169ZcunTpgce8evUqL7zwArVq1eK1114j5Z0MO3XqRMOGDalVqxYzZswAYNiwYdy6dYt69erRq1ev+5bLiWSoZ3Z2p+nl008hNta8rVgxaNUKWreG55+H8uXTd44rV+Cvv+DQIZg6FU6cgL59YcIEKFr0v3IffwzR0TBmDOTPD+PHg3QO5l6DB5srGxmpXj347rtHKnr79m22bt3Kq6++CsCAAQOYPn06VatWZc+ePbz99tsEBwfTtGlTfv/9d5RSzJo1i3HjxjFhwoT7Hvfzzz+nadOmfPbZZ6xbt47Zs2db982ZMwcPDw9u3bpFo0aN6NKlC2PGjOGHH37gYIrfRVrlihQp8iS/EZuS5J9dHTgAAwbA/v3Qti2MGgVHj8KWLbBlC+POLabRaDC4VDN3yBoMmMJMjN89niFNhmDwNFgPlWr7v67mJpxDhzBdCWF8tWsM2Q2Gf4CqVWHjRkxVnQj5ew6BRQNTx/TVVxATY/5gKFDA/I1AiAx0p5YdHh5OjRo1eP7554mJiWH37t1069bNWi7eMgDh/PnzBAQEcPHiRRISEh46zn379u0sX74cgPbt2+Pu7m7dN3nyZFasWAHAuXPnOHnyZJpJ/VHLZXta62z/aNiwobYb0dFav/++1g4Oeqxvfh08Z7jWycnW3cGng/XYHWN0cPAcXfRzNx3crJzWzs46eO5nuui4onrC7gm66LiiOvh0sLW8dfvn+XVwJQdz+daVddHhLnrClx100dGFdPC+IK2Tk63lB6weYD3G3ece+15DHVwRrSdMSL1v59is+R2JTHH06FFbh6Dd3Ny01lrHxsbqpk2b6kmTJumoqChdsmTJNMs/99xzetWqVVprrU0mk37uuee01lrPnTtXDxw48J7ydevW1aGhodbX7u7uOiIiQptMJu3j46NjY2OtxzWZTKliunOO+5WztbT+fsA+fZ+8Km3+thIbC6GhsGsX46b3wTRxkLl5p1YtmDgR08B2hL7aBf8r0zD9sw2ALae30HVpV4rkK0pyhfK80XQwHXyv0fINV9qfGkXzxDIcizhG7eK18V3oS5XJVWizoA0VC1Vkw7rvKH8+hja9NZXHlsb3uXPUqezDkVrF8KncgvYb+vDcvBZ0XNSRNxq+QRWPKnQ2dmbRX4u4mXgTU5gJ/yB/GpV5mkbvjsG/tzOm7z+EokUxNSmN/8wXCP3pO0y9feDll2H2bLhpft+4XeNs+qsWOU++fPmYPHkyEyZMIF++fHh6erJ06VLAXGH9888/AYiKiqJMmTIAzJ8//6HHbd68Ob/88gsA69evJzIy0nocd3d38uXLx/Hjx/n999+t73FyciIxMfGh5XIaafbJKhcvMi7ofRodjsSw7ghYOrJMFSHUC8bXgPkroXSl6qwZ/Spjzy3G97Yv1YpU4/mfn8c5jzO3km4B8Nqa11Id2lQE0LAl+k/cDoTiVqQk7q7uhEaGUtqtJHlOniQ6MgqXUqUoUciB01H/UDxfcc7eOMvfV//mVtItEm8nsv3sdgC+3PGl9dg9l/cEQKGo5F6Jb3//lvIFy9OpYW/8HBfif7MYK11Os/RkQ4iNxd9rD8ZNxzHMn49p0mD8O9/G2GZWpv96Re5Tv3596tSpw6JFi1i4cCFvvfUWo0ePJjExke7du1O3bl1GjhxJt27dcHd3p2XLloSFhT3wmCNGjKBHjx7UqlWLJk2aUN7SZ+br68v06dOpUaMG1apVo3Hjxtb3DBgwgDp16tCgQQPmzJlz33I5jdI5YOSGt7e3zkk3cxm3axyNSjfCEKZh9WrYsgVT7BEWe8HyWoq5ET4UqezF8rz/MCUmGK+CVTkTd4nLcanHG7s5uVGxcEXib8dz6topmpdvToBXAMXyFaO4W3FOXz/NkE1DGNBwADP3z8R49hkMU9Zheqst/p4hvFWlB9NCpmJcojEETsX0wlP4B/nzlvdbTNs3DWNXIwZPg7VW/0bDN5i+bzrfvPANVTyqcPXmVWYdmMXak2vxLuVNqQKlOBt1lnM3znHt1rVUsZbKX4rqRauT3zk/wWHBvOThw/rwbSz9JYmQUsk0qtAEQ//PoWVLiI/H9PcGQs7tIdCzt7kfwcUF6tYFB/kyaivHjh2jRo0atg5DPKG0/n5Kqf1aa+8033C/9qDs9Miubf5jd45Ns118wOLeuuinzjq4IjqmgIue1Ocp7TbSWbf64Wld7ttympFYH06jnHT96fV1r2W9tGGeQTMS/dqq1/SV2Cs6OUUb/PDg4Wm25d/Ttj/8eV10CDq4c32tixXTwV5uuujoQg/uC7jPce93bq21XndinXYf4679jf7a7Us33ebnNrrxrMa68JjCqa7Pc0J53eyzcjrfJ+ivfdCX3dDBFdFFh6AHvGh+rrE8SpTQwW+11WN/7Kt1TExW/AlFCtmhzV88ucdt87d5Yn+UR3ZN/vck4NCt2mOUmx7u66r9eiidZ6RKlQjLTCij/Rb56RZzW2hGot/59R2ddDsp1bFSJtr7Jfg7natpffC0XdBWB383SGultK5cWetjx/7bfr/yaX2ArR7w0A+FtPZtDd2qPcZ66J7Lemq3L930c3Of01UmV0n1e8gzUum2X9bUgd+21e6j3PTWRV9p/fPPOrh/C100UJk/EFxctG7XTutp07QOD8/wv524lyT/nO1xk780+9wtOdn8cDR3h1ibcO4aOhlyIYRAn0CWHV1G/1X9qZK/PIeuHuG2Zei7m2M+iucvQdj1MAJqBTCxzURKFShlbWJJ2fQC4B/kf08zTOfqnenu1f2+536gw4ehXDkoVOiJfg0Pum4gzX2LDy9m+fHl91yHsauRBqUa8O76d/n50M9U9ajKlZtXiIwzd7YpFNWLVufcjXPMaDuNc39up9FfVzGsOAinT5uP71udkCYVCOzyLdSoIXMMMoE0++Rsj9vsI8k/peho6NjRPLv1p5+sY+dTJuYtp7fQ1diV9lXbc+zKMf749w/r26teUwwo68dzPYYRlRBNj2U9Mj/JZyMP+sBoVLpRqg+9xV0WUzJ/SXae3cmUkCn8dfkv63vKFSxHxM0Ihjf7lPcKvUDImh/xvzkf46Ik83yEKlWgUycYMgSKF8/y68ytJPnnbJL8n9C4raNoNHEJhg1/m2vMZ85gGupPSNs6VC/uRe8VvSlTsAx/X/kbjSaPykOTQl48dfA8y0pc5Z0rnkyvEoUxIAjI3Un+cd39AZryNZh/V282fJMpIVPo4dWDsOthbD29lYTkBBxwII9DHt5p9A6Blfvw0+ovaLTnHIaVf4KXF2zfjikiJNf/DrOCJP+cTTp8n8S1azrYt5q5o3T+SK1jYvSy957X+T9Cew/Kpx0/d7S2V9ebXk8H7fhRX+/TzdxxOdRBB8/8+NEmSNnpJKj7dow/oF8hLjFO917eWzMS7T7GXTMSrUYqXfOHmtrtSzf90/wPtHZw0ME9Gt/TGS2ejLT552wyyeshxu0ahynM9N+Gq1cxBTxDyM1TzKz9MR3Dv6HCzJp08dhMjAtcc4ijy9FkCqu8fPrsMM7/ewKP1wdRyLiakD4tMfZbg+G1L0EpDJ4GjF2NVPaonKp2D2DwNNhtzTTQJzDN30dlj8rWbwN3thm7Ggm5EMLuc7vZcGoDw5sPJ49DHmZ3nM3nLT7H1cmV2MRY+oZ9S8kRbnTw/J0p/zYkJHxv6r8r/DfBLCEBhg8HPz/z5LOrV7Ps2sXjubOk853HmDFj7lt25cqVHD161Pr6s88+Y8uWLemO4fr160ydOvWx3zdy5Ei++eabh5bLnz9/ppz/sd3vUyE7PTKy5p9qpMrly3rD85664EfoJuNraKdRTtYafrM5zfQfF//QW/cu0UU/cTKPQClUyFzb/9RZB+9akGExidQeNJpIa63/ifxHv/DTC9a/lRqBrjuqtM7/VX696tiq1O/Z+bPW3t5ag9Zlyph/Ojpq7eur9ezZWl+9arPrzG4ep+Z/v29z6f12m3IphYfp16+fXrp0abrOl5awsDBdq1atx37fiBEj9Pjx4x9a7mHX+KTnt0nNXyn1vlLqiFLqsFJqkVLKVSnlqZTao5Q6pZRaopRytpR1sbw+ZdlfMSNieFQGTwNLuiyh8+JONBxdnrZNwrjhAqdVJH7V/CjsWphPm33KsSvHiLwVyb6EfzC+sh7DuxOgXj0MMzZhfGUDITo8K8O2KyEXQu77jQDgdORpDvx7gOHNh+Pu6k7fyPLcjLhATEIMfkv8qDmlJi8teYnF7gMI+exVTHHHYPly89LY+/djGtKVcU574dVXoUQJ6NoVbt605SXnOHc68O9827Iu/1G6Uaacb9iwYdSsWZM6derwf//3f+zevZvVq1czZMgQ6tWrR2hoKC+//DJBQeY+t4oVK/LRRx9Rr149vL29OXDgAG3atKFy5cpMnz4dgJiYGFq1akWDBg2oXbs2q1atsp4rNDSUevXqMWTIEADGjx9Po0aNqFOnDiNSLGj45Zdf8tRTT9G0aVP+/vvvNGMPCwvj2WefpXbt2nz66afW7Y96/vuVS7f7fSo86gMoA4QBeS2vjcDLlp/dLdumA29Znr8NTLc87w4sedg5nqTmn6pmcumS1lu36tXjX9N+Qytorw/zWmuNXuM89YaTG/Tm0M0PrG2K7CHNbwVji+qtvtV0iKeLfnpiLevftsJg9MuvF9MeXxW+9+8aulXrffu0/vBD85yIAQNseVnZQsqa46D1g/Rzc5974KPOtDraaZSTLj+xvHYa5aTrTKvzwPKD1g96aAwODg66bt261sfixYv1lStX9FNPPaWTLQscRkZGaq3vrfmnfF2hQgU9depUrbXWgwcP1rVr19Y3btzQly9f1sWLF9daa52YmKijoqK01lpHREToypUr6+Tk5Htq3hs3btSvv/66Tk5O1rdv39bt27fXv/32m963b5/28vLSsbGxOioqSleuXDnNmn+HDh30/PnztdZa//DDD9aa/6Oe/37lHvT3u4MH1Pwzam0fRyCvUioRyAdcBFoCPS375wMjgWmAn+U5QBDwg1JKWQLNMI1KN8Lf2JUli5Nwi7jByBawoQqQF56KdSV/shNvV+nBnH9/xTmP8wNrm3e3VwvbSfPv1M1ISJ1gGn0+k9P/HuWj0x78UOUaHgVLMK/wJVSCwnehLy9Vf4mtYVtTvZ+GDSFPHhg3Dtq0gc6dbXh1OYu7q7t1yY/yhcrj7ur+8Dc9xJ07eaWUlJSEq6srr776Ki+++CIvvvjiIx2rY8eOANSuXZuYmBgKFChAgQIFcHFx4fr167i5ufHxxx+zfft2HBwcCA8PT/NmMJs2bWLTpk3Ur18fMNfYT548SXR0NC+99BL58uVLdb677dq1i2XLlgHQp08fhg4dCpgr3o9y/vuVK1my5CP9Hu4n3clfax2ulPoGOAvcAjYB+4HrWus7N309j/kbApaf5yzvTVJKRQFFgAffSPMxGTwNTG4ziRdu9bVOvGpfrhUdanfl023DWW1JAL53DUO8+xiS+LOXtDrN7/yN/P2mYFyYD8MVZ54P+Ab/02P4uc3PHIs4xuQ9k1lyZAnuru4cunSIHWd30Kx8M/N7v/gCgoMxjexHiMs+Atvnztv2PY7vfL97aJk7TT3Dmw9n2r5pjHhuRKb8f3F0dGTv3r1s3bqVoKAgfvjhB4KDgx/6PhcXFwAcHBysz++8TkpKYuHChURERLB//36cnJyoWLEicXFx9xxHa81HH33EG2+8kWr7d494YxoAlcakxEc9/6OWe1zpbvNXSrljrs17AqUBN8A3A447QCm1Tym1LyIi4omO0bV2ABU8zDd3CGwSyNpXthCVcOOB7ckiZwq5EIIxYBmGTSfh+HEMXT7E2NXIhegLtK7UGlcnV7rU6EJ0QjSDNw5mzM4xtPulHQsOLQBnZ0wTB+HfNoZGM9aa76AmHijlXI1RhlEYuxpT9QFkpJiYGKKiomjXrh0TJ060LudcoEABoqOjn/i4UVFRFC9eHCcnJ0wmE2fOnEnzuG3atGHOnDnExMQAEB4ezuXLl2nevDkrV67k1q1bREdHs2bNmjTP4+Pjw+LFiwFzIn/c89+vXHplRLNPayBMax0BoJRaDvgAhZVSjpbaf1ngTg9pOFAOOK+UcgQKAfeMvdNazwBmgHmS15MEtvPsTm7E37DWTHyr+N635ig1/Jztgd8I7ppg1tnYGe/S3pjCTPRZ0YdRv43icuxlVlQehmHcGPP9ij/7LKsvIUfJrGbSO3fyusPX15dBgwbh5+dHXFwcWmu+/fZbALp3787rr7/O5MmTrR29j6NXr1506NCB2rVr4+3tTfXq1QEoUqQIPj4+eHl50bZtW8aPH8+xY8d49tlnAfNQzQULFtCgQQMCAgKoW7cuxYsXp1GjtDu7J02aRM+ePRk7dix+fn6Pff6hQ4emWS7d7tcZ8KgP4BngCOa2foW5ff9dYCmpO3zftjwfSOoOX+PDzvEkHb4PGy4o7MODhiReuHFBN5vTLNXw3l9fb6HHNEUHr56c5ntyM5nklbNl+VBPrfUezB23B4C/MDclzQCGAh8opU5hbtO/c6fk2UARy/YPgGHpjSEtDxsuKOzD/SaYBfoEcvzKcY5dOcZQn6G4Oblx7Mox2pXZxsynHfH7fRBbDq0EMn8YoxC2IGv7CLt0v/WG+tfrz+o/jfwde4Y8Gl6s1I5dl/amOSAgt5G1fXK2x13bR27jKOzSg74ZHvkglJVf92Ng5EJWhf1K0ViI6NeFsaXy8jRlMOStCW5u5o5i14uEOEUQ6NQCnJ2hdWvz8NEcSmud5sgUkb09SSVeav5CpMF0eiv+izrTAk9WJB7ittJUisvLlTzxrNjsQcvTYCpxE//2NzEuBcOdW8e6uMCKFdC2rU3jfxJhYWEUKFCAIkWKyAdADqK15urVq0RHR+Pp6Zlqn9T8hXgMpjAT/su6Y+y50noPh85LOhNfoCA3osNp8+J1AmoFsDF0I8auazGMNZiHh16+DO3bmxeQMxrN9xzIQcqWLcv58+d50qHVwnZcXV0pW7bsY71Hav5C3OV+N6XZE76HEm4leH/j+0TFR1HZvTIru6/Eq7jXf2++fh18fWHfPli4EAICsv4ChLCQm7kIkUHudAzXK1GPrWFbAfAu7c2QJkPoVqubuVB0NKbePoRcO0zga3OhXz8bRizs2YOSv92t5y/Ek0o5Qmhz382sCFiBq6Mr+y/sxz/In74r+hKTEIPpyj78fS7QqER9ePll+PFHW4cuxD2kzV+IR3T3CCG/6n6s67mODac2sPfCXn4+9DMrjq9AoVgZsBLDe00griu8+SbEx8N779n4CoT4jzT7CJFBXln1CnMPzgWgeYXmNCjZgI6V22IYOt08AmjRIkzPlJD7DYssI80+QmQyU5iJNSfW8EmzT8jvnJ+DFw8yac8k2hn9WDOqFzRtiumzvvgv6iwzhUW2IMlfiHRK2RcwuuVoVndfjVMeJzo+1ZH4pHj8lnWj42tu+He5jdGYjCG5vK1DFkKSvxDpldZs4aXdltKkfBP2D9hPmYJlWPPPRvJ5lKR0tDLPBYiMtHHUwt5J8hcinR60eNz1uOvEJcXRvmp7zt28iFffGF545gSbXm0BCQnW8qYwE+N2jcviyIU9k+QvRCZJ2Ry0tudalnZbSp48jmyupGlX+xDff9gMtJZVQ4VNyFBPITLJ3c1BXWp2wSOvBz/9+RO//rmU94ruZf6osvyTN56l3Zbm+lVDRfYiyV+ITPKgu8bFtJ1M65Ge7HG9QMFYZ5IT0n9PViEehzT7CGEDIRf2Eequ6B9fk1gSaL24HQMX9+GL7V/ccx9c6Q8QmUGSvxBZzNoX0M3InK+OsLbycPImwtTjC5i6YyIvLXnJ+gEg/QEis0jyFyKL3d0X4Nt3FOvaL+TNM8VwuxxJVHwU7X9px0dbPkp1tzEhMpIs7yBEdpGQQGzg+ww7MZUfnjFver2YLzPa/AAVKoCjdNGJxyPLOwiREzg74/bdFDr3+JxC8Yr88TDz8gbeHlSF5Hx5oXp16NgR/vc/W0cqcgFJ/kJkI6YwE/7/fs+K17YQ9sZRmhTyYlojKPeRK0ufLQR795rvEHb5snQEi3SR5C9ENmLtD6jUkqIVarBz8CE+aPwB/+a5SUDFEL6Y3BWiojAN7iQdwSJdpM1fiBzg+JXjdFzUkZPXTtLQoSxnos9jLPcBhncn2Do0kY1Jm78QOVz1otU5/PZhGpdpzP7k8zg6OVNuzFQ4ccLWoYkcKkOSv1KqsFIqSCl1XCl1TCn1rFLKQym1WSl10vLT3VJWKaUmK6VOKaUOKaUaZEQMQuR2u87u4lTkKQJqBXDJOZE6/eNY+n/tIDHR1qGJHCijav6TgA1a6+pAXeAYMAzYqrWuCmy1vAZoC1S1PAYA0zIoBiFyrZSLxC3uupiFnReS6JQH/4ahdBrtRXxSfKqy0hEsHibdyV8pVQhoDswG0FonaK2vA37AfEux+UAny3M/4Cdt9jtQWClVKr1xCJGb3T0xrEftHqzr/Su1491Z5XCCOt9W5XTkaZkRLB5Zujt8lVL1gBnAUcy1/v3AICBca13YUkYBkVrrwkqptcAYrfVOy76twFCt9b67jjsA8zcDypcv3/DMmTPpilOIXOnGDUb3qcBn9a/j7OiCi6OL+ebxMiNYkPkdvo5AA2Ca1ro+EMt/TTwAaPMnzGN9ymitZ2itvbXW3sWKFcuAMIXIhQoW5NPAtbyzF+Jvx3Mj/gY7z+4kWSfbOjKRzWVE8j8PnNda77G8DsL8YXDpTnOO5edly/5woFyK95e1bBNCPAFT6QQWNc7PsP/lwSUJPtv2GbWn1GLN32tSl5O+AJFCupO/1vpf4JxSqpplUyvMTUCrgX6Wbf2AVZbnq4G+llE/jYEorfXF9MYhhD2ydgT3Wc3X00/x69lmuCXAsYjjvLS4E3P/mJuqnPQFiDsyZJKXpd1/FuAMnAb6Y/5gMQLlgTOAv9b6mqX9/wfAF7gJ9L+7vf9uMslLiLSN2zWORqUbpWrjNy2fwJKVX7KkdCTXXaFruRfYdu2ArA5qhx7U5i8zfIXIjZKSCJ/yNc+cG0F4AU3zIg3YNnAf5rqXsBcyw1cIe+PoyImOTYkv5k6dq45sv3qAljN9uJl409aRiWxCkr8QuZC1L8A/iIOvhDDgYB62Xfgf9afVI/yGjK8QkvyFyJVSTgpT9erxY78gvgiG0KunqDu9Lvsu/NeMKqOA7JMkfyFyoUCfwNSdu5068Wn7MUxfo4m8dY0ms5tgPGKUUUB2TDp8hbAXWsPLL7N8708EdM9DErfJ55SPNd3X0LJSS1tHJzKBdPgKIUApmDGDzkV8GLLbvOlm4k0WHV5E4m1ZGdTeSPIXwp64uGD6/gNmNtB8ut+NfHlcmfXHLDos6kB0fLStoxNZSJK/EHbEFGbCf8sbGJ+fyRe7XFi70o38jvnYHLqZaj9Uw3jEeE956QzOnST5C2FHrKOAWr4Cv/2GIdyJ1Usd6V++I5FxkfRY1oPZB2YDsiREbicdvkLYs5MnoVUriI7m4JJJtD74AdduXaNXnV5sOLVBloTI4aTDVwiRtqpVYccOKFKEel0G8kedqRTNV5QFhxbwXIXnJPHnYpL8hbB3FSqYPwDKl+fUoN4kJ8RRtmBZlh1bxuANg20dncgkkvyFEFCqFKaFX+LfRbN03k3+NnnR+HZpJu2ZRL+Z7dE3btg6QpHBJPkLIQAIiT2BscdyDC36k+/MBXaMi6DNSfjpwq90eqMQtyuUg4ULARkFlBtIh68QIm1JSejTp+mx/lWWXN/JcxFubFzixO5tP+P/a3/pDM4BHtTh65jVwQghcghHR9RTT7H4qR24r32L6funU6Mn3AgKYGnftZL4czhp9hFCPNS0F6fR8amOhHlA3uhbNHCrbOuQRDpJ8hdCPJQpzMTu87vpUrIl5/Nrnp7WgCs3r9g6LJEOkvyFEA9kvTFMVyNBb2zlqwvVOaGv0ujHhlyMvmjr8MQTkuQvhHiglDeGAfjozYV8swnO3jiH90xvzlw/Yy0ro4ByDkn+QogHuufGMA0a8GFxPyaZXLgYfZFGMxtx8upJWQsoh5HRPkKIxzdiBO80WIWz4UXevLmO+j/WxzmPM8v8l8kooBxCav5CiMdXvz506sSAiTt4o/bLxCbGknA7gZL5S9o6MvGIJPkLIZ7MiBGY3KMIOrSEt7zf4mbiTZrMacLRiKO2jkw8ggxL/kqpPEqpP5RSay2vPZVSe5RSp5RSS5RSzpbtLpbXpyz7K2ZUDEKIrGMqFIl/b2eMQYqpPl8x128uN+Ju0GR2Ew5fPmzr8MRDZGTNfxBwLMXrscBErXUVIBJ41bL9VSDSsn2ipZwQIocJuRCCscUUDEdi4ZNP6OfZibmd5nIr6RY+c3z469Jf1rIyCij7yZDkr5QqC7QHZlleK6AlEGQpMh/oZHnuZ3mNZX8rS3khRA4S6BOIodVr0LcvTJ0KxYrRd+gvzHLrSUx8DE3nNOXPf/+UUUDZVEaN9vkOCAQKWF4XAa5rrZMsr88DZSzPywDnALTWSUqpKEv5VNMFlVIDgAEA5cuXz6AwhRAZbu5ceOMNWLECVqygz8aNqDrQ76UbPPtjI1xd88sooGwo3TV/pdSLwGWt9f4MiMdKaz1Da+2ttfYuVqxYRh5aCJGRHBygSRMYP958W8i//qJ3ty94+3QRbpFIfPxNShUoZesoxV0yotnHB+iolPoHWIy5uWcSUFgpdeebRVkg3PI8HCgHYNlfCLiaAXEIIWxNKfDywtTLh8W1lfkDIDken9lNOHXtlK2jEymkO/lrrT/SWpfVWlcEugPBWutegAnoainWD1hleb7a8hrL/mCdE24qIIR4JCnXApry7gZmr4bIW5E0md2Ef67/Y+vwhEVmjvMfCnyglDqFuU1/tmX7bKCIZfsHwLBMjEEIkcVSrQXk7U3/Fu/z4xqIvHWNZ2c9y/kb561lZRSQ7cidvIQQmSsmBry8mFY7joHelyldoDQhr4dw/Mpx6zcE6QzOHHInLyGE7eTPD1On8lb79iR5t+W96PXUmV4HrTVLuy2VxG8jsryDECLztWsH3bvz7ldb6VPBjys3r+Ccx5mGpRvaOjK7JclfCJE1vvsOUzVn1p9cR/daAVyMuUizuc24lXjL1pHZJUn+QogsYbp5FP9uYPwliUWRrfik2SccunQIw3wDibcTbR2e3ZHkL4TIEiEXQjD2WoWhYgsYMoTRlV9n8DOD2RO+h5dXvUyyTrZ1iHZFOnyFEFki0CfQ/GRGOWjUCF54gYk7d1LcrTgfB39MTHwMK7uv5M5SX6YwEyEXQv57n8hQUvMXQmStqlVh7Vo4dw7atmVYnbcJqBXA6hOr6b28N4AsBpcFpOYvhMh6TZtCUBD4+aFeeolF69YRmxDLL4d/IeJmBH/8+4eM/89kUvMXQthGu3Ywbx6YTKiePVnZNYgaRWuw+fRmmpVvJok/k0nyF0LYTq9eMHkyrFzJ9kGdiIiNoEKhCqw4voLxu8bbOrpcTZK/EMK23n0X04i++BfYgPGqgT/fOEgl90oEbgnkx30/2jq6XEuSvxDC5kJa18SY6Idh3FIKfT+DHf13UMKtBB9u+pDQa6G2Di9XkuQvhLC5wKZDMUxYDgEBMHQopVcF89vLv6G1pvnc5lyKuWQtKyuBZgxJ/kKI7MHBAebPB4MB+ven2h9nGf/8eC7EXKDZ3GZEx0fLENAMJEs6CyGyl6goaN4cTp+G337j69iNfBz8MZXcK3Ej/oYMAX0MD1rSWWr+QojspVAhWL8ePDygXTs+KhNAh6c6cDryNOUKlqNFxRa2jjBXkOQvhMh+SpeGDRsgIQFTv+b87+xumpZryh///sHra163dXS5giR/IUT2VKMGpvmf498kHONvxdjuv55Wnq2Y/cdsPgn+xNbR5XiS/IUQ2VaIxy2MtUdh2HgC1akT615cRJ0SdRi7cyzbz2y3dXg5miR/IUS2FegTiKH3cPMyENu24dLqBbb5LsYjrwftF7bn+JXj1rIyBPTxSPIXQmR/ffqYVwI9eRJ3Qzu+qzOU2MRYDPMMXIq5JENAn4AM9RRC5Bz79kH79pCUxNSprzDw+DeULlCahKQEjN1kCOjdZKinECJ38PaG3bvBw4O3X/6BboWbcCH6AiXyl5AhoI8p3clfKVVOKWVSSh1VSh1RSg2ybPdQSm1WSp20/HS3bFdKqclKqVNKqUNKqQbpjUEIYUcqV4ZduzC1qIApfDctnatxJOII/Vf2t3VkOUpG1PyTgA+11jWBxsBApVRNYBiwVWtdFdhqeQ3QFqhqeQwApmVADEIIO2KKPYK/4QrGs0+z5eO/ecG1FvMPzeez4M9sHVqOke7kr7W+qLU+YHkeDRwDygB+wHxLsflAJ8tzP+AnbfY7UFgpVSq9cQgh7EfIhRCM3ZZimL8d1aoVqz87hpdreb7e9TV7w/faOrwcIUPb/JVSFYH6wB6ghNb6omXXv0AJy/MywLkUbztv2SaEEI8k0CfQ3Lnr4gIrV+JSvxHBYy5SwCEvvgt8ORf1X4qRIaBpy7Dkr5TKDywDBmutb6Tcp81Dih5rWJFSaoBSap9Sal9ERERGhSmEyG3y54dff6VY2WpMWBPP9bhIDPMNxCTEyBDQB8iQ5K+UcsKc+BdqrZdbNl+605xj+XnZsj0cKJfi7WUt21LRWs/QWntrrb2LFSuWEWEKIXIrDw/YuJH+EWX5epcroZGhNJzREP8gf1kF9D4yYrSPAmYDx7TW36bYtRroZ3neD1iVYntfy6ifxkBUiuYhIYR4MqVLw+bNDD3sTptzzpy4eoLqRapL4r+PjKj5+wB9gJZKqYOWRztgDPC8Uuok0NryGuBX4DRwCpgJvJ0BMQghBFSqhGnBaPYXSaTeNWd2ntvJCNMIW0eVLTmm9wBa652Aus/uVmmU18DA9J5XCCHuZgoz4b9/KMZnv6VJ1/dpMLwYo7aPolzBcrzW8DVbh5etyAxfIUSuEXIhxNzG32kwLn6d2TItlqKuHgzZMoQrN6/YOrxsRZK/ECLXsA4BBfjqK0pdiefXf1tzK/EW3ZZ2I/F2om0DzEYk+Qshcqdq1eC112g0ZQUznxnNtn+20W1pt1RF7HkOgCR/IUTuNWIEODnRZ85+utXsxqq/V/Hhxg8B7H4OQLo7fIUQItsqVQo++ABGj+aXD34n7HoY3/7+LdfirrH2xFq7ngMgNX8hRO42ZAgUKYLjx5+yqfcm3F3dmXdwHn3q9LHbxA+S/IUQuV3BgjB8OGzZwsF1swBwVI58v/d7NodutnFwtiPJXwiR+735JqbGJfA/8DHLui1l/kvzSUpOwm+xH6Ywk62jswlJ/kKI3M/FhZBeBoyLkjD8fometXsy+JnB3Eq6xYwDM2wdnU3IPXyFEPYhORkaNoTr1+H4cRIdHWj9c2tCwkP4/bXfqVOijq0jzHByD18hhHBwgHHj4J9/YPBgnBwcWdJ1CY4Ojvgu8CXyVqS1qD2M/5fkL4SwH88/D0OHwvTpMHEiJfOX5KtWX3Ex5iJtF7YlWSfbzfh/GecvhLAvX30FoaHwf/8Hnp6889I7nLp2ikl7JmGYZ+DolaN2Mf5fav5CCPvi4AA//QRPPw29ekFICBPbTKROiTpsP7udNpXb5PrED5L8hRD2KG9eWL0aSpSADh3Ytmcx4TfCKZavGL/89QvGI0ZbR5jpJPkLIexT8eKwbh2mojH4r+7D0nZz2dF/By6OLvRa3ivXTwCT5C+EsF81axIypCdGo8bwwfdUK1SJeX7zSEpO4ovtX9g6ukwlyV8IYdcC+83A8PFM2LwZBg8mwCuAgY0GsuPsDlYcW2Hr8DKNJH8hhHjlFfPon6lTYc4cJrwwgbIFytJ7RW9Cr4Vai+Wm8f+S/IUQAuDrr6F1a3jrLVwO/Mm458dxK/EWbRa04VbirVw3/l+WdxBCiDuuXgVvb0hMhP37+frEHD4O/pgGJRtw9sbZHDf+X5Z3EEKIR1GkCKxYAdeuQbdufNT4/2hStgkH/j1A8/LNc1TifxhJ/kIIkVK9ejBrFuzYgWloACeunaBswbKsOL6CBYcW2Dq6DCPJXwgh7tazJ6YhXfF3XIGxwCvsemUXbk5u9F/Vn42nNto6ugxhs+SvlPJVSv2tlDqllBpmqziEECItIR0aYjxVH8PgSZQ/cYkFnRfkqvH/Nkn+Sqk8wBSgLVAT6KGUqmmLWIQQIi2BzYZhmLHJvATE66/jV60j7z39HrvO7WLl8ZW2Di/dbFXzfxo4pbU+rbVOABYDfjaKRQgh0la0KHz2Gfz5J2zbxrjnx1G6QGn6LO/D2aiz1mI5cfy/rZJ/GeBcitfnLduslFIDlFL7lFL7IiIisjQ4IYSw6tULihWDiRNxcXRhbOuxxCbG0m5hOxJvJ+bY8f/ZtsNXaz1Da+2ttfYuVqyYrcMRQtgrV1d46y1YswZOnKB3nd580uwTjkQcocW8FvgH+ee48f9gu+QfDpRL8bqsZZsQQmQ/b78Nzs4waRIAX7T8gvol67P7/O4cu/6/rZJ/CFBVKeWplHIGugOrbRSLEEI8WIkS5uafefPg2jVMYSbORp2lSN4iLPprUY7sALZJ8tdaJwHvABuBY4BRa33EFrEIIcQjef99uHkT0/Sh+Af5s7TbUjb32YyDgwMBQQEEnw62dYSPxWZt/lrrX7XWT2mtK2utv7RVHEII8Uhq14bWrQnZuQRjp4UYPA3UL1Wfca3HkXA7ge/3fm/rCB9Ltu3wFUKIbOf99wlcH41h738jEAc1HoRvFV82hG7gyOWc04AhyV8IIR6Vry9Urw4TJ4JlRWQH5cA8v3kUdClIj2U9iEuKs3GQj0aSvxBCPCoHBxg8GPbvh507rZtL5C9B+6rt+evyXwzdPNS6PTtP/pLkL4QQj6NPH/DwgG+/Tb25Th/yOuZl8t7JrDuxLttP/pLkL4QQjyNfPvOkr1WrIPS/WzwaPA2sCFhBHpUH/6X+dF3aNVtP/pLkL4QQj2vgQHB0hMmTU21uU6UNrzd8nZtJN/HI60GLii1sE98jkOQvhBCPq1Qp6NkTpk2DtWutm01hJoKOBtGmchtOXTvFBxs/sGGQDybJXwghnsTEiVC3LnTuDKtXW9v4jV2N/NrrV7xLe/Pdnu+Yf3C+rSNNkyR/IYR4Eu7usHkz1K8PXboQsmmutY3fQTmwqvsqCjgX4JPgT0i8nWjraO8hyV8IIZ5U4cKwaRN4exP4ziIMB65Zd5UuUJo5fnMIjw7Plnf/kuQvhBDpUagQbNwITz8NAQFgNFp3da3Zlb51+/Llji/537n/2TDIe0nyF0KI9CpYEDZsgGefNXcEL15s3fV92+8p5FKIrsauxCTEWLfbegKYJH8hhMgIBQrA+vXg42Ne/vn33wEo6FKQEc+N4ELMBQKWBgBkiwlgSlvWp8jOvL299b59+2wdhhBCPNyNG+DlZf4wOHAAXFwA6BnUk0VHFuFfy5/gsOAsmQCmlNqvtfZOa5/U/IUQIiMVLAjTp8PRo/D119bN816aRwm3EhiPGOlXt5/NZ/5K8hdCiIzWrp257f+rr+CIeZnnXWd3kXA7AQcc+H7v9za/+YskfyGEyAzffWf+FvDqq5hCt+Af5M8y/2WMbjmahNsJdFrSCVOYyWbhSfIXQojMUKyY+Ybve/YQEjTZ2sY/xGcIjcs2RmvN1rCtNgtPkr8QQmSWnj2hbVsCR23FoDwBcHRwZH6n+dzWtzn470FsNehGkr8QQmQWpcydvw4O8MYb1rt/PVXkKca0HsO6k+uYe3CuTUKT5C+EEJmpfHnzqJ9Nm+Dnn62b33n6HSq5V+KdX9/hzPUz1u1ZNflLkr8QQmS2t982T/4aPBguXQLM9/4dbRhNXFIcLy15iWSdnKWTvyT5CyFEZnNwgJkzISrK3Als0aN2Dz589kP++PcP2i9sb10SOivmAEjyF0KIrFCjBnToALNmQXy8dfO458dR2b0yG0I30N2re5ZN/kpX8ldKjVdKHVdKHVJKrVBKFU6x7yOl1Cml1N9KqTYptvtatp1SSg1Lz/mFECJHefttiIiAoCDrpm3/bCMyLhLnPM78uO9Htp7OmuGf6a35bwa8tNZ1gBPARwBKqZpAd6AW4AtMVUrlUUrlAaYAbYGaQA9LWSGEyP1at4aqVWHKFOC/Bd6CugUxtd1UEpMTs2zyV7qSv9Z6k9Y6yfLyd6Cs5bkfsFhrHa+1DgNOAU9bHqe01qe11gnAYktZIYTI/RwczLX///0P/viDkAsh1jb+V+q/wguVXyDpdhIbQzdmfigZeKxXgPWW52WAcyn2nbdsu9/2eyilBiil9iml9kVERGRgmEIIYUP9+kHevDB1KoE+gdY2fqUUMzvMxCmPEyEXQjJ98tdDk79SaotS6nAaD78UZT4BkoCFGRWY1nqG1tpba+1drFixjDqsEELYlru7eb3/hQshMjLVrvKFyvPNC98QHBbMjP0zMjWMhyZ/rXVrrbVXGo9VAEqpl4EXgV76v4+qcKBcisOUtWy733YhhLAfAwfCrVswf/49u15v8DpVPKoweOPgTJ38ld7RPr5AINBRa30zxa7VQHellItSyhOoCuwFQoCqSilPpZQz5k7h1emJQQghcpx69cy3fJw6FZKTU+1SSjGqxSjikuLoYuyC1jpTJn+lt83/B6AAsFkpdVApNR1Aa30EMAJHgQ3AQK31bUvn8DvARuAYYLSUFUII+zJwIJw8CVu23LOrR+0eDHpmEPsv7qfjoo6ZMvlLbuMohBC2EB8P5cpBkyawcuU9u5N1MpUmVeJM1Bk+bfYpX7T84rFPIbdxFEKI7MbFBV57DdasgbNn79n92z+/EZ0QzVCfoUzfPz3Dx/5L8hdCCFt5803zz+nTU21OOflrTOsxGLsa8Q/yz9APAEn+QghhK+XLp7neT8rJXwAGTwPGrkZCLoRk2KmlzV8IIWxp0yZo08Z8z99BgzL00NLmL4QQ2VXr1uZO38GDzd8CwsKy5LSS/IUQwpYcHGDbNhg/HkwmqFkTRo2CuLjMPW2mHl0IIcTDOTnB//0fHD8OHTvCiBHg5QUbNmTaKSX5CyFEdlG2LCxZAps3g6MjtG0LAQH3zALOCI4ZfkQhhBDp07o1/PknTJwIMTHmpqEMJslfCCGyIxcXGJZ5NzuUZh8hhLBDkvyFEMIOSfIXQgg7JMlfCCHskCR/IYSwQ5L8hRDCDknyF0IIOyTJXwgh7FCOWNJZKRUBnHlowfsrClzJoHByCnu7Znu7XpBrthfpueYKWutiae3IEck/vZRS++63pnVuZW/XbG/XC3LN9iKzrlmafYQQwg5J8hdCCDtkL8l/hq0DsAF7u2Z7u16Qa7YXmXLNdtHmL4QQIjV7qfkLIYRIQZK/EELYoVyd/JVSvkqpv5VSp5RSmXdXhGxCKVVOKWVSSh1VSh1RSg2ydUxZRSmVRyn1h1Jqra1jyQpKqcJKqSCl1HGl1DGl1LO2jimzKaXet/y7PqyUWqSUcrV1TBlNKTVHKXVZKXU4xTYPpdRmpdRJy0/3jDhXrk3+Sqk8wBSgLVAT6KGUqmnbqDJdEvCh1rom0BgYaAfXfMcg4Jitg8hCk4ANWuvqQF1y+bUrpcoA7wHeWmsvIA/Q3bZRZYp5gO9d24YBW7XWVYGtltfplmuTP/A0cEprfVprnQAsBvxsHFOm0lpf1FofsDyPxpwQytg2qsynlCoLtAdm2TqWrKCUKgQ0B2YDaK0TtNbXbRpU1nAE8iqlHIF8wAUbx5PhtNbbgWt3bfYD5luezwc6ZcS5cnPyLwOcS/H6PHaQCO9QSlUE6gN7bBxKVvgOCASSbRxHVvEEIoC5lqauWUopN1sHlZm01uHAN8BZ4CIQpbXeZNuoskwJrfVFy/N/gRIZcdDcnPztllIqP7AMGKy1vmHreDKTUupF4LLWer+tY8lCjkADYJrWuj4QSwY1BWRXlnZuP8wffKUBN6VUb9tGlfW0eWx+hozPz83JPxwol+J1Wcu2XE0p5YQ58S/UWi+3dTxZwAfoqJT6B3PTXkul1ALbhpTpzgPntdZ3vtUFYf4wyM1aA2Fa6witdSKwHGhi45iyyiWlVCkAy8/LGXHQ3Jz8Q4CqSilPpZQz5s6h1TaOKVMppRTmduBjWutvbR1PVtBaf6S1Lqu1roj5bxystc7VNUKt9b/AOaVUNcumVsBRG4aUFc4CjZVS+Sz/zluRyzu5U1gN9LM87wesyoiDOmbEQbIjrXWSUuodYCPmkQFztNZHbBxWZvMB+gB/KaUOWrZ9rLX+1XYhiUzyLrDQUrE5DfS3cTyZSmu9RykVBBzAPKrtD3LhUg9KqUVAC6CoUuo8MAIYAxiVUq9iXtreP0POJcs7CCGE/cnNzT5CCCHuQ5K/EELYIUn+QghhhyT5CyGEHZLkL4QQdkiSvxBC2CFJ/kIIYYf+H3bct8WkOb6+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "n_epoch = 3000 # epoch size\n", "a, b, c = 1.0, 1.0, 1.0 # initial parameters\n", "epsilon = 0.0001 # learning rate\n", "\n", "N = np.size(t)\n", "\n", "for i in range(n_epoch):\n", " for j in range(N):\n", " a = a + epsilon*2*(y[j] - a*t[j]**2 - b*t[j] - c)*t[j]**2\n", " b = b + epsilon*2*(y[j] - a*t[j]**2 - b*t[j] - c)*t[j]\n", " c = c + epsilon*2*(y[j] - a*t[j]**2 - b*t[j] - c)\n", "\n", " L = 0\n", " for j in range(N):\n", " L = L + (y[j] - a*t[j]**2 - b*t[j] - c)**2\n", " \n", " if i % 500 == 0:\n", " print(\"epoch %4d: loss = %10g, a = %10g, b = %10g, c = %10g\" % (i, L, a, b, c))\n", " \n", " \n", "y_est = a*t**2 + b*t + c \n", "\n", "\n", "plt.plot(t, y, 'r-', label='Real data')\n", "plt.plot(t, y_est, 'g-x', label='Estimated data')\n", "plt.legend()\n", "plt.savefig(\"fig-res-missle_est.pdf\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. 如何使用sklearn求解线性问题?\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X: (100, 1)\n", "Y: (100, 1)\n", "a = 2.773500, b = 5.311307\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnF0lEQVR4nO3deZhU9Z3v8feXpoHGhQZhRBoIoIgSUVpbI+LGEkHMKJc4xmzXyZhwZyaLJiMTTGYm88yjA9H7JM69k0kuT3Suc8fHZaJBJ6IoAm6JTppFkRCQTaFFJSooi7J97x+n2mq6lq6uOnXOqarP63l8pE51Vf3O6a7f97d+j7k7IiIiHfWIuwAiIpI8Cg4iIpJBwUFERDIoOIiISAYFBxERydAz7gKEYeDAgT5ixIi4iyEiUlFWrFjxB3cflO25qggOI0aMoLW1Ne5iiIhUFDN7LddzGlYSEZEMCg4iIpJBwUFERDIoOIiISAYFBxERyVAVq5VERHJZuKqN2xev541d+xnS2MCcaWOY2dwUd7EST8FBRKrWwlVt3PzQGvYfPAxA26793PzQGgAFiC5oWElEqtbti9d/HBja7T94mNsXr4+pRJVDwUFEqtYbu/Z367ikKTiISNUa0tjQreOSpuAgIlVrzrQxNNTXHXWsob6OOdPGxFSiyqEJaRGpWu2Tzlqt1H0KDiJS1WY2NykYFEHDSiIikkE9BxGRhIpzA5+Cg4hIEYqtuAt9Xdwb+GIfVjKzOjNbZWa/Sj0eaWYvmtlGM7vfzHrFXUYRkY7aK+62Xftx0hX3wlVtob0u7g18sQcH4AZgXYfHPwR+7O6nAO8B18dSKhGRHIqtuLvzurg38MUaHMxsKHAF8PPUYwMmA79I/cjdwMxYCicikkOxFXd3Xhf3Br64ew53AH8NHEk9PgHY5e6HUo+3A1qDJiKJUmzF3Z3Xxb2BL7bgYGafAd529xVFvn62mbWaWevOnTtDLp2ISG7FVtzded3M5ibmzRpHU2MDBjQ1NjBv1riaWK00EbjSzGYAfYDjgX8CGs2sZ6r3MBTIOsPj7guABQAtLS0eTZFFRIrfed3d18W5gc/c469XzexS4CZ3/4yZ/QfwoLvfZ2Y/A15293/J9/qWlhZvbW2NoKQiItXDzFa4e0u255K4z+G7wH1mdguwCrgz5vKISBXTneKyS0RwcPflwPLUvzcD58VZHhGpDWFtNKvGABP3aiURkdiEsdGs2A1xSafgICI1K4yNZnHvZC6XRAwriUj5VePQR6mGNDbQ1o0NaNnEvZO5XNRzEKkB1Tr0UaowNprFvZO5XBQcRGpAtQ59lCqMjWZx72QuFw0ridSAah36CEOpG82q9VakCg4iNSCMsXXJLeqdzFHMH2lYSSRCC1e1MXH+UkbOfZSJ85dGNuZfrUMftSiq+SMFB5GIxDkpHHcSNwlPVPNHGlYSiUi+L3UUlXScSdwkPFHNH6nnIBIRTQpLGKJaOqvgIBKRal0PL9GKav5IwUEkIpoUljBENX+kOQeRiFTreniJXhTzRwoOIhHSpLBUCg0riYhUsjLdzVM9BxGREkWa8fbIEWhthUWLgv+++EW44YbQP0bBQaQM4k6PHffn15Kw7iaX17vvwuLF8Nhj8PjjsHMnmMH558OgQeF8RicKDiIhi6SyKPHzFTzCU5bNje6wenW6d/DCC0GP4YQTYPp0mDEDLrsMBg4s/QRyUHAQCVncO6G7+vy4g1e1yZbQEIrY3Lh7NyxZEgSDxx6DHTuC4y0t8Dd/A5dfDueeC3V1+d8nJAoOIiGLeyd0V58fd/CqJgtXtWFAtinhLjc3usPatUEgWLQInnsODh2Cfv1g2rSgdzB9Opx4YjmK3iUFB5GQxZ0eu6vPjzt4VZPbF6/PGhgMsm9u3LMHli5NDxdt2xYcP/NMuOmmICBMmAA946+atZRVJGRx74Tu6vOVxiM8uQKqkxqic4cNG+COO4I5ghNOgKuugnvuCYaLFiwIAsRLL8G8eXDRRYkIDBBjz8HM+gDPAL1T5fiFu//AzEYC9wEnACuAL7v7gbjKKdJdce+E7urz50wbc9ScAyiNR7Gy9dJ6H/yIK97dAN9aHPQONm0Knjj9dPjmN4PewYUXQq9eMZS4cOZl2kDR5QebGXCMu+8xs3rgOeAG4DvAQ+5+n5n9DHjJ3X+a771aWlq8tbW1/IUWqRJarRSO9sn9E3a2MWlzK5M2tTLh9TU0HPoIGhpg8uQgGFx+OYwcGXdxM5jZCndvyfpcXMHhqEKY9SUIDn8BPAoMdvdDZjYB+Ht3n5bv9QoOIhKpAwfg2Wdh0SI+eOgRjtu6EYDtA4bw0bTpnHzdNXDxxUGASLB8wSHWwS0zqyMYOjoF+AmwCdjl7odSP7IdyNqcMbPZwGyA4cOHl7+wIt2k1nnxEnnttm9PryxasiSYXO7Vi+MuvRRu/AbMmMHQ0aPjLWOIYg0O7n4YGG9mjcAvgdO68doFwAIIeg5lKaBIkSphL0EiK2BKu3ahntPBg/Cb36RXFq0JysDw4fClLwXDRZMnwzHHFPf+CZeIaXF332Vmy4AJQKOZ9Uz1HoYC0dyBXSSHYiqcpO8lSHLwKvbahXJOb74ZpKdYtAieeCLYmNazZzCBfNttQUAYOzZIXVHl4lytNAg4mAoMDcCngR8Cy4CrCVYsXQc8HFcZRYqtcOLaS1BoIEty8Cr22hV1TocPw29/m+4drFgRHD/pJLj66iAYTJ0Kxx/f7fOodHH2HE4C7k7NO/QAHnD3X5nZ74D7zOwWYBVwZ4xllBpXbCUax0a47gSyXBVtrlQQUSr22hUcVP7wh6OT2L3zDvToEWw+u/XWICCcdVZN9A7yiS04uPvLQHOW45uB86IvkUimYluxcewl6E4gy1UBG0GQibP3UOy1y3VOTcf3DnoE7b2DF18MNqcNHBgEgvYkdgMGhH4uhUjq3E8i5hxEkqrYVmwcG+G6E8jmTBvDt+9fnZH6wSH2oaVir13HoHL8h3u4cOtqpm5ZwYy21fD9VIrrc8+FH/wg2HfQ0hL0GGKU5LkfBQeRPErpAUR9S9DuBLKZzU3ceP/qrO+ThBxL3b527syse4fR7z/Dgf98lHFbX6GnH+HA8f3oNePyoHcwbRr80R+Vr9BFSPLcj4KDSB5xp8Loju4GsqaYEwSW7IMP4Kmn0sNFbW18EqC5GW6eCzNm0OtTn0pMrqJskpwEMblXTSQhou4BdEfn8erPntPEst/vLCiQVVyOJXdYvx4WLeLt+39J/xUvUH/4EHt792XXhZfS9A//EKS4HjIk7pIWLO4MvvkoOEhVS+pkXxiyjVc/uKKNebPGFXSOFdEr2rcPli9P9w62bAFg96BP8NA5V7L85BZWNJ1Ozz59mNc8jpkVFBgg2QE6EbmVSqXcSpJN58oTgi9eoZVnOcsVRoU8cf7S7KtzGht4fu7kMIpasFCD8KZN6WCwfDl8+CH07QtTpsCMGcza0o+VlrnvII7zDkOcDZjE5laSylJprfAkTvaFuTolKePVJZ/TRx/BM8+kA8KGDcHxU0+FP//zYDL5oougTx8AVs19NOvbJGGcvhhJHbZUcJCCJHnJXS5JqTw7CjNgxT1e3d5YyFaGLs/p9dfTSeyeegr27oXevWHSJPjGN4KlpqeckvWlYZ13pTV2oqbgIAVJYiu8K3FWnrkqnjADVpzj1dmG7Do76pwOHoTnn0/3DtauDY6PGAHXXRf0DiZNCoaPuhDGeVdiYydqCg5SkCS2wrsSV+WZr+IJM2DFOaGcrbHQ2Zk99sJddwXB4Mkn4f33ob4+uM/BV74SBITTTsubpiJf676U867Exk7UFBykIHEPYRQjrsozX8UTdsCKcry6Y0WdbRlLjyOHGf/GBiZtbmXKllbGvpm6PWZTE3zuc0EwmDIFjjuu4M/L17ov5bwrsbETNQUHKUiSl9zlE8dkX76KpyKWj2aRaxhpwL7dXLJ5BZM2t3LRllX0//ADDlkPtp8+Hm6YFwSEceOKSmJXztZ9JTZ2Oiv3nImCgxSkUiu1OHRV8SR1dUo+7RW1+RHGvbmRSZtambS5lTN3vEoPnJ3HNLLklE+xfNQ5PDuymeMGDyp5WWk5WvcdJ9ENjuoBVUJjp10UcyYKDlKwSqzU4lCpvayc3nuP5t8s5jubW7lk8woG7tvNEYzVQ07lxxd+gWUnn8vaE0fhlk5i90EIwzNht+47V6gOHweIpgpr7EQxZ6LgIBKyiu9lucNLL6WXmv761/zzkSO81+c4nh51NstGtfDMyLN5r28/6sw4nGUjbRjDM2EH2WwVantgqLTNc1HMmSg4iJRBxfWy3n8fliwJgsFjj8EbbwTHzzkHvvc9nj75XP7y1Tr2dqhbG+rr+Ow5TTy4oq0svaSwg2w1TUJHMWei4CBSi9xh3br0voNnn4VDh4LbYU6bFkwkT58OgwcDcAlwa44J0JZPDChbLynMIFsNk9Dtohi6VG4lkVqxdy8sW5YOCK+9FhwfNy59R7QJE4K9CFUoqbm2ihXGaiXlVhLphqpKq/Dqq+lg8PTTQR6jY46BqVPhe98L0lQMG5bxsqq6BikVPxfUSbmHLtVzEOmg4luXH34YBIH2gLBxY3D8tNPSvYMLLwzyGOVQ8ddACqaeQ4WoxtZaOYR1nbK9T0WmVdi69egkdvv3BxlMJ0+GG28MegejRhX8dhV5DSR0Cg4JoURghQnrOuV6n1z5ghK1ouXAAXjuuXTvYN264PioUXD99UHv4NJLoaG4idZqWtUjxVNwSIiktdaS2osJ6zrlep9yrtsvSVtbunfw5JOwZw/06gWXXAKzZwe9g1NPLSpNRWf9GurZtf9gxvHYr4FEKrbgYGbDgH8DTiTYi7LA3f/JzAYA9wMjgK3ANe7+XlzljEqSWmtJ7sWEdZ1y/fxhdxrq60JZIlhSgD10CF54Id07eOml4PiwYfDFLwa9g8mT4dhju12ursq898ChjOP1Paxyd3hLUXp0/SNlcwj4K3cfC5wPfN3MxgJzgafcfTTwVOpx1cvVKoujtZavdR63XNfDCW6buXBVW0nv09hQz7xZ42hqbMAIds8WMxHbHmDbUhlM2wNs3vK99RbcfXeQwXTQoODuZ7fdBo2N8MMfwpo1wfLTn/0Mrrwy9MAAwe/+4OHMntOxfXrG3jCQaMXWc3D3HcCO1L8/MLN1QBNwFXBp6sfuBpYD342hiAUptnXY+XWTThtUtp2m3SnHnGljEtWL6Szb5p923enhzJk2hjn/8RIHjxxdEba3mktNp1DQ8Nfhw9Damu4dtK+4GzwYZs0Khoo+/Wno16+ksnRHrt/xrn2Zw0yVKKnDpUkUZ8/hY2Y2AmgGXgROTAUOgDcJhp2yvWa2mbWaWevOnTujKWgnRbUOc7zuwRVtfPacppJbrGGUv19D9k1QSRhzntnc9HHLPptCezgzm5s4tk9m2+jgYQ+lh5Srkt234y2491748peDIHD++XDLLcH8wS23wMqVwfzCnXfC1VdHGhggWT3YsBX7fa1VsU9Im9mxwIPAje7+vnWYUHN3N7OsGzHcfQGwAIJ9DlGUtbNiJ0dzvW7Z73dGmgAsVzn61PcIbdy9HNo3/4yc+2jWm84U2sPJ1RoOo4fUnqrB/Ahj39rMpM2tXLp5Bc1vrIf/dQQGDgzSU8yYAZddBiecUPJnhqHqMsp2UMpihlrsccQaHMysniAw3OPuD6UOv2VmJ7n7DjM7CXg7vhLmV+zwS1KGbfINIfz4c+MT/2UoNVdO2XLt7N7Nj+teZdvjD3LRxlb+aG+wnuLlIafy6le/xWnXXwstLVBXV9rnFKGrSq7adhF3VOz3LskLNMopztVKBtwJrHP3H3V46hHgOmB+6v8Px1C8ghRbuSQlAVi+clRCVtFSW7mhtZLdYe3a9NzB889z3qFDjD+uH8+c3Mxjw5pZf9YFfPXqCbFe00IruUr43Rej2O9d0paZRyXOnsNE4MvAGjNbnTr2PYKg8ICZXQ+8BlwTT/GOlq3FVWzlkpSue1LKUawwWrl96nt8fP6NDfX8/ZWfLOz1e/YEu5Hb9x5s2xYcP+ssmDMHZsyg1/nnM7VnT6Z2+8zKo1YruXbF/r0npacftThXKz1HcCOmbKZEWZau5GpxzZs1jnmzxnW7ckpK1z0p5ShFsa3cbPmDPjp0JPcL3GHDhnTv4Jlngp3Kxx0XJLH7u78LVhc1Jffa1Wol167Yv/ek9PSjpsR7BZg4f2nWP45KvIOUBAr6ne7fD8uXpwPC5s3B8bFj00nsJk4MVhpVAP0dF6eaExEq8V6JarnFVa2rNHL97uq2boGf/CQIBkuXBllOGxpgyhS46aagdzBiRLSFDcmcaWOY84uXjtrkVl+nnc9dKbTHUW3fFQWHAtRqt7KaV2m0/057HTrIudvXMmnTb5m0eQUnv7s9+IFTTglyFs2YEeQv6tMn3gKHpfNAQeUPHBSluxV5V8OX1fhdUXAoQKVP3BYryROYJbXStm3jn/f8lnd++TATtqzimIMf8lFdPf/1iTPZe/3XOPNrn4fRo8t7AjG4ffH6jB3hB494In6fUSpHRZ7k70qxFBwKUA0Tt8VI6nBat7/cBw/Cr3+dXlm0Zg3NwL7BTTx+9mU8OnQ8W888j2/+8XguSsjvtBxDFEn9fUatHBV5NV5bBYcCVeva73ySOpxW0Jd7xw54/PEgGDzxBLz/PvTsGSSzu/12mDGDvqefziwzZkVY9kIq/XINUST19xm1clTk1XhtFRwkQ3sF1rZrP8bRw9JJGE7L9iXuceQwJ65dCX+7LAgIK1cGTwwZAtdcE0wkT50Kxx9f8OeE3XovtNIv1xBFrQ6PdlaOirwar62CgxylcwXm8HGAaErIcFr7l7v/vt1cvGUlkza3csnmlfT/8APo0QMuuAD+8R+DyeQzzyzqBjhh3nGuPcD0yHIjoWyVfrmGKGp1eLSzclTk1XhtFRzkKNlare2BIfa18EeOwMqV/Py1RRz4z18xrm0DPXD+0LcfT48+j6YvfJZz/8e10L9/yR8VRuu9c4DJdoc5yKz0yzlEUYvDo52VqyKvtmur4CBHSdzE2nvvBbfFXLQomFB++21ON+PdT47nrtP+O//ZdBbvjBnHTZefzrkhfjHDuA7ZAkw2nSv9crRsq20NfqmqrSIvBwWHKhDmFz/2iTX34I5n7buSf/3r4KY4/funU1xPm8aAQYP4KvDVMhUjjOtQSCDJVumH3bKtxjX4Un4KDhWulC9+mMkES/LBB7BkSbp30Ja6+UpzM8ydGwSE884LVht1IaxAGcZ1yBVg6sw44p63fGG2bKtxDb6UX83mVqqWbnax+XLy5YuBMk+sucPvf5/uHTz7bLAX4fjjg9tizpgR9BKGDOnW24adAyfbrVyX/X5nwdclKTl5ct0UyYAt86+IrBySPPlyK9VkcEjKlzYMub74EASIXBVZ5EnY9u2DZcvSAWHr1uD4GWekk9hdcAHUZ79FaSHKeU7F/s0koRGihHuSixLvdVJN3excQxcGHx/PNtQUycTzxo3pYLB8OXz0EfTtG+w3mDs32HswfHhoH1fOcyr2byYJE5/VuAZfyq8mg0OSVuQU27LMt1Gt82PIrMjKMvH84YfBfQ7aA8KrrwbHTz0V/uIvgt7BxRdD797Ff0Ye5ZxMT9LfTHdV4xp8Kb+cwcHMFgF/6e5boytONGJfkZNS7GRyVxvVsp0bHF2RhdaafO21dM6ip54Kho/69IFJk+Bb3wp6Byef3L33LFI5W8hJ+ZspVrYeTBKGvCS5euR57l+BJ8zs+2ZW/EBwAs2ZNoaG+qNv7h5HNzvfUEV3X9dxo1pTjgqrY0U2s7mJebPG0dTYgKVeW9Ccy4EDwdzBnDnBfMGIEUGvYM0a+MpX4NFH4Z13gmDxjW9EFhgg85z6962nd88efPv+1Uycv5SFq9qKfu+k/M2Epb2B0bZrP066YVLKNZLqkrPn4O7/YWaPAX8LtJrZ/wOOdHj+RxGUryyS0s0udqiiq9cV2oIueDz8jTfSvYMnnwyWntbXB0NEf/ZnwXDRmDFFpakIW/s5hb22Pyl/M2Gppnm3ztQjCkdXcw4HgL1Ab+A4OgSHSpeEicJihyq6el3JFdmhQ/Dii+m5g9Wrg+NDh8LnPx8MFU2ZEtw/OaHKUfnl+pupxMqokudQ8tGGv/Dkm3OYDvwIeAQ42933RVaqGlHsGHkhr+t28Hv7bVi8OAgGixcHaSvq6oJ7JM+fH/QOzjij6CR2UVeeUVV+lVoZVfocSi7V3COKWr6ew/eBP3H3tVEVptYU28IPZYjjyBFobU33Dlpbg81pJ54IV10VBINPfxoaG4s9PSC+yjOqyq9SK6NqXd5arT2iOOSbc7goyoLUqmKHt4p63bvvBje+WbQouBHOzp24Ga8MO50nJ36BV8ZfyJV/egUzzxnW7fLkElflma/yC7MnU6mVUbXNobSr1h5RHGLd52BmdwGfAd529zNSxwYA9wMjgK3ANe7+XlxlrGju8NJL6d7Bb34T9BhOOAGmT6d17Pl8851B7Kg/9uOX/GbhWujRI5RKYuGqtoKW1ZZDrsoPKHr5cLaKtJIroyTMu4WtWntEcYg1fYaZXQzsAf6tQ3C4DXjX3eeb2Vygv7t/N9/7FJNbqWrt3n10ErsdO4Lj55yTTlNx7rlQVxd5uolcnxHlnEQx59xVHqpqScVSLSpxgUBcEps+w92fMbMRnQ5fBVya+vfdwHIgb3Coae7wu9+lewfPPResNurXDy67LJ3EbvDgjJdGnW6iXceWXNRzEsWcc76hsfaAosooOaqxRxSHJKbPONHdU81d3gROzPZDZjYbmA0wPMT8PBVh715YujQdEF5/PTh+5plw001BQDj//C6T2MWRbgKgT31672XUcxLFnHNXAUWVkVSjfDukY+fBmFfWcS93X+DuLe7eMmjQoIhLFjF32LAB7rgj6A0MGABXXgn//u9w9tmwYAFs2xbML8ybBxddVFB203Lu+s1X2b637+DHu3GjntAt5pxznUslzCuIFCuJPYe3zOwkd99hZicBb8ddoFjs3w9PP53uHWzaFBw//fQgLcWMGXDhhSUlsSvnipVsE4MdtfcOop7QLeacNckptSiJweER4Dpgfur/D8dbnAht2ZJOU7F0aRAgGhpg8mT4zneCnckjR4b6keUaEulYCedbsfTjz42PvOLt7jlX67JPkXziXq10L8Hk80DgLeAHwELgAWA48BrBUtZ3871Pxa5WOnAguAta+8qideuC46NGwRVXBL2DSy4JAkQF62qFkFaXiMQjyauVPp/jqSmRFiRK27enewdLlsCePdCrVxAEZs8OAsLo0YlIYheWroZlNKErkjxJHFaqLocOBZvP2ucOXn45OD5sGHzpS0EwmDQJjj02//tUMA3LiFQeBYdyePPNID3FY48FSex274aePYMJ5NtuCwLC2LFV1TvoinoHIpVFwaEA2cbEId0SHnp8L25p2s8lm1KJ7FasCF44eDB89rNBMJg6NdiYJl3SHIRI/GKdkA5LOSeks6VOqK8zGve+z4TNK5i0qZVLtqxkwP738R49sAkTglVFM2bA+PEF9Q5UGablS1VRTO4jEcktsRPSlaB9B6/5Eca+tZlJm1uZtKmV8Ts2UOdHeKfheJaNOoflo1p4dfwFPP4PM7v1/pV6P4ByKWbHdLZr+O37V9P62rvcMnNc2cssUo0UHPLZtYvxLzzJjZtauXRLK4P27uIIxssnncL/vuBzLBvVwssnjcYt2GhuB7r/EZV6P4ByCSv3kQP3vPA6LZ8YUJPXUaRUCg4ducMrr6RXFj3/PD85fJjdvY/h6VHnsGxUC8+MPJt3jmnM+vJidvVW6v0AyiXM3EcONRtkRUql4LBnDzz1VDogbN8eHB8/Hr77XZ455Vz+cmM9e7JngThKMbt6K/l+AOVQTKqKXNcQajfIipQq0Yn3yu6BB4IkdjNnwr33wnnnwc9/HgSIVavg1lu5+CszueXq8TQ1NmBAXY4J5v5964tqoZYz+V0lmtncxLxZ4z6+3k2NDV1ORs+ZNoZc0/61GmRFSlXbPYfmZrjhhmBl0cSJwU7lLDqu0c+1muYHf/zJooqgDWKZisl91Prau9zzwutHpfAtJMhqlZNIdlrKWgRVKMnU3d9LsctmRapFvqWsCg5Ss8p5m1SRSpAvONT2nIPUNK0UE8lNwUFqlu7wJpKbgoPULK0UE8mttlcrSU3TSjGR3BQcstBqpNqhVOIi2Sk4dKJEeCIimnPIkC8RnohIrVDPoYOFq9qUoyehNNQnEi0Fh5T24aRctLwxPhrqE4mehpVSsg0ntdPyxnhpqE8kegoOKfmGjZRrJ17aySwSvcQGBzObbmbrzWyjmc0t9+flGjZqamxQYOhk4ao2Js5fysi5jzJx/lIWrmor6+dpJ7NI9BIZHMysDvgJcDkwFvi8mY0t52dqt2xh2sf/23btx0mP/5czQOh3IxK9RAYH4Dxgo7tvdvcDwH3AVeX8wGJuMlOL4hj/1+9GJHpJXa3UBGzr8Hg78KmOP2Bms4HZAMOHDw/lQ7Vbtmtxjf/rdyMSraT2HLrk7gvcvcXdWwYNGhR3cWqGxv9FakNSg0MbMKzD46GpYxIzjf+L1IakDiv9FhhtZiMJgsK1wBfiLZKAMpmK1IrE3ibUzGYAdwB1wF3ufmuun9VtQuOjtBYilSvfbUKT2nPA3RcBi+Iuh+SmtBYi1Supcw5SAZTWQqR6KThI0ZTWQqR6KThI0bSsVaR6KThI0bSsVaR6JXZCWpJPy1pFqpeCg5REaS1EqpOGlUREJIN6DiXSJjARqUYKDiXQJjARqVYaViqBNoGJSLVSz6EESd8EpiEvESmWeg4lSPImsDhu5yki1UPBoQRJ3gSmIS8RKYWGlUqQ5E1gSR/yEpFkU3AoUVI3gQ1pbKAtSyBIwpCXiCSfhpWqVJKHvEQk+dRzqFJJHvISkeRTcKhiSR3yEpHk07CSiIhkUHAQEZEMCg4iIpJBwUFERDLEEhzM7E/MbK2ZHTGzlk7P3WxmG81svZlNi6N8IiK1Lq7VSq8As4D/0/GgmY0FrgU+CQwBlpjZqe5+OPMtRESkXGLpObj7OnfPluTnKuA+d//I3bcAG4Hzoi2diIgkbc6hCdjW4fH21DEREYlQ2YaVzGwJMDjLU99394dDeP/ZwGyA4cOHl/p2IiLSQdmCg7tPLeJlbcCwDo+Hpo5le/8FwAKAlpYWL+KzREQkh6QNKz0CXGtmvc1sJDAa+K+YyyQiUnPiWsr638xsOzABeNTMFgO4+1rgAeB3wOPA17VSSUQkerEsZXX3XwK/zPHcrcCt0ZZIREQ6StqwkoiIJICCg4iIZFBwEBGRDAoOIiKSQcFBREQyKDiIiEgGBQcREcmg4CAiIhkUHEREJIOCg4iIZFBwEBGRDAoOIiKSQcFBREQyKDiIiEgGBQcREcmg4CAiIhkUHEREJIOCg4iIZFBwEBGRDLHcQ7pSLFzVxu2L1/PGrv0MaWxgzrQxzGxuirtYIiJlp+CQw8JVbdz80Br2HzwMQNuu/dz80BoABQgRqXoaVsrh9sXrPw4M7fYfPMzti9fHVCIRkegoOOTwxq793TouIlJNYgkOZna7mf3ezF42s1+aWWOH5242s41mtt7MpsVRPoAhjQ3dOi4iUk3i6jk8CZzh7mcCG4CbAcxsLHAt8ElgOvAvZlYXRwHnTBtDQ/3RH91QX8ecaWPiKI6ISKRiCQ7u/oS7H0o9fAEYmvr3VcB97v6Ru28BNgLnxVHGmc1NzJs1jqbGBgxoamxg3qxxmowWkZqQhNVKfwbcn/p3E0GwaLc9dSyDmc0GZgMMHz68LAWb2dykYCAiNalswcHMlgCDszz1fXd/OPUz3wcOAfd09/3dfQGwAKClpcVLKKqIiHRStuDg7lPzPW9mfwp8Bpji7u2VexswrMOPDU0dExGRCMW1Wmk68NfAle6+r8NTjwDXmllvMxsJjAb+K44yiojUsrjmHP4Z6A08aWYAL7j7n7v7WjN7APgdwXDT1939cJ73ERGRMoglOLj7KXmeuxW4NcLiiIhIJ5Ye7q9cZrYTeK3Ilw8E/hBicSqBzrk21OI5Q22ed7Hn/Al3H5TtiaoIDqUws1Z3b4m7HFHSOdeGWjxnqM3zLsc5K7eSiIhkUHAQEZEMCg6pjXQ1RudcG2rxnKE2zzv0c675OQcREcmknoOIiGRQcBARkQw1HRzMbHrqpkIbzWxu3OUpNzMbZmbLzOx3ZrbWzG6Iu0xRMbM6M1tlZr+KuyxRMLNGM/tF6qZa68xsQtxlKjcz+3bq7/oVM7vXzPrEXaZyMLO7zOxtM3ulw7EBZvakmb2a+n//Uj+nZoND6iZCPwEuB8YCn0/dbKiaHQL+yt3HAucDX6+Bc253A7Au7kJE6J+Ax939NOAsqvzczawJ+BbQ4u5nAHUENw6rRv+X4GZoHc0FnnL30cBTqcclqdngQHAToY3uvtndDwD3EdxsqGq5+w53X5n69wcEFUbV37DCzIYCVwA/j7ssUTCzfsDFwJ0A7n7A3XfFWqho9AQazKwn0Bd4I+bylIW7PwO82+nwVcDdqX/fDcws9XNqOTg0Ads6PM55Y6FqZGYjgGbgxZiLEoU7CLIAH4m5HFEZCewE/jU1lPZzMzsm7kKVk7u3Af8TeB3YAex29yfiLVWkTnT3Hal/vwmcWOob1nJwqFlmdizwIHCju78fd3nKycw+A7zt7iviLkuEegJnAz9192ZgLyEMMyRZaoz9KoLAOAQ4xsy+FG+p4pG6P07JexRqOTjU5I2FzKyeIDDc4+4PxV2eCEwErjSzrQRDh5PN7N/jLVLZbQe2u3t7r/AXBMGimk0Ftrj7Tnc/CDwEXBBzmaL0lpmdBJD6/9ulvmEtB4ffAqPNbKSZ9SKYvHok5jKVlQU3z7gTWOfuP4q7PFFw95vdfai7jyD4HS9196puUbr7m8A2MxuTOjSF4B4p1ex14Hwz65v6O59ClU/Cd/IIcF3q39cBD5f6hnHd7Cd27n7IzL4BLCZY2XCXu6+NuVjlNhH4MrDGzFanjn3P3RfFVyQpk28C96QaPpuBr8RcnrJy9xfN7BfASoJVeauo0jQaZnYvcCkw0My2Az8A5gMPmNn1BLcvuKbkz1H6DBER6ayWh5VERCQHBQcREcmg4CAiIhkUHEREJIOCg4iIZFBwEAlZKvvtFjMbkHrcP/V4RMxFEymYgoNIyNx9G/BTgrXnpP6/wN23xlYokW7SPgeRMkilKVkB3AV8DRifSusgUhFqdoe0SDm5+0EzmwM8DlymwCCVRsNKIuVzOUH66DPiLohIdyk4iJSBmY0HPk1wx71vt2fMFKkUCg4iIUtlBf0pwf0yXgduJ7gRjUjFUHAQCd/XgNfd/cnU438BTjezS2Isk0i3aLWSiIhkUM9BREQyKDiIiEgGBQcREcmg4CAiIhkUHEREJIOCg4iIZFBwEBGRDP8fU4TJRX/2GGsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "from sklearn import linear_model\n", "import numpy as np\n", "\n", "# load data\n", "# generate data\n", "data_num = 100\n", "X = np.random.rand(data_num, 1)*10\n", "Y = X * 3 + 4 + 8*np.random.randn(data_num,1)\n", "\n", "print(\"X: \", X.shape)\n", "print(\"Y: \", Y.shape)\n", "\n", "# create regression model\n", "regr = linear_model.LinearRegression()\n", "regr.fit(X, Y)\n", "\n", "a, b = np.squeeze(regr.coef_), np.squeeze(regr.intercept_)\n", "\n", "print(\"a = %f, b = %f\" % (a, b))\n", "\n", "x_min = np.min(X)\n", "x_max = np.max(X)\n", "y_min = a * x_min + b\n", "y_max = a * x_max + b\n", "\n", "plt.scatter(X, Y)\n", "plt.plot([x_min, x_max], [y_min, y_max], 'r')\n", "plt.xlabel(\"X\")\n", "plt.ylabel(\"Y\")\n", "plt.savefig(\"fig-res-sklearn_linear_fitting.pdf\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. 如何使用sklearn拟合多项式函数?" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([800., 90., -20.])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Fitting polynomial functions\n", "\n", "from sklearn.preprocessing import PolynomialFeatures\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.pipeline import Pipeline\n", "\n", "t = np.array([2, 4, 6, 8])\n", "\n", "pa = -20\n", "pb = 90\n", "pc = 800\n", "\n", "y = pa*t**2 + pb*t + pc\n", "\n", "model = Pipeline([('poly', PolynomialFeatures(degree=2)),\n", " ('linear', LinearRegression(fit_intercept=False))])\n", "model = model.fit(t[:, np.newaxis], y)\n", "model.named_steps['linear'].coef_\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 参考资料\n", "* [梯度下降法](https://blog.csdn.net/u010402786/article/details/51188876)\n", "* [如何理解最小二乘法?](https://blog.csdn.net/ccnt_2012/article/details/81127117)\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.9" } }, "nbformat": 4, "nbformat_minor": 2 }