notes9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 8: Discriptive statistics, Estimation\n",
    "\n",
    "\n",
    " #### [Back to main page](https://petrosyan.page/fall2020math3215)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Box plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we display the blox plot of the following data:\n",
    "\n",
    "$$\\{-30,-1, -5, -0.5, 0.5, 0.6, 0, 2, 3, 4.6, 4, 7, 18, 35\\}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAD4CAYAAABfVMQ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXQc1b3g8d+t7lZLaqu1y2qZLQECDqsX7LA8QvII2GDjTCbLg4SEZWwThxyYAWZ4IRObPRNIeDySPGwDcZ6ZBAIM2McJyXNwHI5fgh0ZHGPCZpYEVCWptbTUbi1Wd935Q61YmNZml6pa1d/POXWsrltVv7vI5kfVrdtKay0AAABOMbyuAAAA8BeSCwAA4CiSCwAA4CiSCwAA4CiSCwAA4KjgGOW8SgIAOJjyugLIb9y5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjiK5AAAAjgp6XQEAACbbjBkyU0TKvK5HAUg2NcmrJBcAgEJQJiIJrytRACpEeCwCAAAcRnIBAAAcRXIBAAAcRXIBAAAcRXIBAMA4mWbxLzs6lszzuh75jrdFAAAYp4aGvou9rsNYLKvmHq27zhSxS0WC8WBw5kN1dbueEBHp7X28PJG45i6tu84WCXaGQrO+X1u7fdNI15ro8UO4cwEAgI+Ull66urZ296caGjKzI5FvfD2d3nN9Z+cVJ4mIJBIrVooY+6urt5xVXHzJjQMDjbd2dX3zuJGuNdHjh7iSXKxatcqNMHkb34+87FPGE4CTTLNoS0vLCVebZslG0zR2NTfX3ZlKram2rPK1phl40bLK1vX2PhUdOra9/eIzDzr3qsFzAzstq/q+gYHdRd61RqS8/IG9odDMgcFPhhZROp1+46j+/m0lWndcUFp69f3h8Hk9VVVP7jSMyuf6+jYsyXWdiR4/nNJaj1Y+auF4KaVkjDiTyuv4fuRlnzKegOeU1xWYqBkzZJ6MsIiWaRZtUSrUFo3e93WtU4Hu7ps2iISaS0sv+1YkcuPetrZ5DxlG7Z+mT3/7h6ZZtCUc/swt1dW//OPQuSLBjmj07hWBQH1/Z+fXHgsGP/bTurrdjw1d37IqVmudnJMrtlJlO2OxxHKn29vcHFtp2y2fE9HFIsV/qa5+9su9veuP7un5yeMNDfapQ8e1tJxwlW23zIvFEtccfI1E4uqZEzk+q6KpSXa4Mudi5cqVboTJ2/h+5GWfMp4AnBYMnrw+ElnWLiKSTN7WqFS4vaLi4VdFRAKBozZnMuaZI50bCp3879OmXdcqItLVdeMW226dObx8MpKHsdTXW7fadvvticTyWQMDL80LBD6y37YTpSJGcvhxSpXsE8lEcl1joscP50py4fVtbK/j+xGPRQD4iWHUtB34FOhTqqT9wOdQv0h6xP+gGkZtfOhnpYJ9tt1b51S9LKtsvdb7cr6dolRkZyy277KR61VtV1U9ubO5uf6Sjo6FlxYVndkoYk/74FF9EZFAKvf5FT0TOX44V5KLhoYGMU3TjVB5Gd+PvOxTxhPAVGJZ5Wu1Ts7NVaZUWWMs1rV0pHNjseTlh18DHbDtrqNKSq58sqfnkUB3981HR6Pf/auIiG23n2gY0TdznVVScuW7Ezl+OFeSC8uy3AiTt/H9yMs+ZTwBTCWjJQ9OS6UerurpeegT0ejdW0OhU/o6O798lm3HFxUVnXtDOHxOr1JVm1OptdeFw4tuSaXun2nbHedHIiu+lOtaEz1+ONa5AADAJ5QK6HT65cva2//xNhExREJNweCpd9XUbH1ORKSi4serEolr7mpvP/ePIoFEKDR3ZXn5A3uHzres8rWGUds4ffre1eM5fsR6uPG2yJw5c2Tnzp1OXGpKxvcjL/uU8QQ856u3ReCoiqYm2eFKcgEA8BWSC4ykoqlJdriyiNayZcvcCJO38f3Iyz5lPAEgv7GIFg4Ji2gBBY07FxiJe3cuAABA4SC5AAAAjnIluWhqanIjTN7G9yMv+5TxBID85kpy4fVrg17H9yMv+5TxBFCIWltP/rJplj5lmsae5uba7451vGWVrTdN42XTNF4a3MK/dqOeIkzoxCFiQidQ0JjQ6YG2tnM/M/j16a+fI5Iprq+P3zza8ZZVtj4QOHZjXd2uJ9yqozChEwAAEdtOqpaW45ebZmiraQZfaG095TLTVK/09Kyr9Lpuw9XUPL+5pub3v1WqOO+TJJb/BgAUtHj81Gttu/WssrJVl4VCp3V3dn7pEZFAorT0is6xzrWsitVaJ+fkKlOqbKcXX7c+XDq95wbTDN6gVPidcPj8+6qqNuxwI64rycXq1avdCJO38f3Iyz5lPAE4padnXWUm87cry8pWLi4ru8UUETGM+q22HT+jr2/TtI6OL64T6Tu2tPTrX6yo+NGHvg3U6+RhNMXFi+8pLb3qLcOYsT+RuPzivr5ND3Z3f3tJNHrHe5Mdm+W/AQAT5Zs5F/H4mRel03++Mhbr+cLQvubmGd8SUXZNzX/eu3//78q6u2/6X8XFX3w4V3JxOCyrbL3W++blKlMqsjMW23dZrrKWlmOu1zpVP9aciw/Hiz4UCBy5ta7ulUcPpb7jVNHUJDtcuXPh9QQ8r+P7ERM6AfiB1slKkaKOoc+ZjBWw7finQ6GTfxwMHp0OBq/o7O6+acTzLat8rdbJubnKlCprHO3r1mOx5OWHVfkJU1pEu5IYMucCAFCwAoGj306nX70umVx1RDB4SndX14qbRAaODAZPGtdditGSB6dlE5+AiG2IaGNgYHeRYdRmAoFY5uBje3ufKevpefC0aPR7OwyjNtPR8dmLtE6eEQ5fdJcbdSW5AAAUrOrqX/6xubn+18nk7RtFgvFg8NhHbbvVjkSudfQRiBPa2s5akcm8e+3Q53j8tCWBwDE/nD79nQdEBu+iGEZt4/Tpe1drnQj29z9/fTw+66MiYitV/HY4vGBFefm977hRV1eSi0WLFrkRJm/j+5GXfcp4AnBSfX3zd0TkOyIi7e0Lz06n33qvqGh+n8fV+pBsEvHASOXD76KUll7RWVp6xeddqVgOTOgEAEyUbyZ0Hqy19eNfzWSsebFY57UiIpZVvkbrnplKFZnB4ImP1dbufHrSKzu1ubeI1uLFi90Ik7fx/cjLPmU8AUwW2+441jAq3hj6HIt1LWtoGPiHWCz1JRKL8WP5bxwS3hYBCppv71zgsLH8t1uqqqpEKZVzk1XlI5aNtlVVVXndLAAAciK5cEFnZ6dorXNuIjJi2WhbZ+eYq9ICAOAJV5ILr29hex3fj7zsU8YTAPKbK8nFmjVr3AiTt/H9yMs+ZTwBIL8xodPr+KvKRVZ1OXtNFzChEyhoTOjESJjQCQDAZGpurv1uS8sx13tdD7ex/DcAAHmutfXkL6fTb39OpO8Ew6jedPA3oo5VPtzgt7GmTheR9OCeUEtDQ/8CJ+vrSnKxceNGN8LkbXw/8rJPGU8AhcYwqlqLiqr/LZ1+/RyRTPFEyw8WDJ56W13dricmp7YuTeicM2eOG2HyNr4fedmnjCcAJ5lm0ZaWlhOuNs2SjaZp7GpurrszlVpTbVnla00z8KJlla3r7X0qKiKyb9/9dZZV+YBpBl8wzaLnWls//oGvTU8krp5pmiVPD55XfZ/WdtiJOtbUPL+5pub3v1WqOOe8kbHK3eZKcjFjxgw3wuRtfD/ysk8ZTwBOs+33Lywvv//KaPT7F9h2x6e7uq57qKTk8z+ord0zX0RUd/dNX7XtpOru/ucHDaPytbq61/5h2rSbv5ZO772ivX3BOSIi6fTeUE/P+h8Hg8dsqKt7e14oNPPXWndckCueZVWsNs1AY67NsipWT3Z70+k9N5hm8AXLivy8o2PJPKevz5wLAEDBCwZPXh+JLGsXEUkmb2tUKtxeUfHwqyIigcBRmzMZ88yuruWniAxUTZ/+9o9ERKLR297v7f35LwYGXrxYRLZ1d990mogO1tbuWqdUWGpqtv3GsiIv54oXiyWWu9a4gxQXL76ntPSqtwxjxv5E4vKL+/o2Pdjd/e0l0egd7zkVg+QCAFDwDKOm7cCnQJ9SJe0HPof6RdKRdPrdGSLpOtMMNA4/ValpjSIimUxznUioRanhT0KKzcmt+cRVVv5s99DPtbWNz1hWdFFf39OfjEbveNSpGK4kF0uXLh37IB/H9yMv+5TxBOCFQOAoa2Bg5/sNDf05H3UEAtPjAwMD07XulwMJRl9MpOxvBx9rWeVrtU7OzXUdpcoaY7EuF/+hU1pEO7p2CSt0jkGp/F8rxos6skIngEJTUfGj3SKBfS0tH1u6f//OsG23G4nEN47v7Lz8FBGRaPTeXSKSicdnfTWTsQJtbed+RuueU3NdKxbrWtrQYM/KteVKLDIZKzAwsLtIxDZEtDEwsLsok7EC4y0f0tv7TFl7+4Jzhsrj8fmLtU6eEQ5ftM3BruJtERwa3hYBUGgMo9qORu++xrbbZ7a1zd/S3Fy3vbd33Z223TFNRCQYPG6gpOTSa9Ppdz7X0nLEnwYGXrlIqar/cCJ2W9tZK+Lx017OZN5bbtvtS+Lx015uaztrxXjKLat8bUvLcctFRLROBPv7n78+Hp/1QkvLEdvT6T2Xh8MLVpSX3/uOE/UcwvLfk3juuK7hwPLfXvQvy38DBS3/b+kehOW/XcPy3wAAwHmuJBexWMyNMHkb34+87FPGEwDymyvJhWl6+yaO1/H9yMs+ZTwBIL+5klysWrXKjTB5G9+PvOxTxhMA8hsTOifx3HFdgwmdUyo2ABFhQidGxoROAADgPJb/BgDARZZVc4/WXWeK2KUiwXgwOPOhoa8/t6yy9VqnTheR9ODRoZaGhv4FI12rt/fx8kTimru07jpbJNgZCs36fm3t9k2uNGQUrty5aGxsHPsgH8f3Iy/7lPEEMJWVll66urZ296caGjKzI5FvfD2d3nN9Z+cVJw2VB4On3nZgxc6REwsRkURixUoRY3919ZaziosvuXFgoPHWrq5vHjf5rRgdj0UAwOcymYxs2rRJbr/9dtm0aZNkMhmvq5RXbDupWlqOX26aoa2mGXyhtfWUy0xTvdLTs65yMuKVlz+wNxSaOTD4ydAiSqfTbxw10ev0928r0brjgtLSq+8Ph8/rqap6cqdhVD7X17dhidN1nihXHovMnTvX0wl4Xsf3Iy/7lPEExi+TyciFF14o27dvl1QqJZFIRObPny+/+c1vJBD40FdPFKR4/NRrbbv1rLKyVZeFQqd1d3Z+6RGRQKK09IrOsc61rIrVWidzfieBUmU7R/pq9ebm2ErbbvmciC4WKf5LNHrX74fK0uk9N5hm8Aalwu+Ew+ffV1W1YUeua/T2/uQYEWWXl3/v3QMxq1+z7ZZ5Y9V7sjHnAgB87Nlnn5Xt27fLvn37RERk3759sn37dnn22Wdl0aJFHtfOez096yozmb9dWVa2cnFZ2S2miIhh1G+17fgZqdSa6u7u//4jEZUWUZlo9O4bI5Fr48PPHyl5GEt9vXWrbbffnkgsnzUw8NK8QOAj+0VEiosX31NaetVbhjFjfyJx+cV9fZse7O7+9pJo9I73Dr6GbSdKRYzk8H1KlewTyUQOpU5O4rHIOCilDmub7HoBwEheeuklSaVSH9iXSqVk165dHtUov6RSq89UKvxWWdl3mob2ad1XoVT09ZKS/9o5fXrzpbHYvq8Egx97Zt+++z7vZGzDqLarqp7cqXWqvqNj4aUiIpWVP9sdDp+fCoVmDtTWNj6jVOTFvr6nP5n7/IoeEXvaB/f2RUQCqVzHu8mVOxcrV650I8ykxXdinYvJMHydC7d5OaZe/z4BU8msWbMkEon8/c6FiEgkEpHTTz/dw1rlD62TlSJFHUOfMxkrYNvxT4dCJ//YMKrtA8f1RQKBGW8efL5lla/VOjk317WVKmvM9fXpOWoRsO2uEeZcKC2ic/4jX1Jy5bs9PY8EurtvPjoa/e5fRURsu/1Ew4h+qJ5uY4VOHBJW6ASmhoULF8r8+fNl2rRpopSSadOmyfz582XhwoVeVy0vBAJHv611clYyueqI3t6novH46atEBo4MBk96U0QkkfhvJ1pW6RPp9N6vhMMX/uXg82OxrqUH3uz44JYrsUilHq6Kx8+8qL9/a6lttxvt7QvOse34omDw+Bd6e58pa29fcM7AwO6iTMYKxOPzF2udPCMcvmhbrrqHw+f0KlW1OZVae11//7aSjo4vzLbtjvOLi5dscLyjJsiVFTobGho8/T6Iw4nPCp25eTmmXv8+AVNNJpORZ599Vnbt2iWnn366LFy48HAnc06557GjrdDZ3Fx/m23HFw2uOXHso+n0q9+qqXlhVlHR/L6hY+LxsxZmMm9/or6++bBunfb0rKvs6rr2Aa17TxQRQyTUFAyeuL6ubtcvenrWVSYSK9aK9H9URGylit8uKvrkv1RX/+oPQ+dbVvlaw6htnD5972qRg9e5CCRCodn3erzORUVTk+xg+e9JPHdc15iiyQXLfwMFzVfJxXDt7QvP7u/fsrKhof+CgYFXQ0OvjA7eUdhzTn39+9+d9MpObRVNTbKDt0UAAMjKZP56rFKlb4iIJJO3ntTfv+kmEZURMfqj0Tu/5XX9pgpXkovZs2e7ESZv4/uRl33KeAKYLLbdcaxhVLwhIlJV9dguEfmyx1Waklx5LDKV8VgEAD7Et49FcNjc+1bUZcuWuREmb+P7kZd9yngCQH5jQucknjuua0zROxdM6AQKGncuMBL37lwAAIDCQXIBAAAc5Upy0dTUNPZBPo7vR172KeMJoJBYVtl60zReNk3jpcEt/Ovh5b29j5dbVuWPTNPYZZpFv4vH5y+aSPmhHjsaV5KLnTt3uhEmb+P7kZd9yngCKDTB4Km3HVhavH/B8LJEYsVKEWN/dfWWs4qLL7lxYKDx1q6ubx433vKJXGu8XEkuLrnkEjfCTEr8qTBx0Is6ejmmXv8+AfAX0yza0tJywtWmWbLRNI1dzc11d6ZSa6otq3ytaQZetKyydb29T0VFRFpajl9mmkW/Nc3Ai6ZZ/Ku2tnPPH7pOd/e3jzTNwI7Ozis/LiKyb9/9daYZfKGjY8m8yap7f/+2Eq07Ligtvfr+cPi8nqqqJ3caRuVzfX0bloynfCLXmgjmXAAACp5tv39hefn9V0aj37/Atjs+3dV13UMlJZ//QW3tnvkiorq7b/qqiIhh1P0tGr3nsvr6xJxQaPYP9+/fdm8q9cNaEZFo9I73gsFT7u3t/fn39+/fXpxM/u+7DaP+6aqqDTuGx7KsitWmGWjMtVlWxepc9Uun99xgmsEXLCvy8+HJSm/vT44RUXZ5+ffeHdqnVPVrtt19/HjKh5vIsWNh+W8AQMELBk9eH4ksaxcRSSZva1Qq3F5R8fCrIiKBwFGbMxnzTBGR2tr//Pt8h9raP/zKNEuW9/Y+eWokcu1zIiJ1dbt+YVkVn2pr++QTIkrX1Gy+5uBYsVhi+UTqVly8+J7S0qveMowZ+xOJyy/u69v0YHf3t5dEo3e8Z9uJUhEjOfx4pUr2iWQiIiJjlQ83kWPH4sqdi9WrcyZirvE6vh952aeMJwCnGUZN24FPgT6lStoPfA71i6QjIiLx+NzPmmbJhqE7DSJ9x9t2V+Xwa4VCs38h0v+xYPCER4e++OxwVFb+bHc4fH4qFJo5UFvb+IxSkRf7+p7+5GC9K3pE7GkfPKMvIhJIjad8uIkcOxZW6MQhYYVOAIUmmbyzYWDgxTuKiy++rb6+dV5DQ2auSPGbIvrvi4r1928t3b9/2y2GUfdEOv3KN3t7Hy8/+DqDczmG3vz44GZZ5WvHronSQzFLSq58V0QHurtvPnqo1LbbTzSM6JvjKR9uIseOxZXkQilvF3PzOr4fedmnjCcAL9h2S4mI6EDgmA4RkXh89udE+j4wH6Gz84u3KDVtT319y7cNo2prIvH1Ww++TizWtfTAmx8f3GKxrqXDj+3tfaZs8OvedxdlMlYgHp+/WOvkGeHwRdtERMLhc3qVqtqcSq29rr9/W0lHxxdm23bH+cXFSzaMp3y4iRw7FiZ0AgAwDuXl//pWIHDkI6nUfY83N0//QybT+jGlIi8Olbe1nfePtt15bkXFv60UEamsfPxurVMnxePzFx9qTK0Twf7+56+Px2e90NJyxPZ0es/l4fCCFeXl974zdExFxY9Xidjh9vZz/9jX98wPQqG5K8vLH9g7nnLLKl/b0nLc8vFea7z4bhGv4zvw3SJe4LtFgII25W4f8t0irnHvu0UWLTqkBb58E9+PvOxTxhMA8psrdy4KnR/vXAAoaNy5wEjcu3OxePEhP27yRXw/8rJPGU8AyG/MufA6/hS9c8GcC6CgcecCI3HvzgUAACgcJBcuUUrl3EYrG22rrKwcIyIAAN5w5btFvL6Fne/x9Sp36uEkL/vU6/EEMCUlRaTC60oUgKSIS3Mu1qxZ4+mSzV7H9yMv+5TxBDw35eZcwF1M6MQhYUInUNBILjAq5lwAAABHkVwAAABHuZJcbNy40Y0weRvfj7zsU8YTAPKbK8nFnDlz3AiTt/H9yMs+ZTwBIL8xoROHhAmdQEFjQidGxZwLAADgKJILAADgKFeSi6VLl7oRJm/j+5GXfcp4AkB+c2XOBQDAV5hzgVHxtggOCW+LAABGwtsiOCS8LQIUNO5cYFRM6AQAAI5yJbmIxWJuhMnb+H7kZZ8yngCQ35jQCQCYKB6LYFSu3LlYtWqVG2HyNr4fedmnjCcA5DcmdOKQMKETKGjcucComNAJAAAcRXIBAAAc5Upy0djY6EaYvI3vR172KeMJAPmNOxcAAMBRTOjEIWFCJ1DQmNCJUXHnAgAAOIrkAgAAOCowxoJEoxZOxHnnnefUpaZkfD/ysk8ZT8BTt3pdAeQ3lv8GAEwUcy4wKh6LAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAAR5FcAAAARymt9ciFSv1aRGomMX6NiLRN4vXzRaG0U6Rw2ko7/YV2Tkyb1nqBA9eBT42aXEx6cKUatdZzPauASwqlnSKF01ba6S+0E3AWj0UAAICjSC4AAICjvE4u1ngc3y2F0k6Rwmkr7fQX2gk4yNM5FwAAwH+8vnMBAAB8huQCAAA4yrXkQin1BaXUK0opWyk196Cyf1ZK7VVKva6UunDY/jlKqZezZf+qlFJu1dcpSqkF2XbtVUrd7HV9DodS6hGlVKtSas+wfVVKqc1KqTezf1YOK8s5rvlOKXWkUup3SqlXs7+z12X3+6qtSqlipdQOpdSfs+28NbvfV+0copQKKKVeUkptyn72XTuVUu9m/83cpZRqzO7zXTsxBWitXdlEZKaInCAiW0Vk7rD9HxeRP4tIWEQ+IiJviUggW7ZDRM4UESUiz4rIQrfq61CbA9n2fFREirLt/LjX9TqM9pwrIrNFZM+wfd8TkZuzP98sIv9nrHHN901EYiIyO/tzmYi8kW2Pr9qa/Xs1LftzSES2i8gn/NbOYe39HyLyMxHZlP3su3aKyLsiUnPQPt+1ky3/N9fuXGitX9Vav56jaImIPKa17tdavyMie0VknlIqJiJRrfUftdZaRP5dRD7rVn0dMk9E9mqt39Za7xeRx2SwvVOS1vp5Eek4aPcSEflp9uefyoExyjmurlT0MGmtLa31i9mfkyLyqojMEJ+1VQ/al/0Yym5afNZOERGl1BEicrGIPDRst+/aOYJCaSfySD7MuZghIu8N+/x+dt+M7M8H759KRmqbn0zXWlsig/9RFpG67H5ftF0pdYyIzJLB/6v3XVuzjwp2iUiriGzWWvuynSLyLyLyP0XEHrbPj+3UIvIfSqmdSqll2X1+bCfyXNDJiymlfisi9TmKbtFabxjptBz79Cj7pxI/tOFQTfm2K6WmichTInK91rp7lCk/U7atWuuMiJyulKoQkaeVUiePcviUbKdSapGItGqtdyqlzhvPKTn25X07s87WWptKqToR2ayUem2UY6dyO5HnHE0utNbnH8Jp74vIkcM+HyEiZnb/ETn2TyUjtc1PWpRSMa21lX2U1ZrdP6XbrpQKyWBi8X+11v8vu9uXbRUR0VonlFJbRWSB+K+dZ4vIJUqpi0SkWESiSqlHxX/tFK21mf2zVSn1tAw+5vBdO5H/8uGxyEYR+SelVFgp9REROV5EdmRv3yWVUp/IviXyVREZ6e5HvvqTiByvlPqIUqpIRP5JBtvrJxtF5GvZn78mB8Yo57h6UL8Jy/6+PSwir2qtfzCsyFdtVUrVZu9YiFKqRETOF5HXxGft1Fr/s9b6CK31MTL4d3CL1vor4rN2KqUiSqmyoZ9F5AIR2SM+ayemCLdmjorIf5HBTLlfRFpE5DfDym6RwZnKr8uwN0JEZK4M/uV4S0R+KNkVRafSJiIXyeDbBm/J4OMhz+t0GG35uYhYIu0W9ZsAAACKSURBVDKQHcurRaRaRJ4TkTezf1aNNa75vonIOTJ4e3i3iOzKbhf5ra0icqqIvJRt5x4R+U52v6/aeVCbz5MDb4v4qp0y+Fban7PbK0P/3vitnWxTY2P5bwAA4Kh8eCwCAAB8hOQCAAA4iuQCAAA4iuQCAAA4iuQCAAA4iuQCAAA4iuQCAAA46v8D3oRK198Z4d0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# set up the figure\n",
    "fig, ax = plt.subplots(1,1, figsize=(7,4))\n",
    "\n",
    "# data\n",
    "data = np.array([-30, -5,  9, 10, 13, 20, 40, 500]).astype(float)\n",
    "# data = np.array([-30,-1, -5, -0.5, 0.5, 0.6, 0, 2, 3, 4.6, 4, 7, 18, 35]).astype(float)\n",
    "\n",
    "\n",
    "def percentile(data, p):\n",
    "    \"\"\"\n",
    "    Compute the percentiles the way we defined in class \n",
    "    data : array of size N \n",
    "    p : percentile\n",
    "    \"\"\"\n",
    "    data = np.sort(data, axis=0)\n",
    "    rank = int(p * (data.shape[0] + 1) - 1) # the rank\n",
    "    assert rank > 0, \"the rank does not exist\" \n",
    "    alpha = p * (data.shape[0] + 1) - 1 - rank # the fractional part\n",
    "    return data[rank] + alpha * (data[rank + 1] - data [rank])\n",
    "\n",
    "def box_plot(ax, data, width=0.4, showout = True, position = np.array([0.4])):\n",
    "    \"\"\"\n",
    "    ax : matplotlib ax\n",
    "    data : the data \n",
    "    width : box width\n",
    "    showout : show the outliers   \n",
    "    position: the y axis of the box plot\n",
    "    \"\"\"\n",
    "    # compute the five number summary \n",
    "    minim = np.min(data)\n",
    "    maxim = np.max(data)\n",
    "    q1 = percentile(data, 0.25)\n",
    "    q2 = np.median(data)\n",
    "    q3 = percentile(data, 0.75)\n",
    "\n",
    "    # interquartile range\n",
    "    iqr = q3 - q1\n",
    "\n",
    "    # inner fences\n",
    "    left_innerfence = q1 - 1.5 * iqr\n",
    "    right_innerfence = q3 + 1.5 * iqr\n",
    "\n",
    "    # compute outliers \n",
    "    outliers = data[np.logical_or(data <left_innerfence, data >= right_innerfence)]\n",
    "    \n",
    "    # whiskers\n",
    "    if showout==True:\n",
    "        low_whisker = np.min(data[data >= left_innerfence])\n",
    "        high_whisker = np.max(data[data <= right_innerfence])\n",
    "    else:\n",
    "        low_whisker = np.min(data)\n",
    "        high_whisker = np.max(data)\n",
    "\n",
    "\n",
    "\n",
    "    stats = [{'iqr': iqr,\n",
    "              'whishi': high_whisker,\n",
    "              'whislo': low_whisker,\n",
    "              'fliers': outliers,\n",
    "              'q1': q1,\n",
    "              'med': q2,\n",
    "              'q3': q3}]\n",
    "\n",
    "    # add the box plot\n",
    "    flierprops = dict(markerfacecolor='black', markersize=5)\n",
    "    ax.bxp(stats, vert = False, widths=width, positions = position, \n",
    "           flierprops=flierprops, showfliers=showout)\n",
    "\n",
    "    # add Tukey's fences\n",
    "    if showout==True:\n",
    "        ax.vlines(q1-1.5*iqr, position-0.2,position+0.2, linestyle=\"dashed\", linewidth=1)\n",
    "        ax.vlines(q3+1.5*iqr, position-0.2,position+0.2, linestyle=\"dashed\", linewidth=1)\n",
    "\n",
    "        ax.vlines(q1-3*iqr, position-0.2,position+0.2, linestyle=\"dashed\", linewidth=1)\n",
    "        ax.vlines(q3+3*iqr, position-0.2,position+0.2, linestyle=\"dashed\", linewidth=1)\n",
    "\n",
    "    #  \n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    ax.set_ylim(-0.1,position+0.3)\n",
    "    ax.set_yticks([])\n",
    "    plt.figtext(1,0.8,\n",
    "                r\"$\\min={:.4}$\".format(minim)+\"\\n\"+\n",
    "                r\"$q_1={:.4}$\".format(q1)+\"\\n\"+\n",
    "                r\"med$={:.4}$\".format(q2)+\"\\n\"+\n",
    "                r\"$q_3={:.4}$\".format(q3)+\"\\n\"+\n",
    "                r\"max$={:.4}$\".format(maxim),\n",
    "                ha=\"left\", va=\"top\",\n",
    "                backgroundcolor=(0.1, 0.1, 1, 0.15),\n",
    "                fontsize=\"large\")\n",
    "    \n",
    "def disp_data(ax, data):\n",
    "    ax.scatter(data, np.zeros(data.shape), zorder=2, s=10)\n",
    "    ax.set_yticks([])\n",
    "#     ax.set_xticks([])\n",
    "    mean = np.mean(data)\n",
    "    ax.scatter(mean, 0, zorder=2, s=20, color=\"red\")\n",
    "    ax.set_ylim(-0.01,0.1)\n",
    "    ax.axhline(y=0, color='k', zorder=1, linewidth=0.5)\n",
    "\n",
    "    \n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "\n",
    "    ax.set_ylim(-0.1,1.5)\n",
    "            \n",
    "box_plot(ax, data, width=0.2, showout=True)\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ordinary least squares regression\n",
    "\n",
    "The following data is taken from [here](https://people.sc.fsu.edu/~jburkardt/datasets/regression/x03.txt). The dataset contains systolic blood pressure for 30 people of different ages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:61% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Age</th>\n",
       "      <td>39</td>\n",
       "      <td>47</td>\n",
       "      <td>45</td>\n",
       "      <td>47</td>\n",
       "      <td>65</td>\n",
       "      <td>46</td>\n",
       "      <td>67</td>\n",
       "      <td>42</td>\n",
       "      <td>67</td>\n",
       "      <td>56</td>\n",
       "      <td>64</td>\n",
       "      <td>56</td>\n",
       "      <td>59</td>\n",
       "      <td>34</td>\n",
       "      <td>42</td>\n",
       "      <td>48</td>\n",
       "      <td>45</td>\n",
       "      <td>17</td>\n",
       "      <td>20</td>\n",
       "      <td>19</td>\n",
       "      <td>36</td>\n",
       "      <td>50</td>\n",
       "      <td>39</td>\n",
       "      <td>21</td>\n",
       "      <td>44</td>\n",
       "      <td>53</td>\n",
       "      <td>63</td>\n",
       "      <td>29</td>\n",
       "      <td>25</td>\n",
       "      <td>69</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SDP</th>\n",
       "      <td>144</td>\n",
       "      <td>220</td>\n",
       "      <td>138</td>\n",
       "      <td>145</td>\n",
       "      <td>162</td>\n",
       "      <td>142</td>\n",
       "      <td>170</td>\n",
       "      <td>124</td>\n",
       "      <td>158</td>\n",
       "      <td>154</td>\n",
       "      <td>162</td>\n",
       "      <td>150</td>\n",
       "      <td>140</td>\n",
       "      <td>110</td>\n",
       "      <td>128</td>\n",
       "      <td>130</td>\n",
       "      <td>135</td>\n",
       "      <td>114</td>\n",
       "      <td>116</td>\n",
       "      <td>124</td>\n",
       "      <td>136</td>\n",
       "      <td>142</td>\n",
       "      <td>120</td>\n",
       "      <td>120</td>\n",
       "      <td>160</td>\n",
       "      <td>158</td>\n",
       "      <td>144</td>\n",
       "      <td>130</td>\n",
       "      <td>125</td>\n",
       "      <td>175</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import pandas as pd \n",
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:61% !important; }</style>\"))\n",
    "\n",
    "pdata=pd.DataFrame.transpose(pd.read_csv('files/blood_pressure.csv').astype(int))\n",
    "display(HTML(pdata.to_html(header=False)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGDCAYAAAAiZbrgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXyU1fXH8c/JSkgCQQhLWFQUEawIGrQVd61o64Jb1ap1Qa11qVq1Sm1dat2XKu5LcavKzyoiWhG3WncxyKJUUHGFSUiQLZCQTJLz++OZ4BCSMEAmk0y+79crL2fu3HmeM0OU473n3mvujoiIiEgySkl0ACIiIiLxokRHREREkpYSHREREUlaSnREREQkaSnRERERkaSlREdERESSlhIdaffM7AQzeyVO137EzP62ie91M9s28vg+M/tLy0bXcZjZVDM7OdFxiEj7Y9pHR9oDM9sDuAnYAagFPgMucPeP4nzfR4CF7v7nTXivA4Pc/csWD0xERGKSlugARDbEzLoALwK/A54GMoA9gapExtXWmZkR/M9MXTN90ty9pgXv2aLXExHZXJq6kvZgOwB3f8rda9290t1fcfc5AGZ2ipm9U985MmV0tpl9YWblZnaNmW1jZu+b2Uoze9rMMiJ99zGzhWb2JzNbYmbfmNkJTQViZoeY2SwzW25m75nZsFg+QPQUWNQ9LzKzUjMrNrNTo/pmmtktZvadmS2OTHtlRV7rZmYvmlmZmS2LPO4X9d43zexaM3sXqAAGNhLLN2Z2qZnNAVabWZqZFZjZs5Hrfm1mv4/qn2Vmj0bu95mZ/dHMFm7G9XY1s6LIn8ViM7st0t7JzP5pZj9Evt+PzKxX1Oc6PfI4xcz+bGbfRr6/x8ysa+S1rSJ//idHvr8lZnZ5LH9GIpKclOhIe/A5UBv5y/ZgM+sWw3sOAnYBfgr8EXgAOAHoD/wEOD6qb2+gB9AXOBl4wMwGN7ygme0MTAB+C3QH7gemmFnmJnym3kDXyD3HAndHfa4bCZK74cC2kT5XRF5LAR4GtgQGAJXAXQ2ufRJwJpALfNvE/Y8HfgnkAXXAC8DsyL32By4ws9GRvlcCWxEkTT8HTtzM690B3OHuXYBtCEbpIPjuuxL8GXUHzop8voZOifzsG4kpp5HvYA9gcOTeV5jZkCa+BxFJckp0pM1z95UEf3E58CBQZmZT6v9vvwk3uvtKd58LfAq84u5fufsKYCowokH/v7h7lbv/F/g38KtGrnkGcL+7fxgZWXqUYPrsp5vwscLAX9097O4vAauAwZHppjOAC919qbuXA9cBx0W+ix/c/Vl3r4i8di2wd4NrP+Luc929xt3DTdx/vLt/7+6VwEgg393/6u7V7v4Vwfd8XKTvr4Dr3H2Zuy8Exm/m9cLAtmbWw91XufsHUe3dgW0j3++MyJ99QycAt0X+PFcB44DjzCx6Kv7qyMjfbIKEa6cmvgcRSXJKdKRdcPfP3P0Ud+9HMCJTANzezFsWRz2ubOR5TtTzZe6+Our5t5HrN7QlcFFkWmW5mS0nGH1orO+G/NCglqUiElM+0BmYEXWPlyPtmFlnM7s/Mm2zEngLyDOz1KhrfR/D/aP7bAkUNPhcfwLqE8mCBv0bu/7GXG8swYjVvMj01CGR9seBacBEMwuZ2U1mlt7IvQpYd6TqW4J6w+jEtyTqcf13KyIdkIqRpd1x93kWrIb6bQtdspuZZUclOwMIRoEa+h641t2vbaH7NmYJQSK2g7svauT1iwimZHZz9xIzGw7MBCyqTyxLKaP7fA987e6DmuhbDPQD/hd53n9zrufuXwDHm1kKcCTwjJl1j3z/VwNXm9lWwEvAfOAfDS4RIkim6g0AagiS2X6IiETRiI60eWa2faRwt1/keX+CmpAPmn/nRrnazDLMbE/gEOBfjfR5EDjLzHazQLaZ/dLMclsqiMgKqQeBv5tZTwAz6xtV35JLkAgtN7MtCOpnNtd0YGWkoDjLzFLN7CdmNjLy+tPAOAsKofsC527O9czsRDPLj3zW5ZH31JrZvma2Y2R0aiXBVFZtI9d/CrjQzLY2sxyCqb3/02ovEWmMEh1pD8qB3YAPzWw1QYLzKcHoRksoAZYRjBQ8AZzl7vMadnL3IoL6mbsi/b8kKIptaZdGrv1BZHrqNYJRHAim67IIRn4+IJjW2izuXgscSlD8/HXk2g8RFAYD/BVYGHntNeAZmlnaH8P1DgLmmtkqgsLk49x9DUGB9jMESc5nwH+BfzZyiwkE01xvRa6/Bjhv4z+5iHQE2jBQOjQz2wf4Z6T2R2JgZr8jSE4aFkGLiLQ5GtERkWaZWR8zGxXZv2YwwUjac4mOS0QkFipGFpENySDYM2hrgpqaicA9CY1IRCRGmroSERGRpKWpKxEREUlaSnREREQkabXrGp2DDjrIX355s1fXiohIcrENd5GOol2P6CxZsiTRIYiIiEgb1q4THREREZHmKNERERGRpKVER0RERJKWEh0RERFJWkp0REREJGkp0REREZGkpURHREREkpYSHREREUlaSnREREQkaSnRERERkaSlREdERESSVtwSHTPrb2b/MbPPzGyumZ0fab/ZzOaZ2Rwze87M8qLeM87MvjSz+WY2Ol6xiUhsJs9cxKgb3mDry/7NqBveYPLMRYkOSURko8RzRKcGuMjdhwA/Bc4xs6HAq8BP3H0Y8DkwDiDy2nHADsBBwD1mlhrH+ESkGZNnLmLcpE9YtLwSBxYtr2TcpE+U7IhIuxK3RMfdi93948jjcuAzoK+7v+LuNZFuHwD9Io8PBya6e5W7fw18Cewar/hEpHk3T5tPZbh2nbbKcC03T5ufoIhERDZeq9TomNlWwAjgwwYvnQZMjTzuC3wf9drCSFvDa51pZkVmVlRWVtbywYoIAKHllRvVLiLSFsU90TGzHOBZ4AJ3XxnVfjnB9NYT9U2NvN3Xa3B/wN0L3b0wPz8/HiGLCFCQl7VR7SIibVFcEx0zSydIcp5w90lR7ScDhwAnuHt9MrMQ6B/19n5AKJ7xiUjTLhk9mKz0dcvkstJTuWT04ARFJCKy8eK56sqAfwCfufttUe0HAZcCh7l7RdRbpgDHmVmmmW0NDAKmxys+EWnemBF9uf7IHembl4UBffOyuP7IHRkzYr0ZZRGRNst+HFBp4Qub7QG8DXwC1EWa/wSMBzKBHyJtH7j7WZH3XE5Qt1NDMNU1lWYUFhZ6UVFRHKIXEZF2rLFSCOmg4pbotAYlOiIi0gglOrKWdkYWERGRpKVER0RERJKWEh0RERFJWkp0REREJGkp0REREZGkpURHREREkpYSHREREUlaaYkOQEREZFN9/jmsWrVu26GHsmtiopE2pnzRIj5ToiMiIu3WqlXQtet6zcsTEIq0PXmgqSsRERFJYkp0REREJGkp0REREZGkpURHREREkpYSHRER6bBKSvJvWLx4qwsSHYfEjxIdERGRduKHHw75WSiUOTUUSpldXJz7WHn5tQVN9V2x4oKBxcW5j4ZCqTNCocxXlyzZ64D618rKdjs0FEqZGfUzOxSy+cuWnbIDgHsVJSUDLg6FUj8MhVI/LCnpf4l71dprL136qxGhUOdnQqHUj0OhrClLlx69S1w/+GZQoiMiItIOVFQ83q2qaupd6emFd+TnzxpplvfpqlXX3d5Y39ra4tTVq++9NzW175u9e5eOzMw88C/V1e/csmLFH7cCyM//8IWCgroR9T/p6SOugvTv8/LunwtQVrbrsXV1iw/o0uXGw7p0ueXQurqyfcvKRh4HUFn5f13XrJl0b3r6iId69y4tTE/f8aE1a567r7Ly2S6t9V1sDCU6IiKS1ILRib7jQqG094PRjawpy5efM6ixvqWlw38VCmW+GgqlTi8uzrt31ao7eta/FgrZ/NLSoSeFQhmvh0JpH5SU9P9jXV251b9eVrbzUaFQp6mhUOpHxcVd/tHcaMumWLXqlp+bdfoiP//dl9PTh1V36/bUne6V269Y8YeBDfuWl/9tIIR75ufPfjglpXtd9+4vfGCW8/GaNU8f3ti1a2o+PyI1td9ks0wAamu/PCItbfsJOTkXL87JubA0LW3IhNraBUcCrF79wAhI+yE//92XU1K61+XnT58CaUvLy/92YEt+3paiREdERJLa0qWH71FXt2xkt27PHNinT8Uu2dm/vSA9fch6mwr+8MOhP62p+fQPWVm/Pj8//9NRZp1D5eWX/z26T23t9z/v1u3xI3Nzrzmirq50/yVL9jwaYMmSvQ8Ih+eelZ191jm9ei38aUpKftGqVdfd1lRMoVBqUVM/ixcPOrOx99TVLRlkljuv/nlm5h6VkPFdOPzRto30tvXbsLq6Fds1bCwvv7bAfdXIzp1/M7m+zX3NoPT0HdbeKz19x3nuayL3cQMaXt/q6kobTR4TTYmOiIgkufQaqM1es+b/BrpXW9eut3+VnX1uWcNe4fCHh6ak9Hm2W7eH/5eePiTcrduTt7qvHl5e/te+9X0yMnZ/MCvr2BW5ueOK09K2e7S2dsEvg/fOOjYtbej9Xbve/lVqap/a/PxZ97lXDmlqVKegoLawqZ9evb54oPHPEe5sllEe3WKWWu5ekd2wZ07OZV9B2tLS0h1Pr6n5Nu2HHw4e5V4+Emo7NexbUTFhjFlOUW7uVQt/bK3rnJLSe+3hGqmpBeVQl+1eRefOp8+EcM+ysp/+sqbm27SyssIxUD0AarMajzuxlOiIiEhS6979hQ/S0rb5Z2Xl81eWlHR7v6Sk1zVr1ry8XnLgXtUzJaVrqP55ZuY+FZC6PBz+pFd9W1radsX1j1NSei1yD/eKvLdvTc3sy+tHZUpK8j4CLBye3YsWk17hHs5ZN+baHLPOqxv2TEvbsqZz59POrqsL7VNaOvDd6up3TzPrPtWsU0nDvrW1C8ekpQ16bt3WlIq6utLsH/uU5EDKarNMOnf+9fJOnQ79XTg859TS0oHv1dR8vZdZ7ntm2etduy3QWVciIpL0evb83+PA46tX/2OLlSsvuGPFirNO79Tpmzui+5hlltbVrVg7AlNV9U4W1Oalp++4uL6tpubzPsCXAHV1iwvM0hcH780oTkvb6d78/A9fiCWeUChlZlOvpaYOvK9Xry/vb9iektLji9raL49YN77qAenpI79s7Dp5effNz8u778T658XF2RNTUrZdJ6FZuvSYnSHcs0uXm6ZFt5t1+iIc/nQI8Enwuedsb9Zp7X222GLyR8DREBQ+L1685WtpaYMmbOhzJ4JGdEREJKktW3bSjsuW/XpYTc23aWlpgyvBqsDqGvZLT9/1xbq64qOWLz99+3D4s/Rly479g1n27NzcKxbV96mufm9sZeWzXcrLb+pdU/P5b1JTB74UvHf4U+Hw7N+uWHHetgBr1ryYU1Y26qCmYope8dTwp7EkByAn58JX3ddst2TJngeGw3Myli379TlmWfO7dr3tq8b6L19+1uBweE5GdfWHnRYvHnyaezi/W7eJk6L7VFe/dURKSvdpmZkHrDMqlJq6zeSamnmnrlr1956rVt3RMxz+32mpqdtM+vHaY4fU1HybtmbNy9llZYWXmqWXdO/+8jtNfd5E0oiOiIgktbq6ZTlVVa+Pq6yc2B9Sqszy3snLe/gfDft17/7i+6Wlw26vqHjsroqKh7uY5c7Mzb3qwug+qan9X1+27PjnoDYnJaXPcz16vPEMQI8eb71WVjYye/XqB/++evU9fSGlPCWl27vAyy31OTp3PmVZZeW/zq2qeu2KsrLht5hlz87JuXRtfIsXb/vburqywj59VpwBsGbNvw+vqHjwGPA0s9yinJyLT01PHxKu7x8Oz8moq1tycGbmwec1vFd+/kcTFy/etv/KlRe/CJCSUvCv/PyPJta/Xlk5+YyKikf2BjDLe6tLl/HntNTnbGnm7omOYZMVFhZ6UVFRosMQEZEE+fhj6Np13ba99mJwPO4VCtn8nJxxP+/S5brv4nF9aXF5ixYxXVNXIiIikrSU6IiIiEjSUo2OiIhIDAoKPC5TYhJfGtERERGRpKVER0REkl4olDJz5cor+iU6Dml9SnRERCTpFRTUjejS5a8LN9wzWF21cuWfBsQ7po1VVrbLEcXFOU/G49orVlwwsLg499Hg0NPMV5cs2euAde+9+8GRA0s/DoU6vbRkyd4HNHWtUChl5ro/9llJScFfAGpqvkwvLt5ifCiU8UYoZPOXLj1818auUVPzZXoolDk1FEp/q76tvPz6Po1ce35p6fanNvfZlOiIiIi0kNra4tREx7CxamuLU1evvvfe1NS+b/buXToyM/PAv1RXv3PLihV/3Apg1aq/9wyH3785M3Ov6/v0qdg5I2OPG6ur37519ep/bNHY9aI3P+ze/a3dwdZkZOw2tf711NReMzp1OuwSSFvvvLF6P/xw0Fiz9KXRbbm544qjr52be8WhQF1W1q9fae7zKdEREZGkFz1KU1KSf0NJSZ8riou7PhAKpX5cXNz5XytX/rk/QHFxzhMAq1bdMCUUSplZVrb7LwCWLNl/n1Ao6/lQKLWouDh74vLlZw3+8doZbyxevN0ZoVDWlMWL+86urS1OLS+/qXdxcbe7QqG0D0Kh1A9LSvpcUd+/rGznoyKjIx8VF3f5R/TBn8EIxdCTQqGM10OhtA9KSvr/sa6u3FasuGBgODzzr+6rRwQjGakttolcefnfBkK4Z37+7IdTUrrXde/+wgdmOR+vWfP04QDh8Ee9IbW8e/dX3jLLpEeP1/4LKZVVVa9ucNRr5coLRkPa0m7dJhYBpKVtG+7Z87NHt9jimRnAertTB/Fc1a+2duHhGRm7N7pDdL2KisfGmOV8FL1zdWOU6IiISIdTV1d6SKdOh97Zq9fCkdD529Wr7/0DQJ8+q04AyMm57LCCgroR+fnvvbRs2alDq6v/e31W1hF/6d27dNe0tCETKyom3BsOf5Zef73a2u8Pyc29+swePT7axSzDy8uvvN+s86Lu3d/cNz9/5p4ZGaP+DbBkyd4HhMNzz8rOPuucXr0W/jQlJb9o1arrbouOrbb2+5936/b4kbm51xxRV1e6/5Ilex7dtevtX6Wnj7jCLHtmMKJRW9jY5yop6XNl/cGi6/9kTWni27BGGq2ubsV2AF273vmpWacFS5bsu19d3Q8pwbSVVWdnnzt/Q99zTc38I1JT+042y9xQ17VWrbr9LxkZu95qlrWmuX7BYaTbPddcH1CiIyIiHVBKSo9XunX75yepqX1q09NHvOC+avum+lZVTf1Vamq/id26PTknJaV7XX5+0WSw6vLyccPr+6SnD3ksN/ePJRkZu1QtX37OMAj3ys8vuikzc4/K9PRh1ZERDMLhWcempQ29v2vX279KTe1Tm58/6z73yiHRozoZGbs/mJV17Irc3HHFaWnbPVpbu+CXsX6u3r2Lry4oqC1s/KfysMbek5Nz2VeQtrS0dMfTa2q+Tfvhh4NHuZePhNpOwXfVvS41ddDk6ur/3lpS0uPT6uq3b83I2PsvmZl7VDYXS3n59X3cV+3aufPJG0xG6gW1QXWpPXq89Vpz/ZYuPWoXqOnepcut05rrB0p0RESkAzLLWlL/OCWlcyXUZTfV172yoLb2u9OiR0cg3Ke2tqTnj9foXVz/uLb2uz6QsSg1tU/t+teq6ltTM/vy+uuUlOR9BFg4PLtXfZ+0tO3WXislpdci93CvhtdpSWlpW9Z07nza2XV1oX1KSwe+W1397mlm3aeadSoB+OGHX/6spmbOJVlZJ57Uu/fKHbKyTjixuvq/1y5ffnqTySFARcVDY8yyZ+TmXhVTEXhV1TtZ1dUf/jE7++xrNtS3uvqdI1JSur+SmblPxQY/Xyw3FxER6ajMskpSUvLv7dXry/ua6bX24MjU1AHF4fD0gtra4tSGyY5ZRnFa2k735ud/+EJTV6qp+bwP8CVAXd3iArP0xQ3v0ZSSkt5X19WVNjpyAxmhgoI1jY4O5eXdNz8v774T658XF2dPTEnZ9rkgni+HmHUp6tbtsU8BunV7/JM1a16cXVX19u7AvKZiqa1dOCY9fYcHNhRzvcrKR7aE6r6rVt3y5KpVtwCeDnW5oVDau7m5V/yqvhanunpGZnAY6S9iOkhUIzoiIiLrSF0SDs/sX/8sM/PAp2trvz1+2bJfD3OvoqrqnawlSw7Yu6rqtUZHgfLy7p4D6aVlZSMvrqp6JyscnpOxdOkxOwOkpw9/Khye/dsVK87bFmDNmhdzyspGHRT9/urq98ZWVj7bpbz8pt41NZ//JjV14EsAKSndf3Cv7l1T82X6+ncN9O5dcmX0yqR1fxpPcgCWLz9rcDg8J6O6+sNOixcPPs09nN+t28RJQcxDPnEv36V+BGf58rFD3FcWpqUNbLJGZ+nSX42AcK8uXW5a7/T2cPiz9HB4TgaAe3V6ODwnw72K3Nxrvuja9Z69u3a98/CuXe88PCNj1OWQtqRr1zsPz84+Z+0o14oVZ/8cUlduscUzHzR1/2ga0REREYmSljb0zqqq124MhVI7pafv9pf8/Pem1tUt/nNl5eQrKiv/bytIWWOWOwNodOVTSkr3utzcK85ateqmP//ww95vgnlKSs8XgY979HjrtbKykdmrVz/499Wr7+kLKeUpKd3eBdYmBKmp/V9ftuz456A2JyWlz3M9erzxDEDXrnd/UFa2y5elpYPfBasrKKj5aUt95jVr/n14RcWDx4CnmeUW5eRcfGp6+pAwwBZbTP6otPQnd1ZUPH5nRcXD3SF1WVraoPu6d5/6LsDixdv+tq6urLBPnxVn1F+vuvqtI1JStnglM/OA1Q3vVVa20zQI9wWoqnp5QlnZy+TmXr1fbu4Vi7Kzf7d2SrGq6uUVgEe3AdTUfHZEampBzAXO5r7BkbA2q7Cw0IuKWmyFnYiItDMffwxdu67bttdetNszqUIhm5+TM+7nXbpc912iY0kCeYsWMV1TVyIiIpK0lOiIiIhI0lKNjoiISBtRUODtdtqtrdKIjoiIiCQtJToiIiKStOKW6JhZfzP7j5l9ZmZzzez8SPsWZvaqmX0R+We3qPeMM7MvzWy+mY2OV2wiIiLSMcSzRqcGuMjdPzazXGCGmb0KnAK87u43mNllwGXApWY2FDgO2AEoAF4zs+3cfb0ttEVERABycmDFivWa8xIQirQ95RDHRMfdi4HiyONyM/sM6AscDuwT6fYo8CZwaaR9ortXAV+b2ZfArsD78YpRRETat+22W79t0SKmt34k0la1So2OmW0FjAA+BHpFkqD6ZKj+ULS+wPdRb1sYaWt4rTPNrMjMisrKyuIZtoiIiLRzcU90zCwHeBa4wN1XNte1kbb1tm129wfcvdDdC/Pz81sqTBEREUlCcU10zCydIMl5wt0nRZoXm1mfyOt9gNJI+0Kgf9Tb+wGheMYnIiIiyS2eq64M+AfwmbvfFvXSFODkyOOTgeej2o8zs0wz2xoYBJpnFRERkU0Xz1VXo4CTgE/MbFak7U/ADcDTZjYW+A44BsDd55rZ08D/CFZsnaMVVyIiIrI5dHq5iIgkm8ZqPqWD0s7IIiIikrSU6IiIiEjSUqIjIiIiSUuJjoiIiCQtJToiIiKStJToiIiISNJSoiMiIiJJS4mOiIiIJC0lOiIiIpK0lOiIiIhI0lKiIyIiIklLiY6IiIgkLSU6IiIikrSU6IiIiEjSUqIjIiIiSUuJjoiIiCSttA11MLMUYCegAKgE5rr74ngHJiIiIrK5mkx0zGwb4FLgAOALoAzoBGxnZhXA/cCj7l7XGoGKiIiIbKzmRnT+BtwL/NbdPfoFM+sJ/Bo4CXg0fuGJiIiIbLomEx13P76Z10qB2+MSkYiIiEgLabIY2cwGmdnzZvapmT1lZn1bMzARERGRzdXcqqsJwIvAUcDHwJ2tEpGIiIhIC2muRifX3R+MPL7ZzD5ujYBEREREWkpziU4nMxsBWOR5VvRzd1fiIyIiIm1ac4lOMXBb1POSqOcO7BevoERERERaQnOrrvZtzUBEREREWlqzOyObWXeC/XK2jzR9Bjzp7kvjHZiIiIjI5mpuefkQ4FNgF+Bzgt2RRwKfmtn2Tb1PREREpK1obkTnGuB8d386utHMjgKuJVh2LiIiItJmNbePzo4NkxwAd38W+En8QhIRERFpGc0lOqs38TURERGRNqG5qaueZvaHRtoNyI9TPCIiIiItprlE50Egt4nXHopDLCIiIiItqrl9dK5uzUBEREREWlpzy8vPMLNBkcdmZhPMbIWZzYkcBSEiIiLSpjVXjHw+8E3k8fHATsBA4A/A+PiGJSIiIrL5mkt0atw9HHl8CPCYu//g7q8B2fEPTURERGTzNJfo1JlZHzPrBOwPvBb1WlZ8wxIRERHZfM2turoCKAJSgSnuPhfAzPYGvmqF2EREREQ2S3Orrl40sy2BXHdfFvVSEXBs3CMTERER2UzNrbraw91rGiQ5uPtqd19lZl3MTEdBiIiISJvV3NTVUWZ2E/AyMAMoAzoB2wL7AlsCF8U9QhEREZFN1NzU1YVm1g04GjgG6ANUAp8B97v7O60TooiIiMimaW5Eh8i01YORHxEREZF2pbnl5ZslspNyqZl9GtU23Mw+MLNZZlZkZrtGvTbOzL40s/lmNjpecYmIiEjHEbdEB3gEOKhB203A1e4+nGD5+k0AZjYUOA7YIfKee8wsNY6xiYiISAcQt0TH3d8CljZsBrpEHncFQpHHhwMT3b3K3b8GvgR2RURERGQzNFujA2BmnQlWVw1w9/qDPge7+4ubcL8LgGlmdgtBkrV7pL0v8EFUv4WRNhEREZFNFsuIzsNAFfCzyPOFwN828X6/Ay509/7AhcA/Iu3WSF9v7AJmdmakvqeorKxsE8MQERGRjiCWRGcbd78JCAO4eyWNJyaxOBmYFHn8L36cnloI9I/q148fp7XW4e4PuHuhuxfm5+dvYhgiIiLSEcSS6FSbWRaRERYz24ZghGdThIC9I4/3A76IPJ4CHGdmmWa2NTAImL6J9xAREREBYqjRAa4k2B25v5k9AYwCTtnQm8zsKWAfoIeZLYxc5wzgDjNLA9YAZwK4+1wzexr4H1ADnOPutT8EPHkAACAASURBVBv9aURERESimHujpTDBi2ZGMI1UAfyUYMrqA3df0jrhNa+wsNCLiooSHYaIiLQtm1peIUloQzsju5lNdvddgH+3UkwiIiIiLSKWGp0PzGxk3CMRERERaWGx1OjsC/zWzL4FVhMMCbq7D4trZCIiIiKbKZZE5+C4RyEiIiISB7EkOk1XK4uIiIi0YbEkOv8mSHYM6ARsDcwnOIBTRESSzOSZi7h52nxCyyspyMviktGDGTNCp/JI+7TBRMfdd4x+bmY7A7+NW0QiIpIwk2cuYtykT6gMB1uZLVpeybhJnwAo2ZF2aaNPL3f3jwGtwhIRSUI3T5u/NsmpVxmu5eZp8xMUkcjmieX08j9EPU0BdgF0mqaISBIKLa/cqHaRti6WEZ3cqJ9M4EXg8HgGJSIiiVGQl7VR7SJtXSw1OlfXPzazFCDH3dfENSoREUmIS0YPXqdGByArPZVLRg9OYFQim26DIzpm9qSZdTGzbIJDN+eb2SXxD01ERFrbmBF9uf7IHembl4UBffOyuP7IHVWILO1Ws4d6ApjZLHcfbmYnENTnXArMaAs7I+tQTxERaYQO9ZS1YqnRSTezdGAM8Ly7h9EmgiIiItIOxJLo3A98A2QDb5nZlsDKeAYlIiIi0hJiKUYeD4yPavrWzPaNX0giIiIiLSOWYuTzI8XIZmb/MLOPgf1aITYRERGRzRLL1NVp7r4SOBDIB04FbohrVCIiIiItIJZDPeur138BPOzus81MFe0iIu2IDuqUjiqWEZ0ZZvYKQaIzzcxygbr4hiUiIi2l/qDORcsrcX48qHPyzEWJDi1QXQ2PPQYXX5zoSCQJxZLojAUuA0a6ewWQQTB9JSIi7UCbPahz+XK48UbYems4+WR45RWoqEhsTJJ0Ykl0HBgK/D7yPBvoFLeIRESkRbW5gzq/+QYuuAD694fLLoOhQ2HqVJg9Gzp3TkxMkrRiSXTuAX4GHB95Xg7cHbeIRESkRbWZgzqnT4djj4VttoG774YjjoCZM+HVV+Ggg0DlnxIHsSQ6u7n7OcAaAHdfRjB9JSIi7cAloweTlZ66TlurHdRZVwdTpsBee8Fuu8G0aUEtztdfB3U5w4fHPwbp0GJZdRU2s1Qixz6YWT4qRhYRaTfqV1e16qqrykp49FH4+9/h889hwIDg8dixkJsbv/uKNBBLojMeeA7oaWbXAkcDf45rVCIi0qLGjOjbOsvJS0vhnnuCqaklS6CwECZOhKOOgrRY/soRaVmxHAHxhJnNAPYn2FNnjLt/FvfIRESk/Zg/H267LRjFqaqCQw8Npqj23FO1N5JQzSY6ZpYCzHH3nwDzWickERFpF9zhrbfg1lvhhRcgMzNYJn7hhbD99omOTgTYQKLj7nVmNtvMBrj7d60VlIiItGE1NfDMM0GCU1QEPXrAlVfC2WdDz56Jjk5kHbFMmPYB5prZdGB1faO7Hxa3qEREpO0pL4eHHoI77oBvv4XttoP77oPf/AayWnmpukiMYkl0ro57FCIi0nYtXAjjx8MDD8CKFcFS8fHj4ZBDICWWXUpEEmdDNTpjgG2BT9x9WuuEJCIibcKsWcH01MSJwX44Rx8NF10Eu+6a6MhEYtZkomNm9wA7AO8B15jZru5+TatFJiIirc892NTvllvg9dchOxvOOQfOPz84k0qknWluRGcvYCd3rzWzzsDbgBIdEZFkVFUFTz4ZjODMnQsFBcGBm2eeCXl5iY5OZJM1l+hUu3stgLtXmGkjBBGRpLN0aVBQfOedUFICw4YFe+Ecdxxk6LQfaf+aS3S2N7M5kccGbBN5boC7+7C4RyciIvGxYAHcfjtMmAAVFTB6dHD21AEHaIM/SSrNJTpDWi0KERFpHe+/H0xPTZoUHMlwwgnwhz/AjjsmOjKRuGgy0XH3b1szEBERiZPaWnj++SDBee+9oObmssvg3HODWhyRJKYT1kREktXq1fDII8Gp4QsWBKumxo+HU0+FnJxERyfSKpToiIgkm5ISuOsuuPfeoNh4t93ghhvgiCMgNTXR0Ym0KiU6IiJt2OSZi7h52nxCyyspyMviktGDGTOib+PtmSuCE8QffxzCYTj88OAE8d13V4GxdFjNbRj4CeBNva5VVyIi8TV55iLGTfqEynAtAIuWVzJu0icUfbuUZ2csCtrdGTD7A7Z46FJYUBScOTV2bHCC+KBBCf4EIonX3IjOIZF/nhP55+ORf54AVMQtIhERAeDmafPXJjn1KsO1PPXh91hNmMPnvc2Z059jh9KvKOucx4MHnMIZT90cnCYuIkAMq67MbJS7j4p66TIzexf4a7yDExHpyELLK9dry61azXGzpnHqjCkUlC/hi+79ufSg85i8w75Up2VwhpIckXXEUqOTbWZ7uPs7AGa2O5Ad37BERKQgL4tFkWSnYGUppxZN4bjZ08itruS9AcO4fPQ5vDlwF9yCE8T75mUlMlyRNimWRGcsMMHMukaeLwdO29CbzGwCwfRXqbv/JKr9POBcoAb4t7v/MdI+LnKvWuD3Oi1dZPM0VcQq7cclowfz+N2T+M17z/LLeW8DMHXoXnx36lnctbzrOtNaWempXDJ6cKJCFWmzNpjouPsMYCcz6wKYu6+I8dqPAHcBj9U3mNm+wOHAMHevMrOekfahwHEEp6UXAK+Z2Xb1Z22JyMZpqogVULLTHtTVwdSpjLnlFsa8+SarMzvzcOHhTN3vGH5z7F6cM6IvfZXIisTE3JtcWBV0CEZyriQ4zRzgv8BfY0l4zGwr4MX6ER0zexp4wN1fa9BvHIC7Xx95Pg24yt3fb+76hYWFXlRUtKEwRDqcUTe8sXbKI1rfvCzevWy/BEQkMVmzBv75z2AH43nzoF8/uOACOP106Np1w++XelpLL2ulxNBnAlAO/CrysxJ4eBPvtx2wp5l9aGb/NbORkfa+wPdR/RZG2tZjZmeaWZGZFZWVlW1iGCLJrbEi1ubaJcGWLIFrroEtt4QzzgiWiD/xBHz1FVx0kZIckc0QS43ONu5+VNTzq81s1mbcrxvwU2Ak8LSZDaTx7LvRoSZ3fwB4AIIRnU2MQySpRRexNmyXNuSLL4LjGR55BCor4Re/CDb422cfbfAn0kJiSXQqG6y6GgVs6v8WLgQmeTBfNt3M6oAekfb+Uf36AaFNvIdIh3fJ6MHr1OhA+y5WTarCand4991geur55yE9HU46KThBfOjQuN02qb5DkY0QS6LzO+DRSK2OAUuBkzfxfpOB/YA3zWw7IANYAkwBnjSz2wiKkQcB0zfxHiIdXv1fYMnwF1vSFFbX1MBzzwUJzocfwhZbwOWXwznnQO/ecb110nyHIpsgllVXs/hx1RXuvjKWC5vZU8A+QA8zW0hQ0DyBYKn6p0A1cHJkdGdupFD5fwTLzs/RiiuRzTNmRN+k+Eusqd2Bb542v318vlWrYMIEuP12+Ppr2HZbuPtuOPlkyG6dLcna/Xcoshk2mOg0XHVlZjGtunL345t46cQm+l8LXLuheESkY2m3hdWhENx5J9x3HyxfHhyseeutcNhhrX6CeLv9DkVaQGuvuhIR2ShNFVC32cLqTz6BU06BrbaCm26C/feH994L6nKOOKLVkxxoh9+hSAuKJdHZxt2vdPevIj9XAwPjHZiICASF1Vnp6yYHba6w2h1efRVGj4Zhw+Bf/4KzzgpWVT3zDPzsZwkNr118hyJx0tqrrkRENkqbLqyuroaJE4MpqTlzgqLia68Nkpwttkh0dGu16e9QJM5i2Rl5OPAoEL3q6hR3nx3/8JqnnZFFJCGWL4f774fx44NanB12CDb2+/WvITMz0dGJdkaWKHFbdSUiknS++SZYPfXQQ7B6NRxwQLCi6sADtcGfSBvVZKJjZn9ooh0Ad78tTjGJiLQt06cH01PPPAMpKXD88cEGf8OHJzoyEdmA5kZ0clstChGRtqauDl54IUhw3n47OG/q4ovhvPOCwzZFpF1oMtGJrK4SEUkKMR+BUFkJjz4anEH1+ecwYEDweOxYyNX//4m0N81NXXUCjgWWAS8AlxBsGrgAuMbdl7RKhCIimymmIxBKS+Gee4Jdi5csgcLCYEXVUUdBWiwLVEWkLWru397HgDCQDVwEfArcBewBPAIcEu/gRERaQrNHIGSVw223wWOPQVUVHHpoMEW1554qMBZJAs0lOkPd/SdmlgYsdPe9I+0vm1nCl5aLiMRqvaMO3Nnt+08585lJMO6jYEn4ySfDhRfC9tsnJkgRiYvmEp1qAHevMbNQg9d04KaItBsFeVksWl5Jal0tv5j3Dqd/NJmdSr5gWXZXuPJKOPts6Nkz0WGKSBw0l+j0M7PxBBsv1T8m8lzbaYpIuzFuj758+tfbOPHDyfRbWcqCLfpy5cHnsctffs9hP9s20eGJSBw1l+hcEvW44fbD2o5YRNq+hQth/HgOuf9+Dlm5kplbD+PqA37LZ7vsxcUHD+EwHYEgkvSaW17+aGsGIiLSYmbNCva/mTgx2A/nmGPgoosYMXIkDyY6NhFpVVozKSLJwR2mTYNbboHXX4fsbDj3XDj/fNhqq0RHJyIJokRHRNq3qip48slgBGfuXCgogBtvhDPPhLy8REcnIgmmREdE2qelS+G+++DOO6GkBIYNC3Y0Pu44yMhIdHQi0kZsMNExs0eB8919eeR5N+BWdz8t3sGJiKxnwYLgBPEJE6CiAkaPhscfh/33b3aDv5iPgBCRpBLLiM6w+iQHwN2XmdmIOMYkIrK+998PpqcmTQqOZDjhhOAE8R133OBbYzoCQkSSUkosfSKjOACY2RZoyktEWkNtbZDYjBoFu+8eFBlfdhl88w08/HBMSQ40fwSEiCS3WBKWW4H3zOyZyPNjgGvjF5KIdHirV8MjjwSnhi9YAFtvDePHw6mnQk7ORl9uvSMgNtAuIsljg4mOuz9mZkXAfgS7Ih/p7v+Le2Qi0vGUlMBdd8G99wbFxrvtBjfcAEccAampm3zZ+iMgGmsXkeTW5NSVmXWJ/HMLoAR4EngCKIm0iYi0jLlzYexY2HJLuO462HtveOedoC7n6KM3K8kBuGT0YLLS171GVnoql4wevFnXFZG2r7kRnSeBQ4AZgEe1W+T5wDjGJSLJzh3+859gg7+pUyErK0h2LrwQBg1q0VvVFxxr1ZVIx2PuvuFebVRhYaEXFenYLZF2JRyGp58OEpxZs4JTw889F373O+jRI9HRSXJoep8B6XCaHNExs52be6O7f9zy4YhI0lqxAh58EO64Izhsc8gQeOihYJl4p06Jjk5EklRzU1e3NvOaExQni4g077vvguTmwQehvBz23TfY0fjggyEllh0uREQ2XXOnl+/bmoGISJKZMSPY4O/pp4Pnxx4LF10EOzc7WCwi0qJiOQIiHfgdsFek6U3gfncPxzEuEWmP6urgpZeCBOfNNyE3Fy64AH7/exgwINHRiUgHFMuGgfcC6cA9kecnRdpOj1dQItLOrFkD//xnkODMmwf9+gXFxqefDl27Jjo6EenAYkl0Rrr7TlHP3zCz2fEKSETajg0ehLlkSbC53113QWkpjBgBTzwBxxwD6emJC1xEJCKWRKfWzLZx9wUAZjYQqN3Ae0SknWv2IMyciuB4hkcegcpK+MUv4OKLYZ99mj1BXESktcWS6FwC/MfMviLYm2BL4LS4RiUiCbfeQZju7PD1HLqd+Df47L1gxOakk4ITxIcOTVygIiLNiCXReQcYBAwmSHTmxTUikSS0wSmgNqj+wMvUulpGf/4+Z0x/jhHF81nWKRcuvxzOOQd6905wlMmlPf6eiLR1sSQ677v7zsCc+gYz+xjQGlGRGDQ7BdSG/xLbJgv2eHsKY4uep/+KxXzdrQ9//vnveH+PX/L6Fb9MdHhJp73+noi0dc3tjNwb6AtkmdkIftxSuwvQuRViE0kK600BAZXhWm6eNr9t/gUWCsGdd/LS3feSUb6Cj/oO5W/7jeXVbXcjMzOD6w/dMdERJqV293si0k40N6IzGjgF6EewS3J9olMO/Cm+YYkkj/opoFjbE+aTT4Ll4U8+CbW1ZBx5JP895CT+VJyjqZRW0G5+T0TameZ2Rn4UeNTMjnL3Z1sxJpGkUpCXxaJG/rIqyMtKQDQNuMNrrwV73rzyCnTuDGedFWzyN3AgewPvJjrGDqJN/56ItGOxHDTTz8y6WOAhM/vYzA6Me2QiSeKS0YPJSk9dpy0rPZVLRg9OUERAdTU89hgMHw4HHsiaj2dx/4GnMfz0fzCqYAyTV2S22K0mz1zEqBveYOvL/s2oG95g8sxFLXbtZNImf09EkkAsxcinufsdZjYa6AmcCjwMvBLXyESSRP1UT5tYTbNsGTzwAIwfH9Ti7LADH191G6es2YaVHvwlu7wFi2BVYBu7NvV7IpJEzN2b72A2x92HmdkdwJvu/pyZzXT3Ea0TYtMKCwu9qKgo0WGItH3ffAO33w4PPQSrV8MBBwQb/B14IKNu/E+jUyZ987J497L9Nuu2o254I27XFmmGdq2UtWIZ0ZlhZq8AWwPjzCwXqItvWCLSIqZPDwqMn3kGUlLg+OODDf6GD1/bJZ5FsCqwFZFEiyXRGQsMB75y9woz604wfSUibVFdHbzwQpDgvP12cKjmxRfDeecFh202EM8iWBXYikiixVKM/DTQB1gJ4O4/uPuc5t8CZjbBzErN7NNGXrvYzNzMekS1jTOzL81sfqQeSEQ2RmUl3HcfDBkCY8bAt98G51F9/z3ceGOjSQ7Etwg2GQtsVVwt0r7EMqJzH8EIzngz+xfwiLvHcgzEI8BdwGPRjWbWH/g58F1U21DgOGAHoAB4zcy2c3cdHiqyIaWlcPfdcM89wWnihYUwcSIcdRSkbfhf8XgWwSZbga2Kq0Xanw0WI6/taNYVOB64HPgeeBD4p7uHm3nPVsCL7v6TqLZngGuA54FCd19iZuMA3P36SJ9pwFXu/n5zMakYWTq0efPgttuCZeJVVXDoocEU1Z576gTxOFFxdbuhfwFkrVimrojU5ZwCnA7MBO4gOOvq1Y25mZkdBixy99kNXupLkDzVWxhpa+waZ5pZkZkVlZWVbcztRdo/d/jvf4OkZsgQePxxOPlk+OwzmDIF9tpLSU4cqbhapP3Z4Li2mU0CtgceBw519+LIS/9nZjEPp5hZZ4LRoMY2G2zsv8yNDjW5+wPAAxCM6MR6f5F2raYmWDl1yy0wYwb06AFXXQVnnw35+YmOrsNQcbVI+xPLiM5d7j7U3a+PSnIAcPfCjbjXNgRL1Geb2TcEZ2h9HDk8dCHQP6pvPyC0EdcWSU7l5UFB8TbbBEvDy8uDguPvvoMrr1SS08qSsbhaJNk1d3r5SOB7d38j8vw3wFHAtwT1M0s35kbu/gnBzsr11/+GH2t0pgBPmtltBMXIg4DpG/lZRJLHwoXB7sX33w8rVwZTUnfeCYccEuyHs4kmz1yUNIXBiZBsxdUiHUFzU1f3AwcAmNlewA3AeQR76jwAHN3chc3sKWAfoIeZLQSudPd/NNbX3eea2dPA/4Aa4BytuJIOadasYP+biROD/XCOOQYuughGjtzsS2vFUMsYM6Kvvi+RdqTJVVdmNtvdd4o8vhsoc/erIs9nufvwRt/YirTqSpKCO0ybFtTfvP46ZGfDGWfA+efDVlu12G20Ykg6EFXky1rNjeikmlmau9cA+wNnxvg+EYlFVRU8+WQwgjN3LhQUBBv7nXkm5OW1+O20YkhEOqLmEpangP+a2RKgEngbwMy2BVa0QmwiyWnp0qCg+M47oaQEhg0L9sI59ljIyIjbbbViSEQ6oiYTHXe/1sxeJzj+4RX/cY4rhaBWR6RDaLEC3gULghPEJ0yAigoYPTrYB2f//Vtl75tLRg9ep0YHtGJIRJJfs1NQ7v5BI22fxy8ckbalRQp4338/mJ6aNCk4kuGEE4ITxHfcMV5hN0orhkSkI4r5CIi2SMXIEm+bXMBbWwvPPx8kOO+9B926wVlnwbnnBrU4IhJPKkaWtVRULNKMjS7gXb0aHnkk2ORvwQLYeutgP5xTT4WcnPgFGiPtoyMiHY0SHZFmxFzAW1ICd90F994bFBvvthvccAMccQSkpq73/kTQPjoi0hFt+harIh3ABrf8nzsXxo6FLbeE666DvfeGd94J6nKOPrrNJDkQ1OZEFyIDVIZruXna/ARFJCISfxrREWlGowW8B27HmGXz4RdnwNSpkJUFp58OF1wAgwYlOOKmaR8dEemIlOiIbMDaLf/DYXj6aRh7UXBUQ8+ecM01QZFxjx6JDnODtI+OiHREmroS2ZAVK4LjGQYOhBNPDHY0fugh+PZb+POf20WSAzp5W0Q6Jo3oiDTlu+/gjjvgwQehvBz23TfY0fjggzfrBPFE0T46ItIRKdERaWjGjGD/m6efDp4fe2xwgvjOOyc2rhagk7dFpKNRoiMCUFcHL70UJDhvvgm5uUFx8e9/DwMGJDo6ERHZREp0pGNbswb++c8gwZk3D/r1C+pxTj8dunZNdHQiIrKZlOhIx7RkCdxzD9x9N5SWwogR8MQTcMwxkJ6e6OhERKSFKNGRjuXzz4PjGR59FCor4Re/gIsvhn32aZUTxEVEpHUp0ZHk5w7vvhtMSU2ZEozYnHRScIL40KGJjk5EROJIiY4kXNwOmqypgeeeCxKc6dNhiy3g8svhnHOgd+/Nv34HpsNBRaS9UKIjCRWXgyZXrYIJE4Ipqm++gW23DWpxTj4ZsrNbKPKOS4eDikh70v52PZOk0qIHTYZCMG4c9O8P558PffvCpEnBaqqzz1aS00J0OKiItCca0ZGEapGDJj/5JFge/uSTUFsLRx4ZbPD305+2UJQSTYeDikh7ohEdSaimDpTc4EGT7vDqqzB6NAwbBv/6V3C45hdfBI+V5MTNJv+ZiYgkgBIdSaiNPmiyujpYGj58OBx4IMyZA9ddB99/D+PHBwdvJqnJMxcx6oY32PqyfzPqhjeYPHNRQuLQ4aAi0p5o6koSKuaDJpctgwceCJKZUAh22AEefhiOPx4yMxMQeetqSwXAOhxURNoTc/dEx7DJCgsLvaioKNFhSDx98w3cfjs89BCsXg0HHBBs8HfggR1qg79RN7zBokZqYPrmZfHuZfslICKRNq3j/MdBNkgjOtI2TZ8eFBg/8wykpAQjNxddBDvtlOjIEkIFwCIim0aJjrQddXXwwgtBgvP228GhmhdfDOedFxy22YEV5GU1OqKjAmARkeapGFkSr7IS7rsPhgyBMWPgu++Czf6+/x5uvLHDJzmgAmARkU2lER1JnNLSYMfie+4JThMvLISJE+GooyBNv5rRVAAsIrJp9LeJtL558+C22+Cxx6CqCg47LKi/2XPPDlVgvLHGjOirxEZEZCMp0ZGYbPYhju7w1lvBAZsvvgidOsEpp8CFF8JgTb+IiEh8KNGRDdqsPVxqaoKVU7fcAjNmQI8ecNVVwdlT+flxjlxERDo6FSPLBm3SIY7l5UFB8TbbBEvDy8uDguPvvoMrr1SSIyIirUIjOm3YZk8XtZCN2sNl4cJg9+L774eVK2GvveCuu+CXvwz2w5FN1lZ+H0RE2hMlOm1UW9ryP6Y9XGbNCva/mTgx2A/nmGOCAuORI1sx0uTVln4fRETaE/0vdhu1SdNFcdLkHi4HbgcvvxwcyzBiBEyeDOeeCwsWBAmPkpwW05Z+H0RE2hON6LRRbWnL/4Z7uGyZk8pt4bns/JuLYe5cKCgINvY780zIy2v1+DqCtvT7ICLSnijRaaPa2pb/Y0b0ZcyATkFB8V13QUkJDBsW7IVz7LGQkZGQuDqKtvb7ICLSXmjqqo1qU1v+L1gQnDc1YAD8+c/BwZqvvhrU5Zx0kpKcVtCmfh9ERNoRjei0UW1iy//33w/2v3nuueBIhhNOgD/8AXbcsfViEKCN/D6IiLRD5u6JjmGTFRYWelFRUaLDSC61tfD880GC8/770K0bnHVWUGRcUJDo6EREYqGzZGQtjehIYPVqeOSRYJO/BQtg662D/XBOPRVychIdnYiIyCZRotPRlZQExcX33gtLl8Juu8ENN8ARR0Bq6obfLyIi0obFrRjZzCaYWamZfRrVdrOZzTOzOWb2nJnlRb02zsy+NLP5ZjY6XnFJxNy5MHYsbLklXHcd7L03vPNOMF119NFKckREJCnEc0TnEeAu4LGotleBce5eY2Y3AuOAS81sKHAcsANQALxmZtu5ey0Ssw0eEeAOb7wR7GA8dSpkZcHpp8MFF8CgQYkLvAPQ8Q0iIokRt0TH3d8ys60atL0S9fQD4OjI48OBie5eBXxtZl8CuwLvxyu+ZNPsEQE/6QlPPx0UGM+aBT17wjXXBEXGPXokMuwOQcc3iIgkTiL30TkNmBp53Bf4Puq1hZE2iVFjRwSkrVrJwsuvgYED4cQToaoKHnoIvv022A9HSU6r0PENIiKJk5BiZDO7HKgBnqhvaqRbo+vezexM4EyAAQMGxCW+9ij6KICClaWcWjSF42ZPI7e6EvbdNzhN/KCDdIJ4Auj4BhGRxGn1RMfMTgYOAfb3HzfxWQj0j+rWDwg19n53fwB4AIJ9dOIYartSkJdFt3mfcMb05/jlvLcBeHHInjy/33E8fOdZCY6uY9PxDSIiidOqiY6ZHQRcCuzt7hVRL00BnjSz2wiKkQcB01sztnarrg5eeonJz15HftH7lGdkMaHwcB4pPJRl3ftw/ZHaxTjRLhk9eJ0aHdDxDSIirSVuiY6ZPQXsA/Qws4XAlf/f3t3Hal3WcRx/f0VMfEAiyAdAmDMx0xSfhpLO1OXjiAhdbCk2W8ungaktcg2rqbU2RqY0nFJg6ZkikrKmMrU0WzEfUFARLZ8Q9Mh8RFFDvv3x+4FHFD0onPvc1/1+bWfnd//uc9j13Q3ss+u6ft+L6imrzwHzIgLgX5n5w8x8JCKuBx6lWtI6kkiMiAAACQZJREFUyyeuPsHbb8M118DkybB4Mf0HDWLRuT/jR70P5Im3e7BLn15c6pM93YLHN0hS43gERLNZsQKmToUrroD2dhg2DM4/H046CXr2bPToJKk78AgIrWNn5GaxZEl1PMOMGbBqFRx/fBVwjjgCwn/TkiR9FINOd5YJ995b9b+5+eZqxubUU+Hcc2GvvRo9OkmSuj2DTne0ejXcdFMVcObPh7594cIL4ayzYKedqi67v7rT/R6SJH0Cg053snIlTJ9eLVE9/TTsvnu1F2fcONh2W8Auu5IkbQy7x3UHy5bBxIkwaBCMHw8DBsDs2bB4MZx55rqQA3bZlSRpYzij00gLF1YHbF57Lbz3HoweDeedB8OHb/BX7LIrSVLnGXS6WibMm1cFnNtvh222qQ7XnDChOpPqE9hlV5KkznPpqqu8+271aPi++8Ixx8DDD8Mll8Bzz8Fll3Uq5EDVZbdXzx4fuGeXXUmSPpozOpvbK6/AlVdWYWbZMl7bfSiXj7mAmYMPpV/05oJnVjGqb+f/OLvsSpLUeXZG3lyeegqmTIGrr4Y334Sjj+afo07j9OV9WbV6zbof69WzB5eO3segIkmbjl1UtY5LV5va/Plw8snVo+FTp1YbjBcsgHnzuOCNnT8QcsAnpiRJ2pxcutoU1qyBW26pNhjfcw/ssEN1PMM558DAget+zCemJEnqWgadz+Ktt2DmzKrB35IlMHhwdX366bD99h/6cZ+YkiSpa7l09Wm0t8OkSbDrrnDGGdC7N7S1wZNPVo+Jf0TIAZ+YkiSpqzmjszEWL4bJk6tZnHfegZEjqwZ/hx3WqRPEfWJKkqSuZdD5JJlw993VAZtz58LWW8Npp1UniA/d+JmYUcMGGGwkSeoiBp0NWb0aZs2qAs7990O/fnDRRdXZU/37N3p0kiSpEww663v99ar3zZQp8OyzsMceMG0anHIK9HLTsCRJzcSgs9bSpVX34mnTqrBz+OFw+eVwwgmwhXu2SzLnwefdJyVJLcKgs2BB1f+mra3ajzNmTLXB+KCDGj0ybQZzHnyeibMXsup/7wHw/KurmDh7IYBhR5IK5FTFpEkwZw6cfXb1eHhbmyGnYL+57fF1IWctu1NLUrmc0bnssqqTcZ8+jR6JuoDdqSWptTijM3iwIaeFbKgLtd2pJalMBh21FLtTS1JrcelKLcXu1JLUWgw6ajl2p5ak1uHSlSRJKpZBR5IkFcugI0mSimXQkSRJxTLoSJKkYhl0JElSsQw6kiSpWAYdSZJULIOOJEkqlkFHkiQVy6AjSZKKZdCRJEnFisxs9Bg+tYh4CXhmA2/3A1Z04XAawRrL0Qp1WmM5unudKzLz2EYPQt1DUwedjxMR92XmgY0ex+ZkjeVohTqtsRytUqfK4NKVJEkqlkFHkiQVq+Sgc2WjB9AFrLEcrVCnNZajVepUAYrdoyNJklTyjI4kSWpxTR90ImJ6RLRHxKIO9/pGxLyIeKL+/vlGjvGziohBEXFXRDwWEY9ExPj6fml1bh0R8yPiobrOn9f3i6oTICJ6RMSDETG3fl1ijU9HxMKIWBAR99X3iqozIvpExKyIWFz/+zykpBojYmj9+a39ej0iJpRUo8rX9EEH+COwfr+EnwB3ZOaXgDvq181sNXBeZn4ZGA6cFRF7UV6d7wBHZua+wH7AsRExnPLqBBgPPNbhdYk1Anw9M/fr8ChyaXX+Frg1M/cE9qX6TIupMTMfrz+//YADgLeAmyioRpWv6YNOZt4NvLze7W8CM+rrGcCoLh3UJpaZyzPzgfr6Dar/TAdQXp2ZmSvrlz3rr6SwOiNiIHACcFWH20XV+DGKqTMiegOHA1cDZOa7mfkqBdW4nqOA/2TmM5RbowrU9EFnA3bMzOVQhQTgiw0ezyYTEUOAYcC/KbDOeklnAdAOzMvMEuucAvwYWNPhXmk1QhVSb4+I+yPiB/W9kurcDXgJ+EO9DHlVRGxLWTV29B3guvq61BpVoFKDTpEiYjvgRmBCZr7e6PFsDpn5Xj1NPhA4OCL2bvSYNqWIOBFoz8z7Gz2WLjAiM/cHjqNabj280QPaxLYE9gd+n5nDgDcpdAknIrYCRgI3NHos0sYqNei8GBE7A9Tf2xs8ns8sInpShZw/Z+bs+nZxda5VLwH8jWr/VUl1jgBGRsTTQBtwZET8ibJqBCAzl9Xf26n2dRxMWXUuBZbWs44As6iCT0k1rnUc8EBmvli/LrFGFarUoHMzMK6+Hgf8pYFj+cwiIqj2ATyWmZM7vFVanf0jok993Qs4GlhMQXVm5sTMHJiZQ6iWAu7MzO9SUI0AEbFtRGy/9hr4BrCIgurMzBeA5yJiaH3rKOBRCqqxg7G8v2wFZdaoQjV9w8CIuA44guo03ReBScAc4HpgV+BZ4KTMXH/DctOIiK8B9wALeX9fx0+p9umUVOdXqTY29qAK4ddn5i8i4gsUVOdaEXEEcH5mnlhajRGxG9UsDlRLPNdm5sUF1rkf1abyrYD/At+j/rtLOTVuAzwH7JaZr9X3ivocVbamDzqSJEkbUurSlSRJkkFHkiSVy6AjSZKKZdCRJEnFMuhIkqRiGXSkJhER34qIjIg9Gz0WSWoWBh2peYwF/kHVaFCS1AkGHakJ1OecjQBOpw46EbFFREyNiEciYm5E/DUixtTvHRARf68P1Lxtbbt+SWo1Bh2pOYwCbs3MJcDLEbE/MBoYAuwDfB84BNadi/Y7YExmHgBMBy5uxKAlqdG2bPQAJHXKWGBKfd1Wv+4J3JCZa4AXIuKu+v2hwN7AvOqYNHoAy7t2uJLUPRh0pG6uPlfoSGDviEiq4JK8f5bUh34FeCQzD+miIUpSt+XSldT9jQFmZubgzBySmYOAp4AVwLfrvTo7Uh1uC/A40D8i1i1lRcRXGjFwSWo0g47U/Y3lw7M3NwK7AEuBRcA0qtPsX8vMd6nC0a8j4iFgAXBo1w1XkroPTy+XmlhEbJeZK+vlrfnAiMx8odHjkqTuwj06UnObGxF9gK2AXxpyJOmDnNGRJEnFco+OJEkqlkFHkiQVy6AjSZKKZdCRJEnFMuhIkqRiGXQkSVKx/g/AlwC2AdD6zgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plot setup\n",
    "plt.figure(figsize=(7,6)) \n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "\n",
    "# load data from a csv file\n",
    "data = np.genfromtxt('files/blood_pressure.csv', delimiter=',')[1:,:]\n",
    "age = data[:,0]\n",
    "sbp = data[:,1]\n",
    "\n",
    "# the least squares line\n",
    "meanx = np.mean(age)\n",
    "meany = np.mean(sbp)\n",
    "varx = np.var(age)\n",
    "# cov = np.cov(age, sbp)[0,1]\n",
    "cov = np.inner(age-meanx, sbp-meany)/age.shape[0]\n",
    "\n",
    "xvals = np.arange(np.min(age)-5, np.max(age)+5, 0.1)\n",
    "yvals = cov*(xvals - meanx)/varx + meany\n",
    "\n",
    "slope = cov / varx\n",
    "intercept = meany - cov*meanx / varx\n",
    "\n",
    "plt.xlabel(\"Age\")\n",
    "plt.ylabel(\"Systolic Blood Pressure (SDP)\")\n",
    "plt.plot(xvals, yvals, c=\"r\", zorder=1)\n",
    "\n",
    "plt.figtext(0.9,0.8, \" slope = {:.4f}\".format(slope) +\n",
    "            \"\\n intercept = {:.4f}\".format(intercept),\n",
    "            ha=\"left\", va=\"top\",\n",
    "            backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
    "plt.title(\"Simple linear regression\")\n",
    "plt.scatter(age, sbp, zorder=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAFkCAYAAACuIwUtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAbT0lEQVR4nO3df3TV9X3H8dcbsoCBoAYSRCCmVJjAXEEyga6s/qKyyYrWH8dWC2117NS1a9fuWFxdPTtDq63H0R+6jm21sNI627nJyRlYpXNWR2ABmUdhTJdmgZJKSEBDUhIS3vvj3sQbTELC/d77vZ97n49zOLn3+733+/348d7v634/n8/38zV3FwAAIRoVdwEAADhbhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYECMze83Mrhhk3RVmdjCi/TxvZndGsS0glxTFXQAgFGbWIGmypB5JxyVtlfQZdz9+ttt097nRlA4oTJyJASPz++4+XtI8SfMl3RNzeYCCRogBZ8HdfynpGSXCTGY2xsweNrNGM3vTzL5jZuck100ysxozO2ZmrWb2MzMblVzXYGbXJB+fY2bfM7OjZrZX0m+l7tPM3MwuTnn+PTNbm3x8fnIfzcn315jZtIHKbmYXm9m/m9lbZnbEzP4xA1UEZAUhBpyFZED8rqQ3kosekjRLiVC7WNJUSV9JrvuipIOSypVojvwzSQPN93afpPcm/10radUIijRK0uOSLpJUKelXkr49yGv/UtJPJJ0vaZqkb41gP0BOIcSAkfkXM2uTdEDSYUn3mZlJ+gNJf+Lure7eJukBSbcm33NS0hRJF7n7SXf/mQ88aektku5PbuOApG8Ot1Du3uLu/+TuHcn93y/pg4O8/KQSYXehu59w9xeHux8g1xBiwMhc7+6lkq6QdImkSUqcYZVI2pVsMjymxKCP8uR7vq7EGdtPzKzezNYMsu0LlQjHXv833EKZWYmZ/Y2Z/Z+ZvS3pBUnnmdnoAV5+tySTtDM5OvJTw90PkGsIMeAsuPu/S/qepIclHVGi+W6uu5+X/HducgCI3L3N3b/o7jMk/b6kL5jZ1QNstknS9JTnlaet71AiLHtdkPL4i5J+XdJCd58g6XeSy22Asv/S3f/A3S+U9IeSHkvtawNCQogBZ2+dpKWSflPS30r6KzOrkCQzm2pm1yYfL08OpjBJbysxRL9ngO09Keme5CCNaZI+e9r6PZI+ZmajzWyZ+jcXlioRpMfMrEyJ/rUBmdnNKYM+jirRPzdQeYCcR4gBZ8ndmyVtlPTnkr6kRJNhbbI57zklzowkaWby+XFJ2yU95u7PD7DJv1CiCfHnSgy8+IfT1n9OiTO5Y5Juk/QvKevWSTpHibPCWiWaMwfzW5J2mNlxSZslfc7df37m/2Ig9xg3xQQAhIozMQBAsAgxAECwCDEAQLAiCTEz+5Pk9SavmtkPzWysmZWZ2bNm9nry7/lR7AsAgF5ph5iZTZX0x5Kq3f03JI1WYqaCNZK2uftMSduSzwEAiExUt2IpknSOmZ1U4mLMQ0rM7n1Fcv0GSc8rMQx5UMuWLfOtW4caGQwAKEDvumi/V9pnYu7+CyVmLWhUYsaBt9z9J5Imu3tT8jVNkioGLJnZajOrM7O6ffv2pVscAEABiaI58XxJKyS9R4m538aZ2e3Dfb+7r3f3anevLi8vP/MbAABIimJgxzWSfu7uze5+UtJTkt4v6U0zmyJJyb+HI9gXAAB9ogixRkmLkrNom6SrJe1TYjqb3vshrZL0dAT7AgCgT9oDO9x9h5n9WNJuSd2SXpa0XtJ4SU+a2R1KBN3N6e4LAIBUOTV3YnV1tdfV1cVdDABAbsnc6EQAAOJCiAEAgkWIAQCCRYgBAIIV1bRTAAJUW9+ijdsb1NjaocqyEq1cXKVFMybGXSxg2DgTAwpUbX2L1tbs1ZG2LpWPH6MjbV1aW7NXtfUtcRcNGDZCDChQG7c3qKS4SKVjizTKTKVji1RSXKSN2xtiLhkwfIQYUKAaWzs0bszofsvGjRmtxtaOmEoEjBwhBhSoyrIStXf29FvW3tmjyrKSmEoEjBwhBhSolYur1NHVrbYT3TrlrrYT3ero6tbKxVVxFw0YNkIMKFCLZkzUvcvnaFJpsZqPd2pSabHuXT6H0YkICkPsgQK2aMZEQgtB40wMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABAsQgwAECxCDAAQLEIMABCsSELMzM4zsx+b2X+b2T4zW2xmZWb2rJm9nvx7fhT7AgCgV1RnYt+QtNXdL5H0Pkn7JK2RtM3dZ0ralnwOAEBk0g4xM5sg6Xck/b0kuXuXux+TtELShuTLNki6Pt19AQCQKoozsRmSmiU9bmYvm9nfmdk4SZPdvUmSkn8rBnqzma02szozq2tubo6gOACAQhFFiBVJukzSX7v7fEntGkHTobuvd/dqd68uLy+PoDgAgEIRRYgdlHTQ3Xckn/9YiVB708ymSFLy7+EI9gUAQJ+0Q8zdfynpgJn9enLR1ZL2StosaVVy2SpJT6e7LwAAUhVFtJ3PStpkZsWS6iV9UomAfNLM7pDUKOnmiPYFAICkiELM3fdIqh5g1dVRbB8AgIEwYwcAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiEGAAgWIQYACBYhBgAIFiRhZiZjTazl82sJvm8zMyeNbPXk3/Pj2pfAABI0Z6JfU7SvpTnayRtc/eZkrYlnwMAEJlIQszMpkm6TtLfpSxeIWlD8vEGSddHsS8AAHpFdSa2TtLdkk6lLJvs7k2SlPxbMdAbzWy1mdWZWV1zc3NExQEAFIK0Q8zMlks67O67zub97r7e3avdvbq8vDzd4gAACkhRBNv4bUkfNrPfkzRW0gQz+76kN81sirs3mdkUSYcj2BcAAH3SPhNz93vcfZq7V0m6VdJP3f12SZslrUq+bJWkp9PdFwAAqTJ5ndiDkpaa2euSliafAwAQmSiaE/u4+/OSnk8+bpF0dZTbBwAgFTN2AACCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgkWIAQCCRYgBAIJFiAEAgpV2iJnZdDP7NzPbZ2avmdnnksvLzOxZM3s9+ff89IsLAMA7ojgT65b0RXefLWmRpD8yszmS1kja5u4zJW1LPgcAIDJph5i7N7n77uTjNkn7JE2VtELShuTLNki6Pt19AQCQKtI+MTOrkjRf0g5Jk929SUoEnaSKKPcFAEBkIWZm4yX9k6TPu/vbI3jfajOrM7O65ubmqIoDACgARVFsxMx+TYkA2+TuTyUXv2lmU9y9ycymSDo80Hvdfb2k9ZJUXV3tUZQHyEe19S3auL1Bja0dqiwr0crFVVo0Y2LcxQJiFcXoRJP095L2ufsjKas2S1qVfLxK0tPp7gsoVLX1LVpbs1dH2rpUPn6MjrR1aW3NXtXWt8RdNCBWUTQn/rakj0u6ysz2JP/9nqQHJS01s9clLU0+B3AWNm5vUElxkUrHFmmUmUrHFqmkuEgbtzfEXDIgXmk3J7r7i5JskNVXp7t9AFJja4fKx4/pt2zcmNFqbO2IqURAbmDGDiAAlWUlau/s6besvbNHlWUlMZUIyA2EGBCAlYur1NHVrbYT3TrlrrYT3ero6tbKxVVxFw2IFSEGBGDRjIm6d/kcTSotVvPxTk0qLda9y+cwOhEFL5Ih9gAyb9GMiYQWcBrOxAAAwSLEAADBIsQAAMEixAAAwSLEAADBIsQAAMEixAAAwSLEAADBIsQAAMEixAAAwSLEAADBIsQAAMEixAAAwSLEAADBIsQAAMEixAAAwSLEAADBIsQAAMEixAAAwSLEAADBIsQAAMEqirsAAID8Ulvfoo3bG9TY2qHKshKtXFylRTMmZmRfnIkBACJTW9+itTV7daStS+Xjx+hIW5fW1uxVbX1LRvZHiAEAIrNxe4NKiotUOrZIo8xUOrZIJcVF2ri9ISP7I8QAAJFpbO3QuDGj+y0bN2a0Gls7MrI/+sQA5LRs9q8gfZVlJTrS1qXSse/ES3tnjyrLSjKyP87EgDOorW/RXZt2afm3fqa7Nu3KWNs+3i3b/StI38rFVero6lbbiW6dclfbiW51dHVr5eKqjOyPEAOGwEE0XtnuX0H6Fs2YqHuXz9Gk0mI1H+/UpNJi3bt8TsbOnmlOFM0VGFzqQVRS39+N2xv4jGRBY2uHyseP6bcsk/0riMaiGROz9v0o+DMxfmljKNnupEZ/lWUlau/s6bcsk/0rCE/BhxjNFRgKB9F4Zbt/Bf2F0B9c8CHGL20MhYNovLLdv4J3hNJKVfB9YtkeDoqw9B5E+/eZzuIgmkXZ7F/BO0LpDy74EFu5uEpra/ZKSpyBtXf2JH9pz4q5ZMgVHERRiEIZVFPwzYk0VwDAu4XSH1zwZ2ISv7QB4HShtFIV/JkYAODdQmml4kwMADCgEFqpOBMDAASLMzG8C9NwAQgFZ2LoJ5QLHAFAIsRwGqbhAhASmhPRTygXOA6F5lCgcHAmFohsTcQZygWOg6E5FCgshFgAsnlgDn3CW5pDgcJCc2IAsjkRZ+gT3uZDcyhwtgqxKZ0QO0vZ/LBk+8AcwgWOg+GuBDiTfD3Q97bYlBQX9WuxycVZNqJEc+JZyHa/S+j9VNkUenNotoVw08Mo5XOfaaE2pWc8xMxsmZntN7M3zGxNpveXDdn+sHBgHr5Q5nvLBfl8QB9MPh/oz3SD33z9wZLR5kQzGy3pUUlLJR2U9J9mttnd92Zyv5kWR/NeyP1U2RZyc2g2hXLTwyjlc5/pUE3p+dzUmOk+scslveHu9ZJkZk9IWiEp6BCLo9+FAzOils8H9MHkc5/pULdOyecfLJluTpwq6UDK84PJZUGjeQ/5oBD7WvP5uztUU/qZmhpDZu6euY2b3SzpWne/M/n845Iud/fPprxmtaTVkjR27NgFc+fOzVh50tHc3Kzy8vK+5+2d3Wpp71JX9ykVF43SxHHFGjemcAZ7nl4fhS7E+mjv7FbTWyc0ykyjTDrl0il3TTl3bNqf5Vyuj7i+u3HWSWNrh7p7XKNSTltOnZKKRltsP1pGUh+7du16xt2XDbQu0//nDkqanvJ8mqRDqS9w9/WS1ktSdXW119XVnfXOMjl0trq6WumULd9QH/3lcn0M9b3I1Hcml+sjLnHWSWqfWGpTY5x9YiOsjwEDTMp8iP2npJlm9h5Jv5B0q6SPZWJH+dxxCZytM30v6GstDPk8OCyjIebu3Wb2GUnPSBot6bvu/lom9pWpjsveX6pvHD6uuzbtypsLI1EY8rlDHyOTrz9YMn6dmLv/q7vPcvf3uvv9mdpPJjouU6+jmXPlDQVxHc1wrV69Ou4i5JRcrY+4OvRztT7iRJ30F1V9ZHRgx0il0yd216Zd7xo623aiW5NKi/XYbQtyZptANvEZRp6wwVbkzbRTmRg6m8/DUlEY8nlI+dnK15krClXehFhU0w0dOHBAV155pWbPnq3/eOgT2rPlB5KkzuNvaevXP6Mff+lG/df6P9XRo0cz8Z+Rc06cOKHLL79c73vf+zR37lzdd999kqTW1lYtXbpUM2fO1NKlSwumPnr19PRo/vz5Wr58uaTcrY9sTcNVVVWlSy+9VPPmzVN1dbWk3KyTbE21dezYMd1000265JJLNHv2bG3fvj0n6yMb9u/fr3nz5vX9mzBhgtatWxdZfeRNc2JUmpqa1NTUpMsuu0zbXmnQimuW6AOffki/2LlFNna8Kq/8mCoPPqtS69RDDz0Ua1mzwd3V3t6u8ePH6+TJk/rABz6gb3zjG3rqqadUVlamNWvW6MEHH9TRo0cLoj56PfLII6qrq9Pbb7+tmpoa3X333QVdH1VVVaqrq9OkSZP6luVinWSreXXVqlVasmSJ7rzzTnV1damjo0MPPPBAztVHtvX09Gjq1KnasWOHHn300ZHUx6DNiXL3nPm3YMECzzVLrlnmy7/0bR9XMd1Xfmurb//fI37o0CGfNWtW3EXLuvb2dp8/f77X1tb6rFmz/NChQ+7uBVcfBw4c8Kuuusq3bdvm1113nbt7QdeHu/tFF13kzc3N/ZblYp1c980X/BPf3eGfenxn379PfHeHX/fNFyLbx1tvveVVVVV+6tSpfstzsT6y7ZlnnvH3v//97j7i+hg0N/KmOTETGhoa1LD/NW36s4+rqPNtbfjMtVo0Y6KmTJmiw4cPx128rOnp6dG8efNUUVGhpUuXauHChXrzzTc1ZcoUSSq4+vj85z+vr33taxqVMv1B3PURdz+PmelDH/qQFixYoPXr10uKv04Gko2pturr61VeXq5PfvKTmj9/vu688061t7fnZH1k2xNPPKGPfvSjkqL7fBBigzh+/LhuvPFGrVu3ThMmTIi7OLEaPXq09uzZo4MHD2rnzp169dVX4y5SbGpqalRRUaEFC3JnZF8u3FLlpZde0u7du7VlyxY9+uijeuGFF7K275HIxkCX7u5u7d69W5/+9Kf18ssva9y4cXrwwQcj236ourq6tHnzZt18882RbpcQG8DJkyd144036rbbbtNHPvIRSdLkyZPV1NQkKdFvVlFREWcRY3Heeefpiiuu0NatWwu2Pl566SVt3rxZVVVVuvXWW/XTn/5Ut99+e6z1kQv3yLrwwgslSRUVFbrhhhu0c+fOnPyMZGOgy7Rp0zRt2jQtXLhQknTTTTdp9+7dOVkf2bRlyxZddtllmjx5sqTojqmE2GncXXfccYdmz56tL3zhC33LP/zhD2vDhg2SpA0bNmjFihVxFTGrmpubdezYMUnSr371Kz333HO65JJLCrY+vvrVr+rgwYNqaGjQE088oauuukrf//73Y62PuC8FaW9vV1tbmyTp315t1Hd+8M968n+lcy5eqL/8q7+WlFufkUUzJuqx2xao5rNL9NhtCyIfqXnBBRdo+vTp2r9/vyRp27ZtmjNnTsF+Z3r98Ic/7GtKlKI7pjI68TQvvviilixZoksvvbSvz+OBBx7QwoULdcstt6ixsVGVlZX60Y9+pLKysljLmg2vvPKKVq1apZ6eHp06dUq33HKLvvKVr6ilpaUg6yPV888/r4cfflg1NTWx1kfcFzTX19frhhtuUEdXj355rF1Vl39I1dffodaWVr20/ss6p+uoZr33PQX1GdmzZ0/fyMQZM2bo8ccf7/v+hPSdiWqC6I6ODk2fPl319fU699xzJWmk35lBRycSYkDgcmWG8rjDFNHKlc9VUv7P2AEUqmxd0HwmcTdrIlq50Nc6HIVzF0cgj+XCDOWVZSXvOhPL9ztF57PG1g6Vjx/Tb1ku/ijhTAxAJJinMb9k45q6KBBiQB6I+2JnKXeaNRGNUH6UMLADCFyOdcAjj0Q1OjECgw7soE8MCBx3b0am5EJf65nQnAgEjlGBKGSEGBC4UDrggUwgxIDAhdIBD2QCIQYEjlGBKGQM7ADyQAgd8EAmcCYGAAgWIQYACBYhBgAIFiEGAAgWIQYACBajE4EY5NCcdEDQOBMDsqx3wt4jbV0qHz9GR9q6tLZmbywzzwOhI8SALAvljrlACAgxIMuYsBeIDiEGZBkT9gLRIcSALGPCXiA6hBiQZUzYC0SHIfZADJiwF4gGZ2IAgGARYgCAYBFiAIBg0ScGYEBMjYUQEGJADHI9IHqnxiopLuo3NRajKJFraE4EsiyEuROZGit31da36K5Nu7T8Wz/TXZt25dTnJg6EGJBlIQQEU2ONTLaCJYQfQNlGiAFnEPUBKoSAYGqs4ctmsITwAyjbCDFgCJk4QIUQEEyNNXzZDJYQfgBlGyEGDCETB6gQAoKpsYYvm8ESwg+gbGN0IjCExtYOlY8f029Zugeo3oDoPzpxVs4FBFNjDU9lWYmOtHWpdOw7h9NMBcvKxVVaW7NXUuJz2N7Zk/wBNCvyfYWCEAOGkKkDFAGRP7IZLKH8AMomc/e4y9Cnurra6+rq4i4G0Cf1eqnUAxRNa0iV69f95QEbdAUhBgyNAxQQu0FDjOZE4Axo+gNyF6MTAQDBIsQAAMFKK8TM7Otm9t9m9oqZ/bOZnZey7h4ze8PM9pvZtekXFQCA/tI9E3tW0m+4+29K+h9J90iSmc2RdKukuZKWSXrMzEYPuhUAAM5CWiHm7j9x9+7k01pJ05KPV0h6wt073f3nkt6QdHk6+wIA4HRR9ol9StKW5OOpkg6krDuYXAYAQGTOOMTezJ6TdMEAq77s7k8nX/NlSd2SNvW+bYDXD3hBmpmtlrRakiorK4dRZAAAEs4YYu5+zVDrzWyVpOWSrvZ3rpw+KGl6ysumSTo0yPbXS1ovJS52HkaZAQCQlObFzma2TNKXJH3Q3VNnRN0s6Qdm9oikCyXNlLQznX2h8DBTBoAzSbdP7NuSSiU9a2Z7zOw7kuTur0l6UtJeSVsl/ZG79wy+GaA/7mALYDjSOhNz94uHWHe/pPvT2T4KV+p9vCT1/d24vYGzMQB9mLEDOYk72AIYDkIMOYk72AIYDkIMOWnl4ip1dHWr7US3Trmr7UR38kaDVXEXDUAOIcSQk3rvYDuptFjNxzs1qbSYG1ECeBfuJ4acxX28AJwJZ2IAgGARYgCAYBFiAIBgEWIAgGARYgCAYBFiAIBgEWIAgGARYgCAYBFiAIBgEWIAgGCZu8ddhj5mttXdl8VdDgBAGHIqxAAAGAmaEwEAwSLEAADBIsQAAMEixAAAwSLEAADB+n8F5PUc3apNnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plot setup\n",
    "plt.figure(figsize=(7,6)) \n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "\n",
    "# load data from a csv file\n",
    "data = np.genfromtxt('files/blood_pressure.csv', delimiter=',')[1:,:]\n",
    "age = data[:,0]\n",
    "sbp = data[:,1]\n",
    "\n",
    "# the least squares line\n",
    "meanx = np.mean(age)\n",
    "meany = np.mean(sbp)\n",
    "varx = np.var(age)\n",
    "cov = np.inner(age-meanx, sbp-meany)/age.shape[0]\n",
    "\n",
    "res = sbp - cov*(age - meanx)/varx - meany\n",
    "\n",
    "plt.gca().spines['bottom'].set_position(('data',0))\n",
    "plt.title(\"Residuals\")\n",
    "plt.scatter(age, res, zorder=0, alpha=0.7);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MLE linear regression for synthetic data\n",
    "\n",
    "True valus are $a=-2, b = 1$ and $\\sigma^2=0.04$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAFzCAYAAACq+qpxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXRT1doG8GdnbNOmc4ECAiIOgDJIFYoKKg5cLjhfRdQ6o1YRnHAAUbnApyAOFwUFRxSH64SKgiKCDLdVW0SG4lgRIWlpU9qmDU2aZH9/nFYaaGkLSU9O8vzWyoKmJ8l7GiVv99n72UJKCSIiIqIGOrULICIiovDC5oCIiIgCsDkgIiKiAGwOiIiIKACbAyIiIgrA5oCIiIgCGNQuoC1GjhwpV6xYoXYZREQUXoTaBUQaTY0clJWVqV0CERFRxNNUc0BEREShx+aAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKoGpzIIRIEkK8L4T4SQixXQiRpWY9REREpP7eCs8CWCGlvEwIYQJgUbkeIiKiqKdacyCESAAwDMB1ACCl9ADwqFUPERERKdQcOegJoBTAq0KI/gAKAEyUUtY0PkgIMR7AeADo1q1buxcZLHlFDizO3YGd5S50S7EgO6sHhvRMVbssIiKig6g558AA4GQAC6SUAwHUAHjgwIOklAullJlSysz09PT2rjEo8oocmLGsEGVOD9LjzShzejBjWSHyihxql0ZERHQQNZuDXQB2SSm/rf/6fSjNQsRZnLsDFpMB1hgDdELAGmOAxWTA4twdKldGRER0MNWaAyllMYC/hBDH1981AkChWvWE0s5yF+LM+oD74sx67Cx3qVQRERFR89RerTABwJL6lQpFAK5XuZ6Q6JZiQZnTA2vM/h93jduHbilcnEFEROFH1eZASrkJQKaaNbSH7KwemLFMGRSJM+tR4/bB5fEiO+s4lSujUOIkVCLSKiGlVLuGVsvMzJT5+flql3FYQvlBwQ+h8NMwCdViMgQ0hFNH9+F7QxR8Qu0CIg2bA43jh1B4yllScNClJGetF2lWE+ZfNUjFyogiEpuDIOPeChrHlRDhiZNQiUjLoq85qKtTu4Kg4odQeOqWYkGN2xdwHyehEpFWRFdz4HQCffsCs2ZFTJPAD6HwlJ3VAy6PF85aL/xSwlnrrZ+E2kPt0oiIWhRdzYHHAwwYAEyZApxyClBQoHZFRyzUH0J5RQ7kLCnA6HnrkLOkgKmOrTSkZyqmju6DNKsJpdVupFlNnAdCRJoRnRMSly4FcnKAPXuAe+4BHn0UiI098udVSahWK3CyIxFpBCckBll0NgcAUFEB3Hcf8NJLwLHHAosWAcOHB+e5IwRn3BORRrA5CLLouqzQWFKS0hCsWgX4fMCZZwK33gpUVqpdWdjgZEciouikdnyy+s4+G9iyBXjkEeCpp4Bly4AFC4AxY9SuTHWMfSaicPPLL0B1deB9Y8bgVHWqiUjO3buxnc0BAFgswJw5wOWXAzfeCFxwATB2LPDss0CHDmpXpxrGPhNRuKmuBhITD7q7QoVSIlUSEM2XFZpyyilAfj4wfTrw4YdAnz7Am28CGpqXEUyccU9EFJ04cnAgkwl4+GHg0kuBm24CrrkGeOst4IUXgG7d1K6u3Q3pmcpmgIgoynDkoDl9+gDr1gH/+Q+wdq0SnvT884Dfr3ZlREREIcXm4FD0emDCBGDbNuC004A77gCGDQN++kntyoiIiEKGzUFrdO8OLF8OvP46UFgI9O8PzJwZMRHMRESRqqzs3OFlZecyxKaN2By0lhBAdjawfTtw0UXA1KlAZqYygZGIiMKOy/VGssezfpLHs36Sy/VWktr1aAmbg7bq2BF4910lgrmsDBg8GJg8GXAxGIiIKJxUVU2eYDafPdtsHj63quqeO9WuR0u4WuFwXXihErc8ebKSkfDhh0oU85lnql0ZaVCo9scgimadOtmnN/pyvWqFaBBHDo5EUhKwcCHw9dfK12edBdxyCyOYqU0aNrgqc3qQHm9GmdODGcsKuQMmkcbU1W03Fhd3mGmzmVbbbPqNNlvsUofjvGHNHb9nz4lX2WyWD2w23dbi4vTH2/J9m033Q+BNbC8u7vxww/ft9rQ5NptxvVKH+Ys9ewb8qy3nwuYgGM46C9i8Gbj3XmX0oE8f4JNP1K6KNGJx7g5YTAZYYwzQCQFrjAEWkwGLc3eoXBkRtYXf7zAIEVdstU69ulOnikEm09Bn3O5Vzzid07s0dbxOl7LHZDplgU7X4f22fr9zZ//Ahltq6tqhgKg1mQYvb/i+xXLli+npm8/q3Nl3clzc7bd5vVsn7d17Xd/Wngubg2BpiGD+9lsgLU257HDFFUBJidqVUZjjBldEwVdT83KK3Z60wGYz/M9m02+025NeqK1dERfK1zSbT9/XseMf86zWabt1OqtMS1u1BjDtcrtXNvmhnJa2dmVa2jdfCRHTZPxzS99vUFU16XzAUJ6c/M7fM+QTE+f9ZjT2rl9Sp5OAkF7vL61O8mNzEGwNKxhmzFAmLfbpAyxeHLURzNSybikW1Lh9AfdxgyuiI+Pz/R5vMg19Iz19y/Dk5LfOAuqSKysnjm3qWLs96UWbTZ/f1M1uT3rxcGuoqVmYCriPNhpP+e3wz6RlXu/PF+v1XZYKYQ64v7g44xGbTfdjTc3cFYChNCFh1jetfU5OSAwFoxGYMgW45BIlgvnaa5UI5hdfVDITiBrhBldEwZeQMGsngJ0AYDT2rqyqemADUHvwlk0AMjIqbgn263u9fxqqqu6Zq9N1+Cgx8amiYD9/A6fz/zKkrD7VYrlnyoHf69TJ/pjf7/h3RcUtA+vqfjhVrz/a09rn5chBKPXurUQwz5sHbNigRDDPmwf4fC0/NsTyihzIWVKA0fPWIWdJASe/qehIN7jie0l0sNLS00ba7XFv22yGXJtNn+/z7Riv06X9EbznHzymYTKg3Z64qPH3/H6nKC09eQ6g86SlfTu9uecIBpfrpYuEiCuwWh/d1dT3dbpUf0rK+wVS1nQqL//Hla19Xo4chJpOp8QujxkD3HorcOedwNtv75+4qIKG2fEWkyFgdjx3XFTP4W5wxfeS6GAOx5ghdXXf3xcbe8WkhITZhQBQUtJttdE4cHtTx9vtiYukdGY29T0hrPkZGZU3H3h/evq3nwL49MD7pXRjz55jZgGe1NTUNTcbDN29R3g6h+Tz7brIaOy7sOUjpd7vr+Scg7DTvTvw+efAG28AP/8MDByozEvwtHqUJ2g4Oz5y8L0kOpjP99sJQpjsFsu1RW73Fwmlpf1mAd6U+PgHfm/q+IyMypsbz/5vfGuqMTiUkpLuj0lZc0xKyqe3mkyD3Ieu066vq9tsAvw6QOrq6jabfD67vrXfLy+/fCBQ1zEhYfaKxs9bU/NySmlp1ii3e43F73foHI6Rp/v9paMNhmPzWnsebA7akxDA1VcrEcyXXKJsDZ2ZCXz/fbuWwdnxkYPvJdHB4uLu/QSQBofjvA2VlTkLhUj4E4j53WDoFdINcZzOmZ39/pKxUu7r7XCcvaHhskNp6eAxgDJCUVLS6+/5DWVlQ3NKS/tv8fn+usXvd1xYWtp/S1nZ0JzWft/jWXuxTpfypdl8Tk3jOoTQS693yziHY8Ta4uIO+W73mvsNhn6z0tLWrGrtuQipoVn0mZmZMj+S9jL45BPgttuA4mLgrruA6dOVJZEhlrOkAGVOD6wx+68qOWu9SLOaMP+qQSF/fQoevpcUbTZuBBIPmFY4bBiOV6eaiJS0eze+48iBmi64QNnl8eabgblzgZNO2p+2GELZWT3g8njhrPXCLyWctd762fE9Qv7aFFx8L4koFNgcqC0xEXjhBWD1amXy4ogRSrNQccjMiyNypLPjKXzwvSSiUOBlhXCybx/w6KPAk08quz/On69sD01ERAB4WaEd8LJC2ImNBZ54AvjuO6BDB+Dii4HLL2cEMxERtSs2B+Fo0CBlBcPMmcqkxd69gddfZwRzPYb+EBGFFpuDcGU0Ag89BGzapIQlXXcdMHIksGOH2pWpitsbE1FblZWdO7ys7NzhatehJWwOwt0JJwBr1wLPPQf873/AiScC//lPWEQwq4GhP0TUFi7XG8kez/pJHs/6SS7XW0lq16MVbA60QKcDbr8d2LYNGDYMmDgROP10ZRlklGHoDxG1RVXV5Alm89mzzebhc6uq7rlT7Xq0gnsraEm3bsBnnwFLlgCTJikRzFOnAvffD5hMalfXLrqlWA4K/eH2xkTUnE6d7I03PlqvWiEaw5EDrWmIYC4sVCKYp03bP4ExCjD0h4go9NgcaFWHDsrujp98AuzdCwwZAtx7L+CK7OF1hv4QEYUeLyto3ZgxyjyE++9XIpg/+ghYtAg4+2y1KwuZw93emIiii81m/CY29urbkpNfjb4JWkeIIweRoCGCec2adotgJiIKZ/v2LbUC3g5xcbcWhfZ13k2025Oft9l0m2w20+rS0sGjmzvW6ZzexW5PXGiz6b+32QwbioszpjXegrlBVdUD3W023Ra7PW1Ow32VlZN62u3W1202fYHNZl5ZVjbsnAMfV1qaNcpmi1leX8tX5eWXHvbua2wOIsnw4cDmzcoowquvKvkIS5eqXRURUbvbt++N4wCjzWQaXBvK16moyHkE0HlSU78eGhNzwb11dfmPVVZO6NXUsdXVcx8RwlSenv7DaYmJ8y70+/ee4nCMGHfgcTU18x8RwrKl4Wufz66vqVmwQK/vsqZTpz2nmM3nPezxrH+ysnJyj4ZjHI5RQ+vqCu6Ljb34gU6dKgcmJDwxzmQ6/a/DPS82B5EmNhZ4/HElgrljRyWC+V//UraFJiKKEl7vL8cLEfNHcXHXB5Tf1M1flpdfFtR9zN3u9bFSlp9nsdz4rNl8pisl5f0CnS55VW3txxc2dbyU7qMMhgHLjcZ+nri428p0uvR1fn/5sY2PKS3NGgUYqnS69NyG+5zOGT2Bug7p6T++qtOl+lNTP80TIn5jbe1//34dj+ebOw2Gvs8nJ7/9o05nlfHxd+2Jj79rz+GeG5uDSHXyyUqDMGsW8OmnyijCa68xgpmIooLfX3q8lDX9DYZj8jt0KMrS67t8Ulu7bEZzx9vtSS/abPr8pm52e9KLTT1m375XewDCn5g4e0fDfUKk/uT3Vx3b1PEGw3Gve72bRnk838ZUVz/dwe8vHWYwHLe24fu1tSvi6uo2ToyPf+DxA85GNPF0wu+vPE45V4dOSteJUlan2GzmlTabcW1xccY0j6fAfIgf0SGxOYhkRiPw4IPAjz8CffsC118PnH9+1EcwE1Hkk9J5nF7f/bW0tG++Mhi6e+PiJvwXcB/d1DV+AMjIqLilc2dfZlO3jIyKW5p6jN9fYQF0zsb3CRFbDfjimjrebD77O7+/5tiysiEbq6ruXqfTWbempq78quH7FRU3TtLru79vtU4OGOqNj3+gCDCU79lz0k1e758Gh+Mfp0npPAXwxQCAy/VmGgCjz7f7/MTEueOSkl65UMqq3nv3XpLTxh/b39gcRIPjjwe++UbZAjo3V2kUnn02aiOYiSjySVl7XEzM6BUNX3u921MAnVOvzwjaP3w6XZIL8McH3lsbB+hrDjzW73eKmprnX9HrO3+Znl7YPylpyWAp6xJKSnrdBwAVFTed4Pc7hqamfv7agY81GLp7LZYbcvx+25l79vTc4PFsuEGI1OVCxBQDgF7ftRYAjMY+b8TF3VFqsVyz12g8+VWfb89h7yfBpYzRQqcDbrsNGD0auPVWJWHx7beBl19WmoUIlFfkwOLcHdhZ7kK3FAuys3pwCSRRFHA6Z3YG/PFG48nlDfe53SvP1emSVzf3GLs9cZGUzsymvieENT8jo/LmA++Pjb1+h8v1ir6q6oHuCQmP/wkAfr/jBJ0u4dcDj62tXZoEeDMSE59502jsXWc09q5wuV780OP5dhKAOR7PhsGAp8uePSesUR7htwDQ22yxvTp33ndxUtILPyclvXD1/nrj3tHpen2k1HFp1d69hqBOLGNzEG2OOgpYtkxpDO68U4lgnjJFufzQKIJZ6x+sDbs3WkyGgN0bGZhEFPnc7q+PA+Ctrp49Jibm0jfKyy8Z5vP9NTY+/t6DVgY0aOrDvyVm8+n7hEhZWVOzaKLZPHpKTc2zvf3+8nPi4nKuOPBYi+WavRUVN+6qrLxrnNE46GWP51tLXd3mi4WI/wkAkpJee7eubuNnDcc7nTNvlNLVJTHxqUcAoKLi1uPj4nL+kHKfbu/e7HFS1qUnJ7/zYcPxen3XD+rqCq+uqXl5rU6X5K2r++E6vb5js81QS3hZIRoJAYwbB2zfrqxkePRRJYL5u+8ARMa2yNy9kSh6+Xw7j9fp0pf6fMUnFxcn5Xs8GybExl5+W8Nv98GUlDT/UcBvdjiG5dbWLn3KaMx8JDFx3m+AMhpRUtLr7/kKFkv27X5/8RklJV3z9u69bCUgvAkJT8wCAJNpcG1c3G1lDTchjC5A57ZYrtsLALW1n11YWjpwQ1lZVq7fX5wVH3/v9UZj77qG505NXTNfp0vcUlk5/su9e69YLkRCYUrKJwsO97yE1NDs9czMTJmfn692GZHns8+USw02GzBxIib1vQS2On3A5kbOWi/SrCbMvyqoK4FCZvS8dUiPN0Mn9k/y9UuJ0mo3lk04Q8XKiOhIbNyo5L41NmwYjlenmoiUtHs3vuPIAQH//KeyHfQttwBPP40H7r8cp/xWEHCI1rZF7pZiQY07cN4Rd28kImodNgekSEhQVjN88w2E0Yj75k7Ada/8G5aaKgDa+2Dl7o1ERIePzQEFGjYMO1ZtwHvnXIWhGz7Hv6dcgd65X2nug5W7NxIRHT7OOaAm5RU58PXbK3D5C4+h165f4Th/NFJfXQhkZKhdGhFFMc45CLnwmHMghNALIX4QQixTuxbab0jPVDw05Sr0KtoGPP44UtesVCKYX3mFEcxERBFO9eYAwEQA29UugpphNCq7PP74I3DSScCNNwLnnQcUhXQXVCKioCkrO3d4Wdm5h50WGI1UbQ6EEF0B/BPAS2rWQa1w/PHAmjXAggXAt98qjcLTTzOCmYjCmsv1RrLHs36Sx7N+ksv1VhIAlJePHWC3x71rt8e/abenPuX1/slAwAOoPXLwDIDJAPzNHSCEGC+EyBdC5JeWlrZfZXQwnU7JQ9i2DTjrLODuu4GhQ4GtW9WujIioSVVVkyeYzWfPNpuHz62quudOADCZsmypqWuzMzKqr9bp4v+qqLhqhNp1hhvVuiUhxGgAe6SUBUKIM5s7Tkq5EMBCQJmQ2E7l0aEcdZSyDfTbbwMTJyrbQz/0kBLBbD7sHUKJiIKuUyf79EZfrgeA+PiJe/bfpfMCumZ/QY1Wao4cnAbgAiHEDgDvADhbCPGmivVQWzREMBcWKhHMjz2mNAl5eWpXRkTUKk7no119vpLhSUmvrAnm89bVbTcWF3eYabOZVtts+o02W+xSh+O8Yc3XMb2L3Z640GbTf2+zGTYUF2dMa7y1tM2m+yHwJrYXF3d+uDWPBYDS0qxRNlvMcptNt8lmM31VXn5pi1G3qjUHUsoHpZRdpZQ9AIwF8LWU8uoWHkbhJj0dWLJE2cypqkq5zHDXXUDNQTuWEhGFjdraFXHV1U8+EReXc5/B0Kuu5Ue0nt/vMAgRV2y1Tr26U6eKQSbT0Gfc7lXPOJ3TuzR1fHX13EeEMJWnp/9wWmLivAv9/r2nOBwj/t4kqnNn/8CGW2rq2qGAqDWZBi9vzWMdjlFD6+oK7ouNvfiBTp0qByYkPDHOZDr9r5bOQe05BxQpGiKYb7sNeOYZ4MQTga++UrsqIopSXu+fhpKS7nfZbKavbTaxzWYTPyu32I99Prt+794rnjaZzpqXmPjkH8F+bbP59H0dO/4xz2qdtluns8q0tFVrANMut3tl36aOl9J9lMEwYLnR2M8TF3dbmU6Xvs7vLz+2qWOrqiadDxjKk5PfyW/NYz2eb+40GPo+n5z89o86nVXGx9+1Jz7+rj1NPXdjYdEcSCnXSClHq10HHaGEBOD554G1a5Xtn889F7jhBmDvXrUrI6IoU1Z2xiS/vyzLap0xLjV1XaYQ8blCJK2Mj7/v9vLyC0ZLWdPP41lzu91ufaO0dOio5p7Hbk960WbT5zd1s9uTXmxNLTU1C1MB99FG4ym/NfV9g+G4173eTaM8nm9jqquf7uD3lw4zGI5b29SxXu/PF+v1XZYKYW7xsX6/Qyel60Qpq1NsNvNKm824trg4Y5rHU9Di5DAu34gAeUUOLM7dgZ3lLnRLsSA7q4e6McFnnKHkIkyfDsyeDSxfDjz3HHDpperVRERRw+3+Ks7v35UdHz95jNU6uRgA9PoeX/h8O0YlJEzfBUzfBeDj1jxXRkbFLS0f1Tyv909DVdU9c3W6Dh8lJj7VZECM2Xz2dzU1Cy8vKxuyEYBep0v7KDV15UFDr07n/2VIWX2qxXLPlNY81uV6Mw2A0efbfX5i4txxQiR6Kytvnb937yU5HTv++fSh6g6LkQM6fHlFDsxYVogypwfp8WaUOT2YsawQeUUOdQuLiQFmzQLy84HOnYHLLlOaA7td3bqIKOI5nY9nAqa/EhIe/7PhPin3JQCmsvasw+93itLSk+cAOk9a2rfTmzumpub5V/T6zl+mpxf2T0paMljKuoSSkl73HXisy/XSRULEFVitj+5qzWP1+q61AGA09nkjLu6OUovlmr1G48mv+nx7WgyEYnOgcYtzd8BiMsAaY4BOCFhjDLCYDFicu0PlyuoNGKCEJj3+OPD554xgJqKQk3JvihCGyv1fu+Hz2c41GI5Z3dbnstsTFx28WkC52e2Ji5qvwY09e46ZBXhSU1O/nmAwdPc2dVxt7dIkwJuRmPjMm0Zj7zqLZVyF0dj/Q7+/9KAPcJ9v10UGw3EftfaxsbGXVgGG4raeM8DmQPN2lrsQZw5YtYI4sx47y10qVdQEg2F/BHO/fkoE87nnMoKZiELCYOj7q5SuvhUVN53g8RSY9+w59h4ASE5+7/O2PldGRuXNjVcLNL5lZFTe3NzjSkq6PyZlzTEpKZ/eajINcjd3nMVyzV7AuKuy8q5xPp9dv2/fUmtd3eaLhYj/qfFx5eWXDwTqOiYkzF7Rlsfq9V0/qKsrvLqm5uWUffs+SKir++E6vb5ji00SmwON65ZiQY07MMK4xu1DtxSLShUdwnHHAatXAy+8AHz3HSOYiSgkkpMXb9Xrey5wuV5fVFY2ZJWU+9KTkl68ubnf3oPN6ZzZ2e8vGSvlvt4Ox9kbGkYaSksHjwGU0YiSkl5/z2WwWLJv9/uLzygp6Zq3d+9lKwHhTUh4Ylbj5/R41l6s06V8aTafE7BOvKXHpqauma/TJW6prBz/5d69VywXIqEwJeWTBS2dA7ds1riGOQcWkwFxZj1q3D64PF5MHd1H3UmJLdm1S1n2uGwZcOqpwMsvK8sfiYgOgVs2h1x4bNlMR2ZIz1RMHd0HaVYTSqvdSLOaVG0M8oocyFlSgNHz1iFnSUHzEyO7dgU++USJYC4qUtIVH3kEcDc7+kZERO2EIwcUNIc9ilFWpqQqvvmmMmHx5ZeBIUPar3Ai0gyOHIQcRw4ouA575URaGvDGG8pqBqdTiWCeNAmorm6HqrWr1aM0RERtxOaAguaIV0784x9KBHNODvDss8ochC+/DEGl2he2+RZEFBHYHFDQBGXlhNWqpCmuW6cEKZ1/PnDddUB5eXCL1biwz7cgIk1jc0BBk53VAy6PF85aL/xSwlnrhcvjRXZWj7Y/2emnA5s2AVOmKLs+9u4NvPcew5PqaSLfgijMlZePHWC3x71rt8e/abenPuX1/sktBeqxOaCgCfrKiZgYYMYMJYL5qKOAyy8HLrkEsNmCW7gGaSrfgihMmUxZttTUtdkZGdVX63Txf1VUXDVC7ZrCBbskCqohPVODv4yyf38gL08JTJo2TVnR8OSTStKiEMF9LY3IzuqBGcsKASBgZUh21nEqV0akHfHxExttXazzAjq/etWEF44ckDYYDMB99wFbtgADBwI33wyMGAH8/rvalaki3PItiMLNvn0fJNhs4uf6dMIfbTbT6j17Bl7W1LFO56Ndfb6S4UlJr6wJbg3vJtrtyc/bbLpNNptpdWnp4NHNHWu3W9+w2XRb9u/dYF7Rlufas+fEq2w2ywc2m25rcXH64w3319VtNxYXd5hps5lW22z6jTZb7FKH47xhLdXO5oC0pVcvYNUqYOFCoKBAiWCeOxfwtksqalgZ0jMV868ahGUTzsD8qwaxMSBqxOV6vTeg31u/D0J/o3HgXK9303SX643kxsfV1q6Iq65+8om4uJz7DIZedcGsoaIi5xFA50lN/XpoTMwF99bV5T9WWTmhV3PHGwz9pu/fu8E9si3PpdOl7DGZTlmg03V4v/Hj/H6HQYi4Yqt16tWdOlUMMpmGPuN2r3rG6Zze5VC1szkg7dHplJGDwkJlA6d77wWysoDNm9WujIjChM/3a28hLNsavo6NHfsdAL3Xuy1h/zF2/d69VzxtMp01LzHxyT+C+fpu9/pYKcvPs1hufNZsPtOVkvJ+gU6XvKq29uMLQ/FcaWlrV6alffOVEDEVjR9rNp++r2PHP+ZZrdN263RWmZa2ag1g2uV2r+x7qNdkc0Da1aULsHQp8O67wJ9/AoMGAQ8/zAhmIoLfX9Zbp0veBgD79i21VlfPugeI3Wq1PvZnwzHl5ReMlrKmn8ez5na73fpGaenQUU09l92e9KLNps9v6ma3J73Y1GP27Xu1ByD8iYmzdzTcJ0TqT35/1bHN1ez1br3HZjPk2e1xb5eXX3jqkTxXc2pqFqYC7qONxlN+O9RxnJBI2iaEsophxAglgnnGDOCDD4CXXlKSFokoKvn91b0Bx/k2m/5qwB8nhHVdUtKLNwlh/vuY9PTvPwbwcUvPlZFRcUtLxxz8+hUWQOdsfJ8QsdWAL66p42NixsyxWG74Xafr4qmouOaftbXLXqiqmnphQsKMv9r6XM3xev80VFXdM1en6480VJsAACAASURBVPBRYuJTRYc6liMHFBlSU4HFi4Hly4GaGiUn4c47GcFMFIXq6rYbgdpj4uMfGtO5s+9kk+n0CVLWDNDpkoM6p+BQdLokF+CPD7y3Ng7Q1zR1fHLyW5vN5nNqjMbedenp+UuFiNtYW/vR8MN5rqb4/U5RWnryHEDnSUv7dnqL9bf2iYk0YeRIYOtW4PbblaTFE08EvvhC7aqIqB1VVz9xHCA8VuvDfwFAWtq6LwGTzemcdv7hPJ/dnrho/yqCwJvdnrioqcfExl6/A5D6qqoHujfc5/c7TtDpEn5t3asKCUgRjOeS0o09e46ZBXhSU1O/nmAwdG9xBjebA4o8Viswbx6wfj0QG6s0DNdeCzi47wBRNKir29RHiNhfG19C0OvTv/F6/zj7cJ4vI6Py5v2rCAJvGRmVNzf1GLP59H1CpKysqVk00e1eH1te/q+T/f7yc2JiLjzoMsa+fUutDsfI0+vqNpt8Pru+tHTwGCmdp5jNo9a39rl8Pru+rm6zCfDrAKlreC4AKCnp/piUNcekpHx6q8k0qFWTstgcUOQaOhT44Qclgvmtt5TwpP/+lxHMRBHO7y/pLUT8z43vMxoHrpOy8jTlA7R9JCXNfxTwmx2OYbm1tUufMhozH0lMnPf3REC7PXFRSUmvW6SsMLjdayeVlg7MKynp+q3Xu/Uas3lkTuMVFC09V1nZ0JzS0v5bfL6/bvH7HReWlvbfUlY2NMfpnNnZ7y8ZK+W+3g7H2RsaRjxKSwePOVTtQmroH8rMzEyZn5+vdhkRL6/IgcW5O7Cz3IVuKRZkZ/XQ/hr6H39UEhULCoALLgDmz1dWOxCRpmzcCCQmBt43bBiOV6eaiJS0eze+48gBBYjYrYAbIpjnzFG2ge7TRwlS8jMtlYjoQFzKGIGO5Df/xlsBA/j7z8W5O8Jq9OCwztFgUAKTLrpICVG65Rbg7beBRYuU5MUIF5EjQkQUEhw5iDBH+pu/FrYCPuLRjV69gK+/VkYONm5UIpjnzInoCOaIHREiopBgcxBhGv/mrxMC1hgDLCYDFufuaNXjtbAV8JGeIwAlPKkhgvm884DJk4EhQ5S5CREoKD8zIooabA4izJH+5p+d1QMujxfOWi/8UsJZ663fCrhHCKo9PEEd3WgcwfzXX0BmJjB1KlBbG6Rqw4MWRoSIKHywOYgwR/qbvxa2Ag766EZDBHNhITBuHDBzprIt9IYNQag2PGhhRIiIwgebgwgTjN/8w30r4JCNbqSmAq+/DqxYAbhcwBlnABMmAE5ny48Nc1oYESKi8MGcgwgUDbPSQ36O1dVKeNK8ecBRRwEvvqgkLWpYNPx3QZHvl18O3jJlzBgMVqeaiOTcvRvb2RwQHUpurhKetH07cM01wNNPKyMMRBROhNoFRBpeViA6lKwsJYL54YeVTITevZXJixpqqomI2orNQZTIK3IgZ0kBRs9bh5wlBVzf3hZmMzB9uhK93L07MHasEqS0e7falRERhQSbgyjAAJwg6ddPuczw5JPAypWajmBms0hEh8I5B1EgZ0kBypyev6OQAcBZ60Wa1YT5Vw1SsTIN+/13YPx4JWnxzDM1FcHc0CxaTAbEmfWocfvg8njDbslqY4eaTMmJlgTOOQg6jhxEAQbghMAxxwBffQW89JIyJ6GVEczh8Bu71tISDzXyxVExotBgcxAFGIATIkIoKxkKC5Vlji1EMIfLB5nWmsXGzUzlvjrscNSgqKwGk975AU+v/FlTjQ6RVrA5iAIMwAmxzp2BDz8E3ntvfwTzlCkHRTCHy2/sWmsWG5qZvS4Pfi2phsfrR4xBh6paL37cVQmPN/BcwrnRIdIKNgdRQAuRyMGi2rC9EMBllyl5CFdfDcyaBQwYAKxf//ch4fIbu9aaxYZmZtfefdDrBAw6Ab+UiDcZYDHq8ecBP79wbnSItILNQZQI90jkYAiLYfuUFODVV4EvvgDcbiWC+Y47AKczbH5j11qz2NDMVLu90AnA6/fD5we6JMege6oF++p8mml0iLSCqxUoYoTdqozqamWHx//8B+jaFT89Mhv31nTR1CqBcJFX5MCkd35AVa0X8SYDuiTHICXODGetF0IAqfEmrlaIblytEGRsDihijJ63DunxZujE/n8n/FKitNqNZRPOUK+w3FzgppuAwkKUXngZZo+8FdvrTPwgayMtLsGkdsPmIMh4WYEiRrgM2x8kKwvYuBGYNg3pny3FnIfHYlnaLswfdzI/1NpAa5dDiLSMIwcUMTTxm+WWLcryx++/B0aPBhYsALp2VbsqIq3jyEGQceSAIoYmfrM86STlMsPcucCqVUoE8wsvaDKCmYgiF0cOiNTSOIJ52DAlgvm449SuikiLOHIQZBw5IFJL4wjmH39UNnZ64okWI5iJiEKNzQGRmhpHMI8aBTzwAHDqqcp+DUREKmFzQBQOGiKY338fsNmAU04BHnwQ2LdP7cqIKAqxOSAKJ5deqowiZGcDjz+uRDCvW6d2VUQUZdgcEIWblBTglVeUCGaPR5msePvtQFWV2pURUZRgc0AUrs47T8lFmDRJyUPo2xf47DO1qyKiKMDmgCicxccDTz8N/O9/QEKCEpx01VVAaanalRFRBFOtORBCHCWEWC2E2C6E2CaEmKhWLURhb8gQJYL5kUeA995TwpPefhvQUE4JEWmHmiMHXgD3SCl7AxgC4HYhRB8V6yEKb2Yz8OijSpNwzDHAuHHAmDHAX3+pXRkRRRjVmgMppV1KubH+704A2wF0UaseIs048URgwwblcsPq1cpchAULGMFMREETFnMOhBA9AAwE8G0T3xsvhMgXQuSX8jorkUKvVyYqbt0KDB4M5OQAZ54J/PKL2pURUQRQvTkQQsQD+ADAJCnlQWu1pJQLpZSZUsrM9PT09i+QKJwdfTTw5ZfK0sctW5QI5scfB+rq1K6MiDRM1eZACGGE0hgskVJ+qGYtRJolBHD99Up40ujRSrLi4MGMYCaiw6bmagUB4GUA26WUT6lVB1HEyMhQ4pc/+ACw2xnBTESHTc2Rg9MAXAPgbCHEpvrbKBXrIYoMl1yijCJce61yiaF/f2DtWrWrIiINUXO1wnoppZBS9pNSDqi/fa5WPUQRJTkZePllYOVKZQvo4cOVSYuMYCaiVlB9QiIRhdA55ygTFe+6C3jxRUYwE1GrsDkginRxccBTTykRzImJyqTFceMYwUxEzWJzQBQtBg9W0hUfe0yZuNi7N7BkCSOYieggbA6IoonJBEybpixzPPZY4OqrlZEERjATUSNsDoiiUd++wPr1wDPPAGvWKBs5zZ/PCGYiAsDmgCh66fXAxIlKBHNWFnD77cqqhp9/VrsyIlIZmwOiaHf00cAXXwCvvQZs26bkIsyaxQhmoijG5oCIlAjma69VwpPGjAGmTFESFjduVLsyIlIBmwMi2q9TJ+C994APPwRKSoBTTwXuv58RzERRhs0BER3s4ouVUYTrrgNmz1YuNXzzjdpVEVE7YXNARE1LTgZeegn46ivA5wPOPBO49VagslLtyogoxNgcENGhjRgBbN4M3H03sGiRsgzy00/VroqIQojNARG1LC4OmDsXyM1VRhQuuAC48kpgzx61KyOiEGixORBC3CGESG6PYogozJ16KlBQAEyfDnzwgRKe9OabjGAmijCtGTnoBOB7IcR/hRAjhRAi1EURUXDkFTmQs6QAo+etQ86SAuQVOY78SU0m4OGHgU2blAjma64B/vlPYOfOI39uIgoLLTYHUsqpAI4F8DKA6wD8KoSYJYQ4JsS1EdERyCtyYMayQpQ5PUiPN6PM6cGMZYXBaRAAZdRg/Xrg2WeVlQx9+wLPP88IZqII0Ko5B1JKCaC4/uYFkAzgfSHE7BDWRkRHYHHuDlhMBlhjDNAJAWuMARaTAYtzdwTvRfR64M47lWTFoUOBO+4Ahg0DfvopeK9BRO2uNXMO7hRCFACYDWADgJOklLcBGATg0hDXR0SHaWe5C3FmfcB9cWY9dpa7gv9iPXoAK1YoEcyFhUouwsyZjGAm0qjWjBykAbhESnm+lPI9KWUdAEgp/QBGh7Q6Ijps3VIsqHH7Au6rcfvQLcUSmhdsiGDevh248EJg6lQgMxPIzw/N6xFRyLRmzsE0KeWfzXxve/BLIqJgyM7qAZfHC2etF34p4az1wuXxIjurR2hfuGNH4L//BT76CCgtBQYPBiZPBlwhGLEgopBgzgFRhBrSMxVTR/dBmtWE0mo30qwmTB3dB0N6prZPARddpFxiuOEGYM4c5VLDmjXt89pEdESE1ND65MzMTJnPIUoi7fn6a+Dmm4GiImD8eGW/hsREtauiyMEl9kHGkQMiCr2zzwa2bAHuvVfZr6FPH+CTT9SuioiaweaAiNqHxaJcXsjLA1JTlUmLY8cygpkoDLE5IKL2dcopygqGf/9bmbTYuzfwxhuMYCYKI2wOiKj9mUzKUsdNm4ATTgCys4FRo4A/m1wYRSEQkmhtihhsDoiozYL2wdK7N7BuHfCf/yh/9u0LPPccI5hDLOTR2qR5bA6IIlgofjsM+geLTgdMmKBEMJ9+uvL3M85QwpQoJNolWps0jc0BUYQK1W+HIftg6d4dWL4cWLxY2ZthwABgxgxGMIdAu0ZrkyaxOSCKUKH6EA/pB4sQyhbQhYXAxRcrW0Mzgjno2j1amzSHzQFRhArVh3i7fLB07Ai88w7w8cdAWRkjmINMtWht0gw2B0QRKlQf4u36wXLBBcoowk03KRkJ/foBq1cH/3WijOrR2hT2GJ9MFKEa5hxYTAbEmfWocfvg8niD8iGQV+TA4twd2FnuQrcUC7KzeoT+g2X1aiWC+ffflT9nzwaSkkL7mqQVjE8OMjYHRBFMlQ/xUHK5gEcfBebOBTp1AubPV5IWKdqxOQgyNgdEpD35+cCNNwKbNwOXX67kJHTsqHZVpB42B0HGOQdEpD0NKxhmzACWLlXClBYvZgQzUZCwOSAibTIagSlTlAjm3r2Ba68FRo4EduxQuzIizWNzQETa1hDBPG8e8L//ASeeqFxm8PlafiwRNYnNARFpn04H3HEHsHWrEr08caLyZ2Gh2pURaRKbAyKKHN27A59/rsw/+PlnYOBAZWtoj0ftyog0hc0BEUWWhgjm7duVCOZp05QJjN9/r3ZlRJrB5oCIIlOHDvsjmB0OYMgQ4N57GcFM1ApsDogosjVEMN98sxKedNJJwNdfq10VUVhjc0BEkS8xEXjhBWDNGmXy4ogRyn4NFRVqV0YUltgcEFH0GD5cSVWcPBl47TWgTx8lRImIArA5IKLoEhsLPPEE8O23yryEiy8G/vUvoLhY7cqIwgabAyKKToMGKSsYZs4EPv1UGUV47TVGMBOBzQERRTOjEXjoISWCuU8f4PrrgfPPZwQzRT02B0REJ5wArF0LPP88kJurRDA/+ywjmClqsTkgIgKUVQw5OcC2bcCwYcCkScDppzOCmaISmwMiosa6dQM++wx4803g11+BAQOA6dMZwUxRhc0BEdGBhACuukqJYL7sMuCRR5QJjN99p3ZlRO2CzQERUXPS04G33gI++QTYuxfIygLuuQeoqVG7MqKQYnNARNSSMWOUuQjjxwNPPaVEMK9apXZVhy2vyIGcJQUYPW8dcpYUIK/IoXZJFGZUbQ6EECOFED8LIX4TQjygZi1E1Dx+mECJYF6wQIlgNhiAc84BbrxRGVGAdn5GeUUOzFhWiDKnB+nxZpQ5PZixrDBs6yV1qNYcCCH0AJ4H8A8AfQBcKYToo1Y9RNQ0fpgcYPhw4McfgfvvB15/HejTBz8veF0zP6PFuTtgMRlgjTFAJwSsMQZYTAYszt2hcmUUTtQcOTgVwG9SyiIppQfAOwAuVLEeImoCP0yaEBsLPP64MkGxUyccn3Mdpr42DV3cFWH/M9pZ7kKcWR9wX5xZj53l3Mqa9lOzOegC4K9GX++qv4+I2kFrh8H5YXIIJ58MfPcdXhszHoO2/A8zHroCp637FJAybH9G3VIsqHEHhjvVuH3olmJRqSIKR2o2B6KJ+w4KNRdCjBdC5Ash8ktLS9uhLKLI15ZLBfwwaYHRiO+uGI97p7yO3V2PwQ2v/Bt3PzkBlt1/heXPKDurB1weL5y1XvilhLPWC5fHi+ysHmqXRmFEzeZgF4CjGn3dFYDtwIOklAullJlSysz09PR2K44okjW+VFC5rw47HDUoKqvBpHd+OKhB4IdJy7KzeuD3lC54eOJzeP2a+3F00VbM/ffVmPzTirCLYB7SMxVTR/dBmtWE0mo30qwmTB3dB0N6pqpdGoURIVXagUwIYQDwC4ARAHYD+B7AOCnltuYek5mZKfPz89upQqLINXreOqTHm1G5rw6/llRDrxPQCaDW60fPtLiDPizyihxYnLsDO8td6JZiQXZWD36YHKDxz6i/rMK9S59F8uqVwODBwMsvA337ql1iJGtqJJqOgGrNAQAIIUYBeAaAHsArUsqZhzqezQFRcOQsKUCZ04Mdjhp4vH4YdAJevx8mvR490uKQZjVh/lWD1C5T26QE3n4bmDgRqKwEpkwBHnwQMJnUriwSsTkIMlVzDqSUn0spj5NSHtNSY0BEwdNwqaDa7YVOAF6/Hz4/0CU5Jmwn0mmOEMC4ccrGTf/6F/Doo8oExm+/VbsyohYxIZEoCjVcd06IMaDWq4wY9OoQh5Q4MycbBlt6OrBkCbBsmTKCkJUF3H03I5gprLE5IIpSQ3qm4pmxA9EzLQ490uKQZDFxsmEo/fOfSgTzrbcCTz+tRDB/9ZXaVRE1ic0BURTjzPV2lpAAzJ8PfPMNYDQC554L3HDD3xHMatNKBDSFnqoTEtuKExKJKGLU1gKPPQbMmaNcenjuOeDSS1UrpyH7wmIyIM6sR43bB5fHq5VmkRMSg4wjB0REaoiJAf7v/4DvvwcyMoDLLlOaA7tdlXIYk02NsTkgIlLTwIHKCobHHwc++wzo0wd45RVlKWQ7Ykw2NcbmgIhIbUajssvj5s1Av37KVtDnngsUFbVbCYzJpsbYHBARhYvjjgNWrwYWLFB2fDzxROCpp9olgpkx2dQYmwMionCi0ynLHbdtA84+G7jnHmDoUGDLlpC+LFeuUGNcrUBEFK6kBN55B7jzTiVA6cEHgYceAsxmtSsLN1ytEGQcOSAiCldCAFdeCWzfDlx+OTB9uhLBnJendmUU4dgcEBGFu7Q04M03ldUMTqdymWHSJKC6Wu3KKEKxOSAi0opRo5S5CDk5wLPPKhHMK1eqXRVFIDYHRERaYrUqaYrr1ilzD847D7j+eqC8XO3KKIKwOSAi0qLTTwc2bVImKL7xhhKe9MEHaldFEYLNARFRCIV0M6OYGGDmTCA/H+jSRYlgvuQS1SKYKXKwOSAiCpGGzYzKnB6kx5tR5vRgxrLC4O92OGCAEsH8xBPA8uWqRTBT5GBzQEQUIu26mZHBAEyerEQw9++vSgQzRQ42B0REIaLKZkbHHgt8/TXwwgvKjo/tGMFMkYPNARFRiKi2mZFOB9xyi7Ls8Zxz2i2CmSIHmwMiohBRfTOjrl2Bjz9WIpj/+ENJV3zkEcDtbp/XJ81ic0BEFCJhsZmREMAVVwCFhcDYsUoE88CBQG5u+9VAmsONl4iIosny5colh127gAkTlKWQ8fFqV3WkuPFSkHHkgIgomvzjH8pchNtvB+bNUyYsfvml2lVRmGFzQEQUbaxWpTFYt04JUjr/fOC66xjBTH9jc0DUjkKalkfUVqedpkQwT5kCLFkC9O4NvPcew5OIzQFRe2m3tDyitoiJAWbMUCKYjzoKuPxy4OKLAZtN7cpIRWwOiNpJu6blEbVV//5AXh4wezbwxRdKBPNLL3EUIUqxOSBqJ6qk5RG1hcEA3HefEpY0YABw883AiBHAb7+pXRm1MzYHRO1EtbQ8orbq1UuJYH7xRaCgAOjXD3jyScDrVbsyaidsDojaieppeURtodMB48cr4UnnnquMKGRlKRs7UcRjc0DUTsIiLY+orbp0AZYuBd59F/jzT2DQIODhhxnBHOGYkEhERK3jcAB33w0sXgyccIIyYfG009SuCmBCYtBx5ICIiFonNRV4/XUlgtnlAs44A7jzTqC6Wu3KKMjYHBCFMYYmUVgaORLYulWJYH7uOaBvX2DFCrWroiBic0AUphiaRGGtIYJ5/XrAYlH2bMjOVi49kOaxOSAKUwxNIk0YOhT44Qdg6lTg7beVCOZ332V4ksaxOSAKUwxNIs2IiQH+/W8lE6F7d2DsWOCii4Ddu9WujA4TmwOiMMXQJNKcfv2A3FwlMGnlSiWCeeFCwO9XuzJqIzYHRGGKoUmkSQYDcM89SgTzoEHALbcwglmD2BwQhSmGJpGmHXMMsGoVsGiRMifhpJOAOXMYwawRDEEiIqLQstmUZY9LlyqjCS+/rOwCGTwMQQoyjhwQEVFode4MfPgh8N57wF9/AZmZyuqG2lq1K6NmsDkgIqLQEwK47DJg+3bg6quBmTOVdEUKSwa1CyAioiiSkgK8+ipw5ZXK1tAUltgcEBFR+zvvPLUroEPgZQUiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAJyQSUdTJK3Jgce4O7Cx3oVuKBdlZPZg8SdQIRw6IKKrkFTkwY1khypwepMebUeb0YMayQuQVOdQujShssDkgoqiyOHcHLCYDrDEG6ISANcYAi8mAxbk7VK6MKHywOSCiqLKz3IU4sz7gvjizHjvLXSpVRBR+2BwQUVTplmJBjdsXcF+N24duKRaVKiIKP6o0B0KIOUKIn4QQm4UQHwkhktSog4iiT3ZWD7g8XjhrvfBLCWetFy6PF9lZPdQujShsqDVysBLAiVLKfgB+AfCgSnUQUZQZ0jMVU0f3QZrVhNJqN9KsJkwd3YerFYgaUWUpo5Tyy0Zf5gG4TI06iCg6DemZymaA6BDCYc7BDQCWq10EERERKUI2ciCE+ApApya+NUVK+XH9MVMAeAEsOcTzjAcwHgC6desWgkqJKBow+Iio9YSUUp0XFuJaALcCGCGlbNUaoszMTJmfnx/awogo4jQEH1lMBsSZ9ahx++DyeDnXIHIItQuINGqtVhgJ4H4AF7S2MSAiOlwMPiJqG7XmHDwHwApgpRBikxDiBZXqIKIowOAjorZRa7VCLzVel4iiU7cUC8qcHlhj9v+Tx+AjouaFw2oFIqKQYvARUduwOSCiiMfgI6K2UeWyAhFRe2PwEVHrceSAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKwOaAiIiIArA5ICIiogBsDoiIiCgAmwMiIiIKYFC7ACIiCr68IgcW5+7AznIXuqVYkJ3VA0N6pqpdFmkERw6IiCJMXpEDM5YVoszpQXq8GWVOD2YsK0RekUPt0kgj2BwQEUWYxbk7YDEZYI0xQCcErDEGWEwGLM7doXJlpBVsDoiIIszOchfizPqA++LMeuwsd6lUEWkNmwMiogjTLcWCGrcv4L4atw/dUiwqVURaw+aAiCjCZGf1gMvjhbPWC7+UcNZ64fJ4kZ3VQ+3SSCPYHBARRZghPVMxdXQfpFlNKK12I81qwtTRfbhagVqNSxmJiCLQkJ6pbAbosHHkgIiIiAKwOSAiIqIAbA6IiIgoAJsDIiIiCsDmgIiIiAKwOSAiIqIAbA6IiIgoAJsDIiIiCsDmgIiIiAKwOSAiIqIAbA6IiIgoAJsDIiIiCiCklGrX0GpCiFIAf6pdx2FKA1CmdhFBxPMJb5F2PkDknRPPJ3jKpJQjVXrtiKSp5kDLhBD5UspMtesIFp5PeIu08wEi75x4PhTOeFmBiIiIArA5ICIiogBsDtrPQrULCDKeT3iLtPMBIu+ceD4UtjjngIiIiAJw5ICIiIgCsDkIESHEHCHET0KIzUKIj4QQSc0cN1II8bMQ4jchxAPtXWdrCSH+JYTYJoTwCyGanZEshNghhNgihNgkhMhvzxrbog3no5X3J0UIsVII8Wv9n8nNHBfW709LP2+h+E/99zcLIU5Wo87WasX5nCmEqKx/PzYJIaapUWdrCSFeEULsEUJsbeb7mnp/qHlsDkJnJYATpZT9APwC4MEDDxBC6AE8D+AfAPoAuFII0addq2y9rQAuAbC2FceeJaUcEObLmlo8H429Pw8AWCWlPBbAqvqvmxOW708rf97/AHBs/W08gAXtWmQbtOG/n3X178cAKeX0di2y7V4DcKg8Ac28P3RobA5CREr5pZTSW/9lHoCuTRx2KoDfpJRFUkoPgHcAXNheNbaFlHK7lPJntesIllaej2beHyh1vV7/99cBXKRiLYerNT/vCwEsloo8AElCiIz2LrSVtPTfT6tIKdcCKD/EIVp6f+gQ2By0jxsALG/i/i4A/mr09a76+7RMAvhSCFEghBivdjFHSEvvT0cppR0A6v/s0Mxx4fz+tObnraX3pLW1ZgkhfhRCLBdC9G2f0kJGS+8PHYJB7QK0TAjxFYBOTXxripTy4/pjpgDwAljS1FM0cZ9qy0dacz6tcJqU0iaE6ABgpRDip/rfNtpdEM5HM+9PG54mbN6fJrTm5x1W70kLWlPrRgDdpZTVQohRAJZCGZLXKi29P3QIbA6OgJTynEN9XwhxLYDRAEbIpteM7gJwVKOvuwKwBa/CtmnpfFr5HLb6P/cIIT6CMrSqyodPEM5HM++PEKJECJEhpbTXD+PuaeY5wub9aUJrft5h9Z60oMVapZRVjf7+uRBivhAiTUqp1T0XtPT+0CHwskKIFPSSmAAAAfhJREFUCCFGArgfwAVSSlczh30P4FghxNFCCBOAsQA+aa8ag00IESeEsDb8HcB5UCb+aZWW3p9PAFxb//drARw0MqKB96c1P+9PAGTXz4ofAqCy4XJKGGrxfIQQnYQQov7vp0L5N9nR7pUGj5beHzoUKSVvIbgB+A3KtbdN9bcX6u/vDODzRseNgrKa4Xcow92q197M+VwM5bcCN4ASAF8ceD4AegL4sf62Tevno7H3JxXKKoVf6/9M0eL709TPG8CtAG6t/7uAsgLgdwBbAGSqXfMRns8d9e/Fj1AmLg9Vu+YWzudtAHYAdfX//9yo5feHt+ZvTEgkIiKiALysQERERAHYHBAREVEANgdEREQUgM0BERERBWBzQERERAHYHBAREVEANgdEREQUgM0BkYYIIU4RQmwWQsTUJx5uE0KcqHZdRBRZGIJEpDFCiBkAYgDEAtglpfw/lUsiogjD5oBIY+pz+r8HUAslbtencklEFGF4WYFIe1IAxAOwQhlBICIKKo4cEGmMEOITAO8AOBpAhpTyDpVLIqIIY1C7ACJqPSFENgCvlPItIYQewP+EEGdLKb9WuzYiihwcOSAiIqIAnHNAREREAdgcEBERUQA2B0RERBSAzQEREREFYHNAREREAdgcEBERUQA2B0RERBSAzQEREREF+H9Buv9MuOvHtQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nbi:hide_in\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plot setup\n",
    "plt.figure(figsize=(7,6)) \n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "\n",
    "# true parameters\n",
    "a = -2\n",
    "b = 1\n",
    "sigma = 2\n",
    "numData = 40\n",
    "\n",
    "# generate data\n",
    "x = np.random.rand(numData)*(b-a)+a\n",
    "y = a * x + b\n",
    "eps = np.random.randn(x.shape[0])*sigma\n",
    "y = y + eps\n",
    "\n",
    "# the least squares line\n",
    "meanx = np.mean(x)\n",
    "meany = np.mean(y)\n",
    "varx = np.var(x)\n",
    "vary = np.var(y)\n",
    "cov = np.inner(x-meanx, y-meany)/x.shape[0]\n",
    "\n",
    "# MLE \n",
    "aMLE = cov / varx\n",
    "bMLE = meany - cov*meanx / varx\n",
    "res = y - cov*(x - meanx)/varx - meany\n",
    "varMLE = np.var(res)\n",
    "R2 = 1 -  varMLE / vary\n",
    "\n",
    "\n",
    "# MLE line\n",
    "xvals = np.arange(np.min(x)-0.3, np.max(x)+0.3, 0.1)\n",
    "yvals = cov*(xvals - meanx)/varx + meany\n",
    "\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.figtext(0.9,0.8, r\" $\\hat a$ = {:.4f}\".format(aMLE) +\n",
    "            \"\\n $\\hat b$ = {:.4f}\\n\".format(bMLE) +\n",
    "            r\" $\\hat\\sigma^2$ = {:.4f}\".format(varMLE) + \"\\n\"\n",
    "            r\" $R^2$ = {:.4f}\".format(R2),\n",
    "            ha=\"left\", va=\"top\",\n",
    "            backgroundcolor=(0.1, 0.1, 1, 0.15), fontsize=\"large\")\n",
    "plt.plot(xvals, yvals, c=\"r\", zorder=1)\n",
    "plt.scatter(x, y, zorder=0, alpha=0.7);"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}