{ "cells": [ { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VfX9x/HXhySQsBLCDIEQZgKIEIiIo4qgogiOVqyjitYWu5c/qvhrKyiu4ihWq2LVqq2l/vw5whAUEFddDFGzIOwkQFgJkD2+vz9y6A8xQBKS3PV+Ph553HvPPefez5cT3vfkm885MeccIiISGlr5ugAREWk5Cn0RkRCi0BcRCSEKfRGREKLQFxEJIQp9EZEQUq/QN7MYM3vFzLLMLNPMzjCzWDN728w2eLedvHXNzB41sxwz+8LMRjbvEEREpL7qe6Q/F1jinEsGhgOZwO3AcufcQGC59xjgYmCg9zUNeKJJKxYRkUazE52cZWYdgXVAP3fEymaWDYx1zu0wszhgpXMuycye8u7/8+j1mm0UIiJSL+H1WKcfsBt4zsyGA6uBXwLdDwe5F/zdvPXjge1HbJ/rLfta6JvZNGp/EqBdu3ajkpOTT2YcIiIhZ/Xq1Xucc10bsk19Qj8cGAn83Dn3iZnN5f+ncupidSz7xo8Tzrl5wDyA1NRUt2rVqnqUIiIih5nZ1oZuU585/Vwg1zn3iff4FWo/BHZ50zp4twVHrN/7iO17AfkNLUxERJreCUPfObcT2G5mSd6i8UAGkAZM9ZZNBd7w7qcBN3hdPGOAIs3ni4j4h/pM7wD8HPiHmbUGNgE3UfuB8bKZ3QxsA6Z46y4GJgI5QIm3roiI+IF6hb5z7nMgtY6nxtexrgN+epJ1iYhIM9AZuSIiIUShLyISQhT6IiIhRKEvIhJCFPoiIiFEoS8iEkIU+iIiIUShLyISQhT6IiIhRKEvIhJCFPoiIiFEoS8iEkIU+iIiIUShLyISQhT6IiIhRKEvIhJCFPoiIiFEoS8iEkIU+iIiIUShLyISgNZtL2zUdvX6w+giIuIfdh8s549Lsvif1bmN2l5H+iIiAaCiqoan39vEuAdX8vrnedxyTr9GvY6O9EVE/NzK7ALuWpjBpt3FjE3qyh8mDaFf1/bc0YjXUuiLiPipLXuKmb0og2WZBSR2bsuzN6YyLrn7Sb2mQl9ExM8Ul1fx+Ds5/PX9zUSEGbddlMz3z06kTXjYSb+2Ql9ExE8453jj83zuezOTXQfK+XZKPLddnEz3jpFN9h4KfRERP/BVXhF3pqWzeut+hsVH85frRjGqT6cmfx+FvoiID+09VM6Db2Uz/7PtxLZtzQPfGcaUUb1p1cqa5f0U+iIiPlBZXcOLH23lkWXrKa2o5vtn9eUX4wcSHRXRrO+r0BcRaWEfbNjDrAXpbCg4xLcGduHOyUMY0K1Di7y3Ql9EpIVs31fC7EUZLE3fRUJsW56+IZXzB3fDrHmmcuqi0BcRaWYlFVU8sXIjT723iTAzpk9I4uaz+xIZcfItmA2l0BcRaSbOORZ8sYP7Fmeyo6iMy0b05PaLk4mLjvJZTQp9EZFmkJ5fxKwFGXy6eR9De3bk0WtSOC0x1tdlKfRFRJrS/uIKHno7m5c+2UZM29bc9+1hXJXam7BmasFsqHqFvpltAQ4C1UCVcy7VzGKBfwGJwBbgKufcfqv9jcRcYCJQAtzonFvT9KWLiPiPquoaXvp0Gw+9tZ5D5VXccEYivz5/ENFtm7cFs6EacqR/nnNuzxGPbweWO+fuN7Pbvce3ARcDA72v04EnvFsRkaD00ca9zFqQTtbOg5zZvzN3Th5KUo+WacFsqJOZ3rkMGOvdfx5YSW3oXwa84JxzwMdmFmNmcc65HSdTqIiIv8ndX8J9i7NY9OUO4mOiePJ7I5kwtEeLtmA2VH1D3wFvmZkDnnLOzQO6Hw5y59wOM+vmrRsPbD9i21xv2ddC38ymAdMAEhISGj8CEZEWVlZZzZPvbuSJlRsxg1+fP4hbzu3nkxbMhqpv6J/lnMv3gv1tM8s6zrp1fcS5byyo/eCYB5CamvqN50VE/I1zjiVf7WT2okzyCku55NQ47pg4mPgY37VgNlS9Qt85l+/dFpjZa8BoYNfhaRsziwMKvNVzgd5HbN4LyG/CmkVEWlz2zoPMWpDOvzfuJblHB+ZPG8OYfp19XVaDnTD0zawd0Mo5d9C7fyFwF5AGTAXu927f8DZJA35mZvOp/QVukebzRSRQFZVU8siy9bz48VY6RIZz92VDuWZ0AuFhgfknxutzpN8deM37xUQ48JJzbomZfQa8bGY3A9uAKd76i6lt18yhtmXzpiavWkSkmVXXOOZ/to0Hl2ZTVFrJtacncOsFSXRq19rXpZ2UE4a+c24TMLyO5XuB8XUsd8BPm6Q6EREf+HTzPmampZOx4wCn941l5qVDGRzX0ddlNQmdkSsi4skvLOW+N7NYsC6fntGRPHZtCpcMi/PrFsyGUuiLSMgrq6zmr+9v4vF3NlLjHL8YP5Afn9ufqNb+34LZUAp9EQlZzjneytjF7EUZbN9XysWn9OCOiYPpHdvW16U1G4W+iISknIKDzFqQwfsb9jCoe3te+sHpnDmgi6/LanYKfREJKUWllcxdtoEXPtpC29ZhzLp0KNedHrgtmA2l0BeRkFBd4/ifVduZszSbfSUVXH1aAtMnJBEb4C2YDaXQF5Ggt3rrPmamZfBlXhGnJXbi+cmjOSU+2tdl+YRCX0SC1q4DZdz/Zhavrc2jR8dI5l49gkuH9wyqFsyGUuiLSNApr6rmmQ8289iKHKqqHT89rz8/GTuAdm0UefoXEJGg4ZxjeWYBdy/KYOveEi4Y0p3fXTKYPp3b+bo0v6HQF5GgsHH3Ie5akMG763fTv2s7Xvj+aM4Z1NXXZfkdhb6IBLSDZZU8unwDz324haiIMH4/aQg3nNGHiBBpwWwohb6IBKSaGscra3L545Js9haXc9Wo3vzXhCS6dmjj69L8mkJfRALO59sLuTMtnXXbCxmZEMOzN6Zyaq8YX5cVEBT6IhIwCg6W8ccl2byyOpeuHdrw8FXDuXxEPK1ahW4LZkMp9EXE71VU1fC3f2/m0eU5lFdV86Nz+/OzcQNorxbMBtO/mIj4tXeyC7h7QQab9hQzPrkbv5s0hL5d1ILZWAp9EfFLW/YUc/fCDJZnFdCvSzueu+k0zkvq5uuyAp5CX0T8yqHyKh5bkcMzH2yiTXgYd0xM5sYz+9I6XC2YTUGhLyJ+wTnH65/ncd/iLAoOlnPlqF789qIkunWI9HVpQUWhLyI+90VuITPT0lmzrZDhvaJ56vpRpCR08nVZQUmhLyI+s+dQOXOWZPPy6u10bteGP155KleO7KUWzGak0BeRFldZXcMLH23lT8vWU1pRzQ/O7ssvxg+kQ2SEr0sLegp9EWlR72/YzawFGeQUHOLcQV35/aQhDOjW3tdlhQyFvoi0iG17S5i9KIO3MnbRp3NbnpmayrjkbiH9B018QaEvIs2quLyKv6zM4en3NxPeyvjtRUncfHZf2oSH+bq0kKTQF5Fm4ZwjbV0+9y3OYueBMq5Iiee2i5LpEX3iFszX1+YxZ2k2+YWl9IyJYvqEJC5PiW+BqoOfQl9EmtxXeUXMWpDOZ1v2c0p8Rx67NoXUxNh6bfv62jxmvPolpZXVAOQVljLj1S8BFPxNQKEvIk1mX3EFD76VzT8/3Uantq25/9vDmJLam7AGtGDOWZr9n8A/rLSymjlLsxX6TUChLyInraq6hn98so2H3sqmuKKaG89M5FfnDyI6quEtmPmFpQ1aLg2j0BeRk/LvnD3MWpBB9q6DnD2gC3dOHsLA7h0a/Xo9Y6LIqyPge8ZEnUyZ4lHoi0ijbN9Xwr2LM3nzq5306hTFU9eP4sIh3U+6BXP6hKSvzekDREWEMX1C0smWLCj0RaSBSiuqeeLdjTz17kZamXHrBYP44Tn9iIxomhbMw/P26t5pHgp9EakX5xyLv9zJPYsyyC8qY9KpcdwxcXCzTLtcnhKvkG8m9Q59MwsDVgF5zrlJZtYXmA/EAmuA651zFWbWBngBGAXsBb7rnNvS5JWLSIvJ3HGAmWnpfLJ5H4PjOvLId0dwer/Ovi5LGqEhR/q/BDKBjt7jB4BHnHPzzexJ4GbgCe92v3NugJld7a333SasWUROQkNOfNpfXMHDb6/nH59spWNUBLMvP4VrRic0qAVT/Eu9Qt/MegGXAPcAv7Ha39SMA671VnkemElt6F/m3Qd4BXjMzMw555qubBFpjPqe+FRd43jp09oWzAOllXxvTB9+c8EgYtq29knd0nTqe6T/J+C3wOE+rM5AoXOuynucCxz+jokHtgM456rMrMhbf8+RL2hm04BpAAkJCY2tX0QaoD4nPn28aS8z09LJ2nmQMf1iuXPyUAbHdazr5SQAnTD0zWwSUOCcW21mYw8vrmNVV4/n/n+Bc/OAeQCpqan6KUCkBRzvxKf8wlLuXZzJwi92EB8TxePXjmTisB66CmaQqc+R/lnApWY2EYikdk7/T0CMmYV7R/u9gHxv/VygN5BrZuFANLCvySsXkQY71olP7SPDGffQSpyDX50/kFvO6U9Ua10FMxid8M/LO+dmOOd6OecSgauBFc6564B3gCu91aYCb3j307zHeM+v0Hy+iH+YPiGJqKP66Q04WFbFuORuLL/1XH51/iAFfhA7mT7924D5ZjYbWAs84y1/BnjRzHKoPcK/+uRKFJGmcnje/t7FmRQcLAegR3QkD101nDP7d/FladJCGhT6zrmVwErv/iZgdB3rlAFTmqA2EWliRaWVrMstZG9xBdFREfzmgkFcd3oC4WEn/KFfgoTOyBUJAdU1jpdXbWfO0mz2l1Rw7egEbr0widh2asEMNQp9kSC3ass+Zi5I56u8A5yW2Ik7J4/mlPhoX5clPqLQFwlSO4vKuP/NTF7/PJ+46EgevSaFyafGqQUzxCn0RYJMeVU1f31/M4+/k0NVjePn4wbw47H9adta/91FoS8SNJxzLMssYPaiDLbuLWHC0O787pIh9I5t6+vSxI8o9EWCQE7BIe5amMF763czsFt7/n7z6Zw9UC2Y8k0KfZEAdqCskkeXbeBv/95CVOsw/jBpCNef0YcItWDKMSj0RQJQTY3jlTW5/HFJFnuLK7j6tN7cemESXdq38XVp4ucU+iIBZs22/cxKS2ddbhEjE2J47sbRDOulFkypH4W+SIAoOFDG/UuyeHVNHt07tuGR7w7n8hHxasGUBlHoi/i58qpqnvtwC39evoHKasdPxvbnp+cNoF0b/feVhtN3jYgfeyergLsWZrB5TzHnD+7G7y4ZQmKXdr4uSwKYQl/ED23afYi7F2bwTvZu+nVtx99uOo2xSd18XZYEAYW+iB85VF7Fn1ds4NkPNtMmPIz/njiYqWcm0jpcLZjSNBT6In6gpsbx2to87l+Sxe6D5VyV2ovpE5Lp2kEtmNK0FPoiPrZueyEzF6SzdlshI3rH8PQNqYzoHePrsiRIKfRFfGT3wXLmLM3i5VW5dO3QhoemDOeKlHhatVILpjQfhb5IC6uoquGFj7Ywd9kGyqqqueWcfvxs3AA6REb4ujQJAQp9kRb07vrd3LUgnY27izkvqSu/nzSEfl3b+7osCSEKfZEWsHVvMXcvzGRZ5i4SO7flmampjB/c3ddlSQhS6Is0o+LyKh5/J4e/vr+ZiDDj9ouTuemsRNqEh/m6NAlRCn2RZuCcI21dPvcuzmTXgXK+nRLPbRcn071jpK9LkxCn0BdpYl/lFTEzLZ1VW/czLD6av1w3ilF9Ovm6LBFAoS/SZPYeKufBt9Yz/7NtxLZtzQPfGcaUUb3Vgil+RaEvcpIqq2t48aOtPLJsPaUV1Xz/rL78YvxAoqPUgin+R6EvchI+zNnDzLR0NhQc4lsDu/CHSUMY2L2Dr8sSOSaFvkgjbN9XwuxFGSxN30Xv2CjmXT+KC4Z01x80Eb+n0BdpgJKKKp5cuZEn39tEmBn/deEgfvCtfkRGqAVTAoNCX6QenHMs/GIH9y7OZEdRGZcO78mMicnERUf5ujSRBlHoi5xARv4BZi5I59PN+xgS15G5V6cwum+sr8sSaRSFvojn9bV5zFmaTX5hKT1jovjJ2P5k7TzIPz7ZSnRUBPdccQpXn5ZAmFowJYAp9EWoDfwZr35JaWU1AHmFpfz361/RyuCGMxL59fmDiG6rFkwJfAp9EWDO0uz/BP6RurRvw8xLh/qgIpHmoT+8KULtkX1ddh8sb+FKRJqXjvQlpJVVVvPkuxuP+XzPGHXnSHA54ZG+mUWa2admts7M0s1slre8r5l9YmYbzOxfZtbaW97Ge5zjPZ/YvEMQaTjnHIu/3MH4h97lT8s2MKJ3DG3Cv/7fISoijOkTknxUoUjzqM+Rfjkwzjl3yMwigA/M7E3gN8Ajzrn5ZvYkcDPwhHe73zk3wMyuBh4AvttM9YsA3+y8mT4hictT4utcN3vnQWampfPRpr0k9+jA/GljGNOvc4NeQyRQmXOu/iubtQU+AH4MLAJ6OOeqzOwMYKZzboKZLfXuf2Rm4cBOoKs7zhulpqa6VatWndRAJHQd3XkDtUfp93172NdCu6ikkkeWrefFj7fSITKcWy8YxDWjEwgP06+2JDCZ2WrnXGpDtqnXnL6ZhQGrgQHA48BGoNA5V+Wtkgsc/t8VD2wH8D4QioDOwJ6jXnMaMA0gISGhITWLfE1dnTelldXMWZrN5SnxVNc45n+2jQeXZlNUWsl1p/fhNxcMolO71j6qWMR36hX6zrlqYISZxQCvAYPrWs27revMlW8c5Tvn5gHzoPZIv17VSshoyFRL/jE6b/ILS/l08z5mpqWTseMAo/vGMuvSoQyO69icpYv4tQZ17zjnCs1sJTAGiDGzcO9ovxeQ762WC/QGcr3pnWhgX9OVLMGurhOlZrz6JUCdwd8zJqrOlsvIiDCueuoj4qIjeezaFC4ZFqerYErIq0/3TlfvCB8ziwLOBzKBd4ArvdWmAm9499O8x3jPrzjefL7I0Y43XVOX6ROSiKrjKpeV1TX8YvxAVtw6lkmn9lTgi1C/I/044HlvXr8V8LJzbqGZZQDzzWw2sBZ4xlv/GeBFM8uh9gj/6maoW4LY8aZr6nJ5SjzOOe5emMm+kgoATu0VzePXjqR3bNtmq1MkEJ0w9J1zXwApdSzfBIyuY3kZMKVJqpOQdKzpmmOdKJVTcJBX1+axr6SCQd3bc+fkoZw1oEtzlykSkHRGrvid6ROS6mzBPPpEqQNllcxdtoHn/72Ftq3DuHPyEL43pg8RasEUOSaFvvidw7+sPVb3Tk2N439Wb2fO0mz2Fldw9WkJ/NeFg+jcvo0vyxYJCAp98UuXp8TX2amzeut+Zi1I54vcIlL7dOJvN43mlPhoH1QoEpgU+hIQdh0o44E3s3h1bR49OkYy9+oRXDpcHTkiDaXQF79WXlXNsx9s4bEVG6isdvxkbH9+et4A2rXRt65IY+h/jvitFVm7uGtBBlv2lnDBkO787pLB9OncztdliQQ0hb74nU27D3HXwgxWZu+mf9d2PP/90Zw7qKuvyxIJCgp98RsHyyr584ocnvtwM5HhYfzuksFMPTNRLZgiTUihLz5XU+P43zW5PLAkm73F5UwZ1YvpE5Lp2kEtmCJNTaEvPrV2235mLshg3fZCUhJieGZqKsN7x/i6LJGgpdAXnyg4WMYfl2TzyupcunZow0NThnNFSjytWqkFU6Q5KfSlRVVU1fD8v7cwd/kGyququeWcfvx8/EDaqwVTpEXof5q0mJXZBdy1MINNu4sZl9yN310ymH5d2/u6LJGQotCXZrdlTzF3L8xgeVYBfbu049kbUxmX3N3XZYmEJIW+NJtD5VU8tiKHZz/YTESYMePiZG46qy+tw9WCKeIrCn1pcs45Xv88j/sWZ1FwsJzvjOzFbRcl0a1jpK9LEwl5Cn1pUl/mFnFn2les2VbI8F7RPHX9KFISOvm6LBHxKPSlSew5VM6DS7P516rtdG7XhjlXnsp3RvZSC6aIn1Hoy0mprK7hhY+28qdl6ymtqOYHZ/fl5+MH0jEywteliUgdFPrSaO9v2M2sBRnkFBzinEFd+cOkIQzophZMEX+m0JcG27a3hNmLMngrYxd9Orfl6RtSOX9wN/1BE5EAoNCXeiupqOIv72xk3vubCG9l/PaiJG4+uy9twsN8XZqI1JNCX07IOUfaunzuW5zFzgNlXDaiJzMuHkyPaLVgigQahb4cV3p+EbPSMvh0yz5Oie/IY9emkJoY6+uyRKSRFPpSp33FFTz0Vjb//HQbMW1bc+8Vw/juab0JUwumSEBT6MvXVFXX8PePt/Lw2+sprqhm6pmJ/Gr8IKLbqgVTJBgo9OU//p2zh1kLMsjedZCzB3ThD5OHMKh7B1+XJSJNSKEvbN9Xwr2LM3nzq5306hTFk98bxYSh3dWCKRKEFPohrLSimiff3ciT727EDG69YBA/PKcfkRFqwRQJVgr9EOSc482vdnLPokzyCkuZPLwnMy5OpmdMlK9LE5FmptAPMVk7DzAzLZ2PN+1jcFxHHr5qOKf36+zrskSkhSj0Q0RhSQUPv72ev3+8leioCGZffgrXjE5QC6ZIiFHoB7nqGsdLn27jobeyOVBayffG9OE3Fwwipm1rX5cmIj6g0A9in2zay8wFGWTuOMCYfrHMvHQoyT06+rosEfGhE4a+mfUGXgB6ADXAPOfcXDOLBf4FJAJbgKucc/utts9vLjARKAFudM6taZ7ypS75haXcuziThV/sID4misevHcnEYT3Ugiki9TrSrwJudc6tMbMOwGozexu4EVjunLvfzG4HbgduAy4GBnpfpwNPeLfSzMoqq3n6vU08vjIH5+CX4wfyo3P7E9W66VswX1+bx5yl2eQXltIzJorpE5K4PCW+yd9HRJrWCUPfObcD2OHdP2hmmUA8cBkw1lvteWAltaF/GfCCc84BH5tZjJnFea8jzcA5x9L0ncxelEnu/lImDuvBHRMH06tT22Z5v9fX5jHj1S8prawGIK+wlBmvfgmg4Bfxcw2a0zezRCAF+ATofjjInXM7zKybt1o8sP2IzXK9ZV8LfTObBkwDSEhIaETpArB+10FmLUjnw5y9JHXvwEs/PJ0z+3dp1vecszT7P4F/WGllNXOWZiv0RfxcvUPfzNoD/wv8yjl34Djzw3U94b6xwLl5wDyA1NTUbzwvx1dUUskjy9bz4sdbadc6jJmTh/C9MX0ID2vV7O+dX1jaoOUi4j/qFfpmFkFt4P/DOfeqt3jX4WkbM4sDCrzluUDvIzbvBeQ3VcGhrrrG8fKq7cxZms3+kgquHZ3ArRcmEduu5Vowe8ZEkVdHwOuMXhH/d8LDQq8b5xkg0zn38BFPpQFTvftTgTeOWH6D1RoDFGk+v2ms2rKPyx7/gBmvfsmAru1Z+POzueeKYS0a+ADTJyQRddT1eaIiwpg+IalF6xCRhqvPkf5ZwPXAl2b2ubfsDuB+4GUzuxnYBkzxnltMbbtmDrUtmzc1acUhaGdRGfe9mckbn+fTo2Mkj16TwuRT43zWgnl43l7dOyKBx2qbbHwrNTXVrVq1ytdl+J2yymqe+WAzj7+TQ1WNY9q3+vGT8/rTtrXOqRMRMLPVzrnUhmyj9PBDzjmWZRYwe1EGW/eWMGFod/574hASOjdPC6aIhA6Fvp/JKTjEXQszeG/9bgZ0a8+LN4/mWwO7+rosEQkSCn0/caCskj8v38BzH24hKiKM308awg1n9CGiBVowRSR0KPR9rKbG8cqaXP64JIu9xRVcNao30y9Kokv7Nr4uTUSCkELfh9Zs28+stHTW5RYxMiGG524czbBe0b4uS0SCmELfBwoOlPHAkmz+d00u3Tq04ZHvDufyEfG6CqaINDuFfguqqKrhuQ838+jyDVRWO348tj8/PW8A7dtoN4hIy1DatJB3sgu4e0EGm/YUMz65G7+fNITELu18XZaIhBiFfjPbvKeYuxdmsCKrgH5d2vHcTadxXlK3E28oItIMFPrN5FB5FX9esYFnP9hMm/Aw7piYzI1n9qV1uFowRcR3FPpNrKbG8draPO5fksXug+VcOaoXv70oiW4dIn1dmoiIQr8prdteyMwF6azdVsjw3jHMu34UKQmdfF2WiMh/KPSbwO6D5cxZmsXLq3Lp0r4ND04ZzrdT4mnVSi2YIuJfFPonoaKqhhc+2sLcZRsoq6rmlnP68bNxA+gQGeHr0kRE6qTQb6R31+/mrgXpbNxdzLmDuvKHyUPo37W9r8sSETkuhX4Dbd1bzN0LM1mWuYvEzm15Zmoq45K76WxaEQkICv16Ki6v4i8rc3j6vc1EhBm3XZTM989OpE142Ik3FhHxEwr9E3DOkbYun3sXZ7LrQDnfTonntouT6d5RLZgiEngU+sfxVV4RM9PSWbV1P8Pio/nLdaMY1UctmCISuBT6ddh7qJwH31rP/M+2Edu2NQ98ZxhTRvVWC6aIBDyF/hEqq2v4+8dbeeTt9ZRUVHPTmX355fkDiY5SC6aIBAeFvufDnD3MWpDO+l2H+NbALvxh0hAGdu/g67JERJpUyIf+9n0l3LMokyXpO+kdG8VT14/iwiHd1YIpIkEpZEO/tKKaJ1bm8NR7m2hlxvQJSdx8dl8iI9SCKSLBK+RC3znHwi92cO/iTHYUlXHp8J7MmJhMXHSUr0sTEWl2IRX6GfkHmLkgnU8372NIXEcevSaF0xJjfV2WiEiLCYnQ319cwcNvr+cfn2wlOiqCe644hatPSyBMLZgiEmKCOvSrqmv456fbePCt9Rwqr+KGMxL59fmDiG6rFkwRCU1BG/ofbdzLrAXpZO08yBn9OjPz0qEk9VALpoiEtqAL/bzCUu5dnMmiL3YQHxPFE9eN5KJTeqgFU0SEIAr9sspqnnp3E0+8m4Nz8KvzB3LLOf2Jaq0WTBGRwwI+9J1zLPlqJ7MXZZJXWMolw+KYMTGZXp3a+ro0ERG/E9Chn7XzALPSMvho016Se3Tgnz8cwxn9O/u6LBERvxWQoV9UUsnDb2fz90/5g1BgAAAFc0lEQVS20SEynLsvG8o1oxMID2vl69JERPxaQIV+dY1j/mfbeHBpNkWllVx7egK3XpBEp3atfV2aiEhAOGHom9mzwCSgwDl3ircsFvgXkAhsAa5yzu232haZucBEoAS40Tm3pikK/WzLPu58I52MHQcY3TeWmZOHMqRnx6Z4aRGRkFGf+ZC/ARcdtex2YLlzbiCw3HsMcDEw0PuaBjxxsgXuKCrlF/9cy5QnP2J/SQWPXZvCv6aNUeCLiDTCCY/0nXPvmVniUYsvA8Z6958HVgK3ectfcM454GMzizGzOOfcjoYWVlZZzTMfbOaxFTlUO8cvxg3gR2P707Z1QM1IiYj4lcYmaPfDQe6c22Fm3bzl8cD2I9bL9ZZ9I/TNbBq1Pw0AlJvZV8d7w1vvgVsbWawf6ALs8XURzSiYxxfMYwONL9AlNXSDpj5sruu0V1fXis65ecA8ADNb5ZxLbeJa/IbGF7iCeWyg8QU6M1vV0G0a2+O4y8zivDeNAwq85blA7yPW6wXkN/I9RESkiTU29NOAqd79qcAbRyy/wWqNAYoaM58vIiLNoz4tm/+k9pe2XcwsF7gTuB942cxuBrYBU7zVF1PbrplDbcvmTfWsY17Dyg44Gl/gCuaxgcYX6Bo8PqtttBERkVCg6xaIiIQQhb6ISAhp8dA3s0gz+9TM1plZupnN8pb3NbNPzGyDmf3LzALugjrHGdvfzGyzmX3ufY3wda0nw8zCzGytmS30Hgf8vjtSHeMLmv1nZlvM7EtvHKu8ZbFm9ra3/942s06+rrOxjjG+mWaWd8T+m+jrOhvDO9n1FTPLMrNMMzujMfvOF0f65cA459xwYARwkdfp8wDwiHdph/3AzT6o7WQda2wA051zI7yvz31XYpP4JZB5xONg2HdHOnp8EFz77zxvHIf71491WZVAdfT4oPb78/D+W+yzyk7OXGCJcy4ZGE7t92iD912Lh76rdch7GOF9OWAc8Iq3/Hng8pau7WQdZ2xBw8x6AZcAf/UeG0Gw7w47enwh4jJq9xsE+P4LVmbWETgHeAbAOVfhnCukEfvOJ3P63o/Pn1N7UtfbwEag0DlX5a1y+PINAefosTnnPvGeusfMvjCzR8ysjQ9LPFl/An4L1HiPOxMk+85z9PgOC5b954C3zGy1dykUOOqyKkC3Y27t/+oaH8DPvP33bIBOX/UDdgPPeVOPfzWzdjRi3/kk9J1z1c65EdSesTsaGFzXai1bVdM4emxmdgowA0gGTgNiqb04XcAxs8OX2F595OI6Vg3IfXeM8UGQ7D/PWc65kdReEfenZnaOrwtqYnWN7wmgP7VTrjuAh3xYX2OFAyOBJ5xzKUAxjZyG82n3jvfjyUpgDBBjZodPFgv4yzccMbaLnHM7vKmfcuA5aj/oAtFZwKVmtgWYT+20zp8Inn33jfGZ2d+DaP/hnMv3bguA16gdy7EuqxJw6hqfc26XdzBWAzxNYO6/XCD3iJmDV6j9EGjwvvNF905XM4vx7kcB51P7C4l3gCu91Y68tEPAOMbYso7YKUbtnNtxryjqr5xzM5xzvZxzicDVwArn3HUEwb6DY47ve8Gy/8ysnZl1OHwfuJDasRzrsioB5VjjO7z/PFcQgPvPObcT2G5mh6+qOR7IoBH7zhcXp48DnjezMGo/dF52zi00swxgvpnNBtbi/cIiwBxrbCvMrCu1UyGfAz/yZZHN4DYCf98dzz+CZP91B16r/ewiHHjJObfEzD6j7suqBJpjje9Fr83WUfuX/m7xXYkn5efUfi+2BjZRe5mbVjRw3+kyDCIiIURn5IqIhBCFvohICFHoi4iEEIW+iEgIUeiLiIQQhb6ISAhR6IuIhJD/Aw4qIxjjwjWdAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "plt.xlim(30, 60)\n", "plt.ylim(100, 600)\n", "\n", "X = np.array([39.93, 42.05, 43.18, 44.68, 49.87, 53.57])\n", "Y = np.array([199, 290, 298, 310, 399, 420])\n", "plt.scatter(X, Y)\n", "\n", "x = np.linspace(30, 60, 100)\n", "y = 15*x-360\n", "plt.plot(x, y)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "求解一元线性回归，实际上就是找一条直线，使得它和散点图中的所有点距离最近。直线和散点的距离一般通过损失函数来度量，**损失函数**（Loss Function）有时也称为**代价函数**（Cost Function），常用的损失函数有两种，一种是 **绝对损失（Absolute Loss）**，另一种是 **平方损失（Squared Loss）**：\n", "\n", "绝对损失函数：$$ loss = |y - \\hat{y}| $$\n", "\n", "平方损失函数：$$ loss = (y - \\hat{y})^2 $$\n", "\n", "其中，$ \\hat{y} $ 表示预测值，由于它是一元线性的，我们可以记为：$ \\hat{y} = ax + b $，这条直线也叫 **回归线**。这条回归线就是我们要学习的对象，可以把它叫做**模型**（model），它反映了数据的某种潜在规律，因此也称为**假设**（hypothesis）。\n", "\n", "在一元线性回归的实际求解中，因为平方损失处处可导，可以方便后面的计算，常常使用平方损失函数：\n", "\n", "$$\n", "\\begin{equation}\n", "\\begin{aligned}\n", "loss &= (y - \\hat{y})^2 \\\\\n", "&= (y - (ax + b))^2 \\\\\n", "&= \\sum_{i=1}^{n}(y_i - (ax_i + b))^2\\\\\n", "&= (y_1 - (ax_1 + b))^2 + (y_2 - (ax_2 + b))^2 + ... + (y_n - (ax_n + b))^2\n", "\\end{aligned}\n", "\\end{equation}\n", "$$\n", "\n", "将其完全展开可得：\n", "\n", "$$\n", "loss = (y_1^2 + y_2^2 + ... + y_n^2) + a^2(x_1^2 + x_2^2 + ... + x_n^2) + nb^2 + 2ab(x_1 + x_2 + ... + x_n) - 2a(x_1y_1 + x_2y_2 + ... + x_ny_n) -2b(y_1 + y_2 + ... + y_n)\n", "$$\n", "\n", "令 $ \\bar{x} = \\frac{x_1 + x_2 + ... + x_n}{n} $，表示 x 的平均数，则：\n", "\n", "$$\n", "loss = n\\bar{y^2} + a^2n\\bar{x^2} + nb^2 + 2abn\\bar{x} - 2an\\bar{xy} - 2bn\\bar{y}\n", "$$\n", "\n", "很显然，这是一个关于 a、b 变量的二元二次函数，我们记为：\n", "\n", "$$\n", "F(a, b) = n\\bar{y^2} + a^2n\\bar{x^2} + nb^2 + 2abn\\bar{x} - 2an\\bar{xy} - 2bn\\bar{y}\n", "$$\n", "\n", "那么问题就转变成了求函数 F(a, b) 的最小值。为了求解这个问题，我们先看看一元函数的极值问题。" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4XPV97/H3d7Tv+77LluXdeAE7gEPYHdNgAlkgCSE3UEp723S7TdPm3j637fPc0qS3vW3ThCRAgJASQgiEJJBgNhu8Ihvvlq19sfbF2reZ+d0/ZuQqQrZG0sycWb6v5/GjseZI56sj6aPffM/v/I4YY1BKKRX8bFYXoJRSyjs00JVSKkRooCulVIjQQFdKqRChga6UUiFCA10ppUKEBrpSSoUIDXSllAoRGuhKKRUiIv25s8zMTFNaWurPXSqlVNA7cuRIjzEma77t/BropaWlVFVV+XOXSikV9ESkyZPttOWilFIhQgNdKaVChAa6UkqFCA10pZQKERroSikVIuYNdBF5UkS6ROTUrPf/kYicE5HTIvIN35WolFLKE56M0J8Cdsx8h4jcCOwC1htj1gD/5P3SlFJKLcS8gW6M2Qv0zXr37wOPGmMm3Nt0+aC2S96q7uTb79T6chdKKeUT3UMTfOPX1dR3D/t8X4vtoa8AtovIIRHZIyJXX25DEXlYRKpEpKq7u3tRO9tX28u/v1mL06n3P1VKBZdzHUN8+506OgcnfL6vxQZ6JJAGbAP+AviJiMhcGxpjvmeM2WKM2ZKVNe+Vq3Mqz0pgbMpBx+D4IstVSilr1Pe4RubLshJ8vq/FBnor8DPjchhwApneK+u3lWW6DkRDz4ivdqGUUj5R3z1CQnQEWUkxPt/XYgP9ZeAmABFZAUQDPd4qarbyzEQA6jXQlVJBpqFnhLKsBC7TxPAqT6YtPgccACpFpFVEHgSeBMrdUxl/DDxgjPFZgzsnOYb46Ai/nFRQSilvqu8ZvjQo9bV5V1s0xtx3mae+4OVaLktEKMtM0JaLUiqoTNgdtPaPcffGQr/sL2iuFC3LTKC+WwNdKRU8mnpHMcY1scMfgibQy7MSae0fZcLusLoUpZTyyPQg1F8tl+AJ9MwEnAZa+katLkUppTwy3SYuzYz3y/6CJtCnpy7WadtFKRUk6ruHyUqKISk2yi/7C55Az9K56Eqp4NLQM0J5pn/65xBEgZ4cG0VmYoxOXVRKBY36nhG/nRCFIAp0cPXRdYSulAoGF0cn6RuZvNQu9ofgCvQsDXSlVHCYzip/zXCBIAv0sswEeoYnGRibsroUpZS6oukpi2XacpmbLtKllAoWDT0jRNiEojT/TFmEIAv08izXS5eGHj0xqpQKbA09IxSnxxMd6b+YDapAL06PxyboEgBKqYBX1z3s1xOiEGSBHh1poyg9XpfRVUoFNKfT0Ng7ooE+H12kSykV6NoHxxmfcmqgz6c8M5GGnmG9v6hSKmA1XFqUSwP9ipZlJzA+5aRtYMzqUpRSak61XUMALM/23xx0CMJAX+6e6VLbpTNdlFKBqbZ7mKTYSL/cR3SmoAv0ipwkQANdKRW4ajqHqchO9Mt9RGcKukBPT4gmPSFaA10pFbDquof93m6BIAx0cPWlNNCVUoGof2SSnuFJKrKT/L7voA30mq5hjNGZLkqpwFLrXuJbR+geWp6VyMDYFD3Dk1aXopRSv2W6e6CB7qGKHNeBqnFPDVJKqUBR0zlMbJSNgtQ4v+87KAN9+i9fnfbRlVIBprZ7mGVZidhs/p3hAkEa6LnJsSTGROqJUaVUwKnrck1ZtEJQBrqIsMx9YlQppQLFyISdCxfHLOmfQ5AGOkCFTl1USgWYukszXPw/ZRGCONCXZyfSNTSht6NTSgUMK2e4QDAHuq7popQKMDVdw0TahJIM/912bqagDfTpqYs600UpFShqu1x3KYqKsCZagzbQC9Nc9+qbvipLKaWsVtdlzRou04I20CNswrKsRGo69eIipZT1JuwOGntHLJuyCEEc6OBepEtH6EqpANDYM4rTwDIN9MWpyE6ktX+MsUmH1aUopcJcjUV3KZopqAN9RU4ixuiaLkop653vGLrUCrbKvIEuIk+KSJeInJrjuf8hIkZEMn1T3pVV5iYDUN2hga6UslZ1xxClGfHERkVYVoMnI/SngB2z3ykiRcCtQLOXa/JYcXo8sVE2zmmgK6Usdq5ziMpca64QnTZvoBtj9gJ9czz1L8BXAcvuMhFhEyqykzivM12UUhYanbTT3DdKZU6ypXUsqocuIncCF4wxxz3Y9mERqRKRqu7u7sXs7ooqc5O05aKUslRN5zDGEPgj9NlEJB74OvA3nmxvjPmeMWaLMWZLVlbWQnc3r5W5SXQPTdA3oncvUkpZY7rtuzLYAh1YBpQBx0WkESgEjopIrjcL89SKHNcBrO4YtGL3SilFdccQsVE2itKtWcNl2oID3Rhz0hiTbYwpNcaUAq3AJmNMh9er88D0X0Q9MaqUssq5zkFW5CQRYcFdimbyZNric8ABoFJEWkXkQd+X5bmspBjS4qP0xKhSyjLnOoapzLG23QIQOd8Gxpj75nm+1GvVLIKI6IlRpZRleocn6BmesPyEKAT5laLTKnOSON8xhNNp2QxKpVSYmm73aqB7SWVuMiOTDi5cHLO6FKVUmKnWQPeu6QOpbRellL+d6xgiPSGarMQYq0sJrUDXE6NKKX871zlEZU4SItbOcIEQCfTEmEgK0+J0hK6U8iun03A+ANZwmRYSgQ6uE6Pn9OIipZQftfaPMTrp0ED3tsrcJOq7R5i0O60uRSkVJqavUNdA97KVecnYnYbaLr0lnVLKP6bbvCsC4KIiCKFAX53nWrbyTLu2XZRS/nGmbZCyzAQSY+a9RtMvQibQyzITiIuK4HTbgNWlKKXCxOn2gUuDyUAQMoEeYRNW5iVxuk1H6Eop3xsYm6Klb4zV+RroPrE6L5mzbYMYo0sAKKV866y7vauB7iNr8lMYmrDT0qdLACilfOuMuxuwRgPdN6YP7Jl27aMrpXzrdNsgWUkxZCfFWl3KJSEV6JW5rgXmtY+ulPK1022BdUIUQizQY6MiWJaVoIGulPKpCbuD2q7hgGq3QIgFOrj66Gc00JVSPlTTOYzdaViTn2J1Kb8l5AJ9dV4yHYPj9A5PWF2KUipETQ8aA2mGC4RgoE+/BNK2i1LKV063DZAQHUFJerzVpfyWkAv01fm6BIBSyrfOtA+yKi8Zm836NdBnCrlAT42PpiA1TkfoSimfcDoNZ9oGA+6EKIRgoINrlH5G13RRSvlAU98oI5OOgDshCqEa6HnJ1PeMMDppt7oUpVSICdQTohCigb4mPxlj/mutBaWU8pbTbQNE2oSKnESrS/mQkAz09YWpAJxo1baLUsq7TrQOsDIviZjICKtL+ZCQDPTclFiyk2I00JVSXmWM4UTrxUuDxkATkoEOrlH68daLVpehlAohjb2jDI7bWV8QeCdEIaQDPYX67hEGx6esLkUpFSJOuAeJOkL3s/WFrr+gp7TtopTykuMtA8RG2VgRgCdEIaQD3X1i9IIGulLKO060XmRNfgqREYEZnYFZlRekJ0RTlB536SWSUkothd3h5FTbwKVX/4EoZAMd3CdGW3SErpRaupquYcannGwI0P45hHigbyhM4cLFMV1KVym1ZNOv9tfpCN0aeoGRUspbjrcOkBQTSVlGgtWlXNa8gS4iT4pIl4icmvG+b4pItYicEJGXRCQgX4OsLUhBBJ2PrpRaspOtA6wrTAm4JXNn8mSE/hSwY9b7dgNrjTHrgfPAX3m5Lq9IjIlkeVaijtCVUksyYXdQ3TEYsPPPp80b6MaYvUDfrPe9boyZXsrwIFDog9q8Yn1hKidaBzDGWF2KUipInW0fYsph2BDA/XPwTg/9y8BrXvg8PrG+MIWe4QnaB8atLkUpFaQuXSFaFOQj9CsRka8DduBHV9jmYRGpEpGq7u7upexuUTa4vwHHWrSPrpRanGMtF8lIiCY/JdbqUq5o0YEuIg8AvwN83lyhn2GM+Z4xZosxZktWVtZid7doq/OSiY60cbSp3+/7VkqFhqNN/WwqSUMkcE+IwiIDXUR2AH8J3GmMGfVuSd4VHWljfUEKR5s10JVSC9c7PEFj7yibS9KsLmVenkxbfA44AFSKSKuIPAh8C0gCdovIMRF5zMd1LsnmkjROXRhkfMphdSlKqSBztNnVrt1UHPiBHjnfBsaY++Z49xM+qMVnNpWk8d299ZxuG2BzSbrV5SilgsiRpn4ibRLQa7hMC+krRadN/2U92qQnRpVSC3O0uZ81BSnERgXeLedmC4tAz0qKoTg9niN6YlQptQBTDifHWy6yOQjaLRAmgQ6uPvqR5n69wEgp5bEzbYNM2J1sKgns+efTwibQNxWn0j00QWv/mNWlKKWCxPTsuGCY4QLhFOjub4hOX1RKeepIUz/5KbHkpcRZXYpHwibQK3OSSIiO0D66Uspj0xcUBYuwCfTICBsbilI10JVSHmkfGKNtYDwo5p9PC5tAB1cfrLpjiJEJ+/wbK6XC2vQ052Dpn0OYBfqmkjQcTsNxXahLKTWPI039xEbZWJ2fbHUpHguvQC9OQwTeb9S2i1Lqyt5v7GNDYSpREcETk8FTqRekxEWxKjeZQw29VpeilApgg+NTnG4bYGt5htWlLEhYBTrA1vJ0jjT1M2HXhbqUUnOrauzDaWBbWXCt/RR+gV6WwYTdqfcZVUpd1qH6PqIihI1BNMMFwjDQr3H/xT1Ur20XpdTcDja4+udx0YG/INdMYRfo6QnRrMxN4lBD3/wbK6XCzvCEnVMXBtgWZP1zCMNAB9halk5VYz9TDqfVpSilAsyRpn4cTsPW8uDqn0O4Bnp5BmNTDu2jK6U+5GB9L5E2CaoLiqaFZaBf6qPr9EWl1CyH6ntZV5hCfPS8N3QLOGEZ6JmJMSzPTuRQvfbRlVL/ZXTSzonWAbaWBV//HMI00AG2ladT1diHXfvoSim3o00XsTsN24Kwfw5hHOhbyzIYmXRwqm3Q6lKUUgHiYH0vETZhS6kGelCZPoN9oE776Eopl4P1vazNTyYxJvj65xDGgZ6dFMuKnET21fZYXYpSKgAMjU/xQctFrlueaXUpixa2gQ6wvSKLw419jE/pui5KhbsDdb04nIbtFVlWl7JoYR7omUzanRzWq0aVCnvv1fYQHx3BppJUq0tZtLAO9K1lGURH2Hi3ptvqUpRSFnu3poetZenERAbX+i0zhXWgx0VHsKU0jXdrtI+uVDhr6RuloWckqNstEOaBDnB9RSbVHUN0DY1bXYpSyiLvuSdHbK8I3hOioIHOR91/kXW2i1Lh672aHnKTY1menWh1KUsS9oG+Oi+Z9IRo3j2vga5UOHI4De/V9nB9RSYiYnU5SxL2gW6zCdctz+Td2h6MMVaXo5Tys1MXBhgYmwr6dgtooAOwfXkm3UMTnOscsroUpZSfTc9yC+YLiqZpoOM6MQpo20WpMLS3poc1+clkJsZYXcqSaaAD+alxVOYk8VZ1l9WlKKX8aGB0iiNN/dywIrinK07TQHe7aVU27zf2MTA2ZXUpSik/eed8Fw6n4eZVOVaX4hXzBrqIPCkiXSJyasb70kVkt4jUuN8G372aZrllVTZ2p2HPeb1qVKlw8VZ1FxkJ0VxVFLyX+8/kyQj9KWDHrPd9DXjTGFMBvOn+f1C7qiiN9IRo3jrbaXUpSik/sDucvHOumxtXZhNhC+7pitPmDXRjzF5g9upVu4Cn3Y+fBu7ycl1+F2ETPlaZxdvnuvUuRkqFgaqmfgbGprh5ZbbVpXjNYnvoOcaYdgD325A4IresymFgzHWSRCkV2t6q7iI6wsb2EDkhCn44KSoiD4tIlYhUdXcHdn96e0UmURHCmzrbRamQ98bZTraWpwft3YnmsthA7xSRPAD328smoDHme8aYLcaYLVlZgf2XMCk2iq1lGbypfXSlQlpDzwj13SMh1W6BxQf6K8AD7scPAD/3TjnWu3lVNnXdIzT2jFhdilLKR6YHbaEyXXGaJ9MWnwMOAJUi0ioiDwKPAreKSA1wq/v/IeHmla5v8Bs6SlcqZL15tovKnCSK0uOtLsWr5m0eGWPuu8xTN3u5loBQnBFPZU4SvzndwUPby60uRynlZX0jkxxu7OORG0Lv91uvFJ3DznV5VDX10zmoN71QKtT85nQHDqdh57o8q0vxOg30OdyxPhdj4LWT7VaXopTysldPtlOaEc/qvGSrS/E6DfQ5LM9OYkVOIq+e7LC6FKWUF/WNTLK/rped6/KC/mYWc9FAv4yd6/J4v6mPLm27KBUyXg/hdgtooF/WHevyXG2XUzpKVypU/OpkOyUZ8azJD712C2igX1ZFThIV2Yn8SvvoSoWE/hBvt4AG+hXtXJfH+43adlEqFLx+xtVuuSNE2y2ggX5Fd6x3tV1+fVrbLkoFu1+d7KA4PXTbLaCBfkUr3G2Xnx9rs7oUpdQS9AxPsK+2hzvWh267BTTQ5/XJTQUcaerXtV2UCmI/P9aGw2m4e2OB1aX4lAb6PD65sQAR+NnRVqtLUUot0otHWllfmEJFTpLVpfiUBvo88lLiuH55Ji8evYDTaawuRym1QGfbBznTPsg9mwqtLsXnNNA9cM+mQi5cHONQw+w78SmlAt2LR1qJihDu3JBvdSk+p4HugdvX5JIYE6ltF6WCjN3h5OVjbdy0Mpu0hGiry/E5DXQPxEVHsHNdLq+ebGd00m51OUopD+2t6aZneCIs2i2gge6xezYVMjLp4Dc6J12poPHikQukJ0TzscrQutXc5Wige+jq0nSK0uN4oUrbLkoFg/6RSXaf7eTODflER4ZH1IXHV+kFNptw79XF7K/rpbZr2OpylFLzeOFIC5N2J/ddU2x1KX6jgb4An726iKgI4dmDTVaXopS6AqfT8OzBZq4pS6cyN7Tnns+kgb4AmYkx7FyXx4tHWvXkqFIBbE9NN819o9y/rcTqUvxKA32B7t9WwtCEXdd3USqAPXugiczEGG5fk2t1KX6lgb5Am0vSWJWXzDMHmjBGrxxVKtC09I3y1rkuPndNUdicDJ0WXl+tF4gI928r4Wz7IEeb+60uRyk1y48ONWMT4b6t4XMydJoG+iLsuiqfpJhInjmgJ0eVCiTjUw5+UtXCLauyyUuJs7ocv9NAX4SEmEg+c3URvzrRTmv/qNXlKKXcXjzaSt/IJA9cW2p1KZbQQF+kB68vA+DxdxssrkQpBa51W767p54NRal8pDzD6nIsoYG+SPmpcdy1sYAfv99M38ik1eUoFfZeO9VBc98ov3/DspC+K9GVaKAvwSM3lDM+5eSp/Y1Wl6JUWDPG8J136ijPSuC21TlWl2MZDfQlWJ6dxG2rc3h6fyMjE3qhkVJW2VvTw5n2QR65YRk2W3iOzkEDfcke+dgyBsameO5ws9WlKBW2vvNOLbnJsdx1VWjfM3Q+GuhLtKk4jW3l6Xz/3XrGpxxWl6NU2Hm/sY+D9X08tL0s7C4kmi28v3ov+bNbK+kcnNBeulJ+ZozhH1+rJjsphs9vDa91W+aige4F15Slc2NlFt9+u5aB0Smry1EqbLxV3UVVUz9/fEsFcdERVpdjOQ10L/nqjpUMTdh5bG+d1aUoFRYcTsM3fn2OsswEPrOlyOpyAoIGupesyktm14Z8frCvgc7BcavLUSrk/fzYBc51DvHnt60gKkKjDJYY6CLypyJyWkROichzIhLrrcKC0Z/dWonDafjXN2usLkWpkDZhd/DPu8+zriCFnWvzrC4nYCw60EWkAPgKsMUYsxaIAO71VmHBqDgjns9vLeH591s40zZodTlKhazH322gtX+Mr+6oDOt557Mt9XVKJBAnIpFAPBD2d33401tWkBoXxf98+SROp66XrpS3tfSN8u9v1bBjTS7bK7KsLiegRC72A40xF0Tkn4BmYAx43RjzutcqC1Ip8VH89c5V/PkLx3m+qiWsblAbiAbHp2juHaW5b5SuwXG6hyfoGZpkeMLOyKSd0QkHjhk3KomOsBEfHUF8TCSpcVFkJsaQlRRDXkosxRnxFKbFEROpsymsYozhf79yGpsIf/OJ1VaXE3AWHegikgbsAsqAi8ALIvIFY8yzs7Z7GHgYoLg4PMLt7k0FPF/VwqOvVXPb6hwyEmOsLinkOZ2G+p4RjrVc5Gz7IGfbBznXMUTvrIXTImxCRkI0yXFRJERHEBcdQbTNFdAGw5Td0DE4zsiEnYtjU1ycNQ1VBIrS4lmZm8SqvGTW5CezsTiNrCT9HvvD7jOdvFndxV/vXEl+avitdz4fWext1ETk08AOY8yD7v9/EdhmjPmDy33Mli1bTFVV1aL2F2zOdw6x81/f5ZMbC/jmpzdYXU7IcToNZ9oH2V/Xw4G6Xo42X2RgzBW+sVE2KnOSWJmbTFlWAiXp8RRnxJObHEtafPSCeq6Tdic9wxO0D4zR1DtKY+8odV3DnO0YpKFnhOlfn6L0OLaUpHPtsgyuW56pYeMDo5N2bv3nvSTGRPLLr1wfVjNbROSIMWbLfNsteoSOq9WyTUTicbVcbgbCI609sCIniYe2l/PYnjp2rsvjxpXZVpcU9AbHp9hzrps3z3ay53w3/e7R87KsBHasyWVzSRobi1Mpz0okwksnyqIjbeSnxpGfGsfmkvTfem5s0sHptgGONvdztOkie89389IHFwAoz0rgpspsbl6Vw5bStLAKH1/5P6+e5cLFMV545CN6PC9j0SN0ABH5W+CzgB34AHjIGDNxue3DaYQOrtth7frWPnpHJvj1n3yUTG29LFj/yCSvn+ngVyc72F/bg91pSIuP4sbKbK6vyOS65ZnkJAfGbFmn03Cuc4h9tT3sOd/Nofo+Jh1OUuKiuG11DjvX53HdssywX29kMd4408lDz1Txu9vL+Pod4dc793SEvqRAX6hwC3SA6o5B7vzWPrYvz+TxB7aE7cL7CzE+5WD3mU5e+uACe893Y3caitPj+fi6XG5dlcPG4jSvjcB9aXjCzns13bx+upPdZzoZmrCTEhfFHevzuHtjAZtL0vTnwQNdQ+N8/P+9S3ZyLC//92vD8qS0P1ouygMrc5P5q4+v5G9/cYZnDzVz/zZdQGguxhiOtw7w/PvN/PJ4O0MTdvJSYnlwexmfWJ/PmvzkoAu/xJhIdqzNY8faPCbsDt6r6eEXx9t46egF/vNQMyUZ8Xx6cyGf2lxEbkpgvMoINMYY/uKFEwxP2PnxvVeFZZgvhI7Q/cAYw5d+8D4H63t58fevZW1BitUlBYzB8SleOnqB5w43U90xRFxUBDvX5XHPpgK2lWeE5EUjwxN2fnOqgxePtrK/rhebwMcqs/ncNcXcuDI7KF59+Mtje+p49LVq/n7XGu7/SKnV5VhGWy4Bpmd4gl3f2ofDaXjlD68jO0D6vlY53TbAswebePmDNsamHKwvTOHeq4v5xIY8kmKjrC7Pb5p7R/lJVQs/qWqha2iCgtQ4Pre1mM9eXRT251zeONPJ7/6wip1r8/jW5zYG3Ss0b9JAD0Bn2we55zv7qchO5Pnf+wixUeH18nHK4eT10508vb+Rw419xEbZ2LWhgC9sK2FdYXi/aplyOHnjTCc/PNjE/rpeoiNs/M6GPL50bSnrC1OtLs/vzrYP8qnv7GdZdiLPP/yRsF8aVwM9QO0+08nDP6zijnV5/Nu9G0OypTDbxdFJnjvcwjMHGmkfGKcoPY4vbivlM1uKSIkPn9G4p2q7hnnmQCMvHmllZNLB5pI0vnxdGbevySEyDKbrdQ9NcNd/7MPudPLKH14fMLOYrKSBHsCm+4Jf2FbM3+9aG7IvJeu7h3lyXwMvHrnA2JSDa5dl8KVrS7l5VY72iT0wOD7FC1WtPLW/gZa+MQpS4/jStaV89poikkO0LdU3Msnnvn+Qxt4RXvi9a8P+lds0DfQAZozh0V9X89099SEX6sYYDtT38uR7DbxxtovoCBt3XpXPl68rY3V+stXlBSWH0/DG2U6eeK+Bww19JERH8Jmri/jydWUUpcdbXZ7XTId5Q88ITzxwNddXZFpdUsDQaYsBTET42o6VYOC7e+sRhL/btSaoQ33S7uQXx9t44r0GzrQPkpEQzVduruD+bSW6zskSRdiE29fkcvuaXE5dGOCJ9xr44YEmnt7fyG2rc3lwexlbgnxOe9/IJJ9//BANPSM8/sAWDfNF0hG6hYwxPPpaNd/dW88d6/L45qfXEx8dXH9je4Yn+NHBZp491ET30AQV2Yk8eH0Zd20sCLuTvv7UMTDOMwca+dGhZgbGplhfmMJ/u66UO9blB92VqLVdwzz09Pu0D4zz/S9u4aMrdEnc2bTlEiSMMXz/3Xr+4bVqVucl870vbqEgCBZ2Otk6wNMHGnnlWBuTDic3rMjiy9eX8dGKzKAeKQab0Uk7Lx69wA/2NVDfPUJWUgxf2FrCfdcUBcXU2LfPdfGV//yAmCgbj31hM1tK0+f/oDCkgR5k3q7u4ivPuX6w/+9nruKGAByljE85ePVkO88caOJYy0XioyO4Z1MhD1xbyvLsRKvLC2tOp2FvTTc/2NfInvPdRNqEHWtzuX9bCdeUpQfcH9kph5PH3qnjX944z8rcZL7/QHAMZKyigR6EaruGeOTZo9R2DfOpzYX8rztWB8S0vuqOQX58uIWXPrjAwNgU5VkJ3L+thHs2F4bsbItg1tAzwrMHm/hJVQtD43aWZSVw79XF3L2pICDW5j/dNsBXf3qC022D3Lkhn0fvWRd0rUZ/00APUuNTDv79rRoe21NPekI0f3F7JXdvLPD7/OOe4QleOdbGzz5o5dSFQaIjbNy+Npf7ri7iI8syAm7Epz5sdNLOL0+08+PDzRxtvkikTbhxZTZ3byzgplXZfl8XpW9kku/uqeOJ9xpIjY/m73et4ePr9AbPntBAD3KnLgzw1y+d5ETrAKUZ8fzRTRXceVW+T9eB7hoa5zenO/n1qXYO1vfhcBrWFaTwyY0F3LWxgPSEaJ/tW/nW+c4hfnqklZc/uEDX0ARJsZHcsiqHHWtzuWFFlk9PYPePTPL4e/U8ta+R0SkHd28s5H/esYqZ/7sHAAAHdElEQVQ0/XnymAZ6CDDG8MbZLv5l93nOtA+SmRjNJzbkc/fGQtYWLH31wfEpB8dbLrK3ppt3a3o4eWEAY6A8M4Eda3P55MYCKnKSvPTVqEDgcBr21fbwyvE2dp/pZGBsirioCLaWp7O9IovtFZksz0pc8hXMk3Yne85387Ojrbx5tospp5M71uXxxzdX6M/UImighxCn0/D2uS5+esT1yzHpcJKbHMvVZelcU5bO6rwkCtPiyUqMmfMX0eE0dA2N09AzQmPPKNUdg5fuvTnlMETYhE3FqXy0Iovb1+ZSkZ2oLZUwMOVwcqi+j91nOni3pof6nhEAkmIjuaoolQ2FqSzPTqQ0M4GyjASS4yLn/LmYtDtpHxijuW+UY80XOdzYx9GmfkYmHWQkRHPnVfncd00xKzTIF00DPURdHJ3ktVMd7K/r5XBDL52D/3WDqOhIG6lxUURF2IiOtDHlcDIwNsXQuP23Pkd8dATrC1O4qiiNTcWpbFuWoSc3FS19oxyo6+WDlosca7nIuY5BnDPiIdImJMdFkRzrOoE55TBM2B30jkwyM0ZW5iZxdWk6N67MYntFlt4uzgs00MOAMYaWvjHquodp7R+lpX+MwbEpJh1OphyGSJuQEhdFclwUWUkxlGUkUJoZT15KnK6louY1PuWgpW+Uhp4RmvtG6RuZZGBsisFxO4JrABEVYSM7KYbCtDgK0uJYnZdMarz2xr1NL/0PAyJCcYbrjvZKeVtsVAQVOUna8w4i+lpIKaVChAa6UkqFCA10pZQKERroSikVIjTQlVIqRGigK6VUiNBAV0qpEKGBrpRSIcKvV4qKSDfQtMgPzwR6vFiOt2hdC6N1LYzWtTCBWhcsrbYSY8y8d73xa6AvhYhUeXLpq79pXQujdS2M1rUwgVoX+Kc2bbkopVSI0EBXSqkQEUyB/j2rC7gMrWthtK6F0boWJlDrAj/UFjQ9dKWUUlcWTCN0pZRSVxCwgS4i3xSRahE5ISIviUjqZbbbISLnRKRWRL7mh7o+LSKnRcQpIpc9Yy0ijSJyUkSOiYjP7+qxgLr8fbzSRWS3iNS436ZdZjuH+1gdE5FXfFjPFb9+EYkRkefdzx8SkVJf1bLAur4kIt0zjtFDfqrrSRHpEpFTl3leROTf3HWfEJFNAVLXx0RkYMbx+hs/1FQkIm+LyFn37+Ifz7GNb4+XMSYg/wG3AZHux/8I/OMc20QAdUA5EA0cB1b7uK5VQCXwDrDlCts1Apl+PF7z1mXR8foG8DX346/N9X10Pzfsh2M079cP/AHwmPvxvcDzAVLXl4Bv+evnacZ+PwpsAk5d5vmdwGuAANuAQwFS18eAX/r5WOUBm9yPk4Dzc3wffXq8AnaEbox53RgzfTPMg0DhHJtdA9QaY+qNMZPAj4FdPq7rrDHmnC/3sRge1uX34+X+/E+7Hz8N3OXj/V2JJ1//zHp/Ctwsvr9jthXfF48YY/YCfVfYZBfwjHE5CKSKSF4A1OV3xph2Y8xR9+Mh4CxQMGsznx6vgA30Wb6M66/abAVAy4z/t/LhA2gVA7wuIkdE5GGri3Gz4njlGGPawfUDD2RfZrtYEakSkYMi4qvQ9+Trv7SNe0AxAGT4qJ6F1AVwj/tl+k9FpMjHNXkqkH8HPyIix0XkNRFZ488du1t1G4FDs57y6fGy9J6iIvIGkDvHU183xvzcvc3XATvwo7k+xRzvW/K0HU/q8sB1xpg2EckGdotItXtUYWVdfj9eC/g0xe7jVQ68JSInjTF1S61tFk++fp8co3l4ss9fAM8ZYyZE5BFcryJu8nFdnrDieHniKK7L5YdFZCfwMlDhjx2LSCLwIvAnxpjB2U/P8SFeO16WBrox5pYrPS8iDwC/A9xs3A2oWVqBmSOVQqDN13V5+Dna3G+7ROQlXC+rlxToXqjL78dLRDpFJM8Y0+5+adl1mc8xfbzqReQdXKMbbwe6J1//9DatIhIJpOD7l/bz1mWM6Z3x3+/jOq8UCHzyM7VUM4PUGPOqiHxbRDKNMT5d50VEonCF+Y+MMT+bYxOfHq+AbbmIyA7gL4E7jTGjl9nsfaBCRMpEJBrXSSyfzZDwlIgkiEjS9GNcJ3jnPBvvZ1Ycr1eAB9yPHwA+9EpCRNJEJMb9OBO4Djjjg1o8+fpn1vsp4K3LDCb8WtesPuuduPqzgeAV4Ivu2RvbgIHpFpuVRCR3+tyHiFyDK+t6r/xRS96nAE8AZ40x/3yZzXx7vPx5FniBZ4xrcfWajrn/Tc88yAdenXXW+Dyu0dzX/VDXJ3H9lZ0AOoHfzK4L12yF4+5/pwOlLouOVwbwJlDjfpvufv8W4HH342uBk+7jdRJ40If1fOjrB/4O18ABIBZ4wf3zdxgo9/Ux8rCuf3D/LB0H3gZW+qmu54B2YMr98/Ug8AjwiPt5Af7DXfdJrjDzy891/eGM43UQuNYPNV2Pq31yYkZu7fTn8dIrRZVSKkQEbMtFKaXUwmigK6VUiNBAV0qpEKGBrpRSIUIDXSmlQoQGulJKhQgNdKWUChEa6EopFSL+P8SK5ycmrLNfAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "x = np.linspace(-2, 2, 100)\n", "y = x*x*x*x - 2*x*x + 8\n", "plt.plot(x, y)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "上面画出了一元函数的图像：\n", "\n", "$$\n", "y = x^4 - 2x^2 + 8\n", "$$\n", "\n", "从图中可以很直观的看出，当 x = 0 时，函数有极大值 8，当 x = 1 或 x = -1 时，函数有极小值 7。这里引申出两个数学概念：**极值** 和 **最值**。在数学分析中，某个函数在指定范围内的最大值和最小值被统称为极值，一般通过 **导数** 来求极值和最值：\n", "\n", "$$\n", "y' = 4x^3 -4x\n", "$$\n", "\n", "令 $ y' = 0 $，可以求得：\n", "\n", "$$\n", "x^3 - x = 0 \\to x = 0, x = \\pm 1\n", "$$\n", "\n", "根据极值点将函数的定义域划分为几个区间，并计算各个区间的单调性：\n", "\n", "|x |($-\\infty$, -1)|-1 |(-1, 0) |0 |(0,1) |1 |(1,$+\\infty$)|\n", "|- |---------------|-----|----------|-----|----------|-----|-------------|\n", "|y'|- |0 |+ |0 |- |0 |+ |\n", "|y |$\\searrow$ |7 |$\\nearrow$|8 |$\\searrow$|7 |$\\nearrow$ |\n", "\n", "检查极值点两侧导数的符号，如果是左正右负，则为极大值点，如果左负右正，则为极小值点。显然，由于函数在 ($-\\infty$, -1) 区间内单调递减，在 (1,$+\\infty$) 区间内单调递增，所以比较所有的极小值，最小的那个就是函数在整个定义域内的最小值。\n", "\n", "这里介绍的是一元函数的极值求解方法，再回到我们刚刚的损失函数：\n", "\n", "$$\n", "F(a, b) = n\\bar{y^2} + a^2n\\bar{x^2} + nb^2 + 2abn\\bar{x} - 2an\\bar{xy} - 2bn\\bar{y}\n", "$$\n", "\n", "这是一个关于 a、b 变量的二元函数，类似的，二元函数的极值也可以通过导数来求解，不过这里要用到 **偏导** 的概念，也就是只对其中一个变量求导，其余变量当作常数看待，那么二元函数就有两个偏导，并令其等于 0：\n", "\n", "$$\n", "\\left\\{\n", "\\begin{aligned}\n", "\\frac{\\partial}{\\partial a} F(a, b) &=& 2n\\bar{x^2}a + 2bn\\bar{x} - 2n\\bar{xy} &=& 0 \\\\\n", "\\frac{\\partial}{\\partial b} F(a, b) &=& 2nb + 2an\\bar{x} - 2n\\bar{y} &=& 0 \\\\\n", "\\end{aligned}\n", "\\right.\n", "$$\n", "\n", "结合两个方程可以求得：\n", "\n", "$$\n", "\\left\\{\n", "\\begin{aligned}\n", "a &=& \\frac{\\bar{xy} - \\bar{x}\\bar{y}}{\\bar{x^2} - (\\bar{x})^2} \\\\\n", "b &=& \\bar{y} - a\\bar{x}\n", "\\end{aligned}\n", "\\right.\n", "$$\n", "\n", "通过这个公式，就可以求出 a，b 的解，这种严格按照公式推导计算出来的解叫做 **闭式解**（closed-form solution）或 **解析解**（analytical solution）。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "bar_x = 45.546666666666674\n", "bar_y = 319.3333333333333\n", "bar_xx = 2096.697333333333\n", "bar_yy = 107367.66666666667\n", "bar_xy = 14876.089999999997\n", "a = 14.93439993913547\n", "b = -360.878802561157\n" ] } ], "source": [ "import numpy as np\n", "\n", "X = np.array([39.93, 42.05, 43.18, 44.68, 49.87, 53.57])\n", "Y = np.array([199, 290, 298, 310, 399, 420])\n", "\n", "bar_x = np.sum(X) / X.size\n", "bar_y = np.sum(Y) / Y.size\n", "print(\"bar_x = {0}\".format(bar_x))\n", "print(\"bar_y = {0}\".format(bar_y))\n", "\n", "bar_xx = np.sum(np.multiply(X, X)) / X.size\n", "bar_yy = np.sum(np.multiply(Y, Y)) / X.size\n", "bar_xy = np.sum(np.multiply(X, Y)) / X.size\n", "print(\"bar_xx = {0}\".format(bar_xx))\n", "print(\"bar_yy = {0}\".format(bar_yy))\n", "print(\"bar_xy = {0}\".format(bar_xy))\n", "\n", "a = (bar_xy - bar_x * bar_y) / (bar_xx - bar_x * bar_x)\n", "b = bar_y - a * bar_x\n", "print(\"a = {0}\".format(a))\n", "print(\"b = {0}\".format(b))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "所以，求得一元线性回归的解为：\n", "\n", "$$\n", "y = ax + b = 14.93x - 360.88\n", "$$\n", "\n", "这种通过求解平方损失函数的最小值来求解线性回归的方法，我们称之为 **最小二乘法（Least Squares，简称 LS）**，这是机器学习中最基础的一个算法。\n", "\n", "通常我们会把损失函数写成下面的形式：\n", "\n", "$$\n", "J(a,b) = \\frac{1}{2m} \\sum_{i=1}^{m}(h(x^{(i)}) - y^{(i)})^2\n", "$$\n", "\n", "这个式子也叫做 **均方误差**（mean-square error, MSE），在统计学中，$h(x^{(i)}) - y^{(i)}$ 被称为**残差**，所以这个式子又叫 **残差平方和**（Sum of Squares due to Error，SSE）。其实它就是我们常用的**欧几里得距离**（欧式距离，Euclidean Distance），也就是残差的 $L_2$ 范数。让损失函数的值最小化，求参数 a b 的值可以写成：\n", "\n", "$$\n", "(a^*, b^*) = arg min_{a,b} J(a,b)\n", "$$\n", "\n", "求解 a 和 b 使损失函数最小化的过程，称为线性回归模型的最小二乘**参数估计**（parameter estimation）。\n", "\n", "实际上，上面的损失函数就是三维空间里的一个曲面（图片来自 Andrew Ng 的在线课程）：\n", "\n", "![](../images/cost-function.png)\n", "\n", "损失函数的最小值就是该曲面的最低点。" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }