5 Estimating distributions and modeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Density Estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objectives\n",
    "+ Definitions and Types of Density Estimation;\n",
    "+ Parametric approaches: MLE and Moment based;\n",
    "+ KDE python examples;\n",
    "+ Modelling: Model Specification and Types of Models;\n",
    "+ Model examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction:\n",
    "In statistics, the probability density function (PDF) can be interpreted as providing a relative likelihood that the value of the random variable would equal the sample [1]. Assume a random variable $X$, and its PDF can be expressed as $f$, which represents the distribution of $X$, therefore, the probabilities between event $a$ and $b$ can be formulated as:\n",
    "\n",
    "$$P(a<x<b)=\\int_a^b{f(d)}\\,{\\rm d}x$$ \n",
    "\n",
    "**Density estimation** is to assume the density or PDF from the set of observed data. In other words, **Density Estimation** is the construction of an estimate of the density function from the observed data [2]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Types of Density Estimation:\n",
    "Generally speaking, the density estimation can be classified into three main categories:\n",
    "\n",
    "***1) Parametric density estimation*** If we have the previous knowledge about the particular form for the density distribution of the observed data set (e.g. Gaussian, Unit, or Gamma), hence we need to estimate the parameters for each specific distribution which could be the best fit to obeserved data.\n",
    "* **Method of Moments** chooses the parameters so that the sample moments (typically the sample mean and variance) match the theoretical moments of our chosen distribution [3].<br>\n",
    "* **Maximum Likelihood** chooses the parameters to maximize the likelihood, which measures how likely it is to observe our given sample [3].<br>\n",
    "* **Bayesian Estiation** minimizes the posterior expected value of a loss function, in other words, it maximizes the posterior expectation of a utility function [4].\n",
    "\n",
    "***2) Non-Parametric density estimation*** Unlike the parametric density estimation, the non-parametric has no-assumption about the density and it will try to estimate the density from the observed data set (without any model). However, the computation for this non-parametric method grows with the data set size increasing.\n",
    "* **Histogram** is a simple non-parametric density estimation method, but highly depends on the structure of the data, and suffers from artifacts, or higher-dimensions.<br>\n",
    "* **Kernel Density Estimation (KDE)** estimate the PDF of a random variable. By using the kernal function, it accumulates or integrates the data, and smooth the density estimation. Buy the proper selection of bandwidth $h$ significantly affects the KDE's performance.<br>\n",
    "* **Nearest Neghbour** <br> kNN (k-th nearest neighbor) using fixed points to determine the volume of region from the data.\n",
    "\n",
    "***3) Semi-Parametric density estimation,*** which is a mixture mode of the previous metioned estimation.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parametric approaches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Method of Moments*\n",
    "The method of moments is one of the oldest methods for point estimation. The basic idea is that the sample moments could give good estimates of the corresponding moments [5].\n",
    "Suppose that the parameter set is denoted as: $\\theta=(\\theta_1,\\theta_2,\\cdots,\\theta_k)$. So there can be $k$ population moments can be determined. For each population moments, or $m$-th, the method of moments can be written as:\n",
    "\n",
    "$$\\mu_k=E[X^k]$$\n",
    "\n",
    "Now the corresponding sample moments, denoted as $m_k$ can be determined:\n",
    "\n",
    "$$m_k=\\frac{1}{n}\\sum_{i=1}^{n} X_i^k$$\n",
    "\n",
    "Note that the number of sample moments should be equal to the number of the parameters, for example:\n",
    "\n",
    "$$m_1=\\overline{X},  m_2=\\frac{1}{n}\\sum_{i=1}^{n}X_i^2$$\n",
    "\n",
    "From the previous equations, the sample moments would have the same estimates as the corresponding population moments, which means:\n",
    "\n",
    "$$\\mu_k=m_k$$\n",
    "\n",
    "Finally, the parameter list $\\theta=(\\theta_1,\\theta_2,\\cdots,\\theta_k)$ can be solved according to the equation, and the solution is the moment estimator of $\\hat{\\theta}$.\n",
    "\n",
    "\n",
    "### *Maximum Likelihood Estimation*\n",
    "Suppose that we have the observed data set (or training data): $X=\\{x_1,x_2,\\cdots,x_n\\}$, which is distributed according to a certain distribution with an unknown parameter $\\theta$, which we need to estimate:\n",
    "\n",
    "$$p_i(X_i=x_i|\\theta)$$\n",
    "\n",
    "$p_i(x_i|\\theta)$ is not a conditional density, is a denote that highlight the density $p(x_i)$ is dependent on the parameters $\\theta$.\n",
    "Parameters $\\theta$ is fixed but unknown. What we will do is finding the parameter value $\\theta$ that maximizes the probability of obtaining the sample that has actually been observed.<br>\n",
    "\n",
    "Again, note that the observed data set $X=\\{x_1,x_2,\\cdots,x_n\\}$ is the $n$ independent samples from our data, (iid). Therefore, we can wirte the likelihood function of $\\theta$ in terms of the observed set of samples $X$ as: \n",
    "\n",
    "$$ p(X|\\theta)=\\prod_{k=1}^n p(x_k|\\theta)$$\n",
    "\n",
    "So maximum likelihood estimation with respect to $\\theta$ will maximize this likelihood function $p(X|\\theta)$, which can be mathematical expressed as: \n",
    "\n",
    "$$\\hat{\\theta}=\\mathop{\\arg\\max}_{\\theta}[p(X|\\theta)]$$\n",
    "\n",
    "Usually, it might be a easier task to maximize $\\log[p(X|\\theta)]$ instead. Because the $\\log$, $\\lg$, or $\\ln$ is monotonic function, hence we have:\n",
    "\n",
    "$$\\hat{\\theta}=\\mathop{\\arg\\max}_{\\theta}[p(X|\\theta)]=\\mathop{\\arg\\max}_{\\theta}[\\log{p(X|\\theta)}]$$\n",
    "\n",
    "The maximum likelihood estimate of value of $\\theta$ can be rewritten as:\n",
    "\n",
    "$$\\hat{\\theta}=\\mathop{\\arg\\max}_{\\theta}[\\log{\\prod_{k=1}^n p(x_k|\\theta)}]=\\mathop{\\arg\\max}_{\\theta}[{\\sum_{k=1}^n \\log{p(x_k|\\theta)}}]$$\n",
    "\n",
    "In this way, the maximization of a long product of terms is transferred into a easier way to maxmize the sum of probabilistic terms. It is a computational-saving strategy.\n",
    "The maximizing the, $Y=\\log{p(x_k|\\theta)}$, term can be achieved by using gradient based algorithm. Assume that $\\theta=(\\theta_1,\\theta_2,\\cdots,\\theta_m)^T$, and the gradient operator is:\n",
    "\n",
    "$$\\nabla_\\theta{Y}=[\\frac{\\partial Y}{\\partial\\theta_1},\\frac{\\partial Y}{\\partial\\theta_2},\\cdots,\\frac{\\partial Y}{\\partial\\theta_m}]^T$$\n",
    "\n",
    "In order to find the optimum, we use the first-order condition which is:\n",
    "\n",
    "$$\\nabla_\\theta{Y}=0$$\n",
    "\n",
    "A set of optimum parameters of $\\theta$ can be derived, in order to check that $\\theta$ satisfies the condition is maximum, the second-order condition might be applied. Also check the upper and lower bound of the range of $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KDE and Python Example\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Kernel Density Estimate*\n",
    "KDE is a non-parametric way to estimate the PDF of an observed dataset. First, it is neccessary to introduce a pre-defined kernel function $K$,which satisfied the codition:\n",
    "\n",
    "$$K=\\int_{-\\infty}^{+\\infty}K(x)\\,{\\rm d}x=1$$\n",
    "\n",
    "In the most cases $K$ will be a symmetric PDF, and normal density, for example, the Guassian Kernel. Therefore, the kernel density estimator $K$ can be written as:\n",
    "\n",
    "$$\\hat{f}=\\frac{1}{nh}\\sum_{i=1}^{n}K(\\frac{x-X_i}{h})$$\n",
    "\n",
    "Note that the parameter $h$ is the window width, or bandwidth. The figure below shows the kernel estimate with individual kernels [copied from 6]. \n",
    "\n",
    "![3.JPG](3.JPG)\n",
    "\n",
    "The shape of the kernel function $K$ outline is determined by the kernel function $K$, and the bandwidth is determined by $h$. More intuitive, the kernel estimator acts as a sum of each kernel window, as well as \"bumps\" in some case.\n",
    "\n",
    "The effect of varying window width is shown in the figure below [copied from 6]. \n",
    "\n",
    "![2.JPG](2.JPG)\n",
    "\n",
    "As we can see, a large vale of window width $h$ will over smooth the density estimate so that the structure of the dataset is hard to identify; a small value of $h$ gives a spiky density estimate which is too rough to explain.\n",
    "\n",
    "So it is important to determine a appropriate value of window width $h$ for density estimation. In order to find a value of $h$, we compute and analyze the minimum squared error between the estimated density and the original density:\n",
    "\n",
    "$$ E[(\\hat{f}-f)^2]=E[\\hat{f}-f]^2+var(\\hat{f})$$\n",
    "\n",
    "Therefore, the $h$ selection became an optimization problem:\n",
    "\n",
    "$$h^*=\\mathop{\\arg\\min}_{h}\\{ E[\\int(\\hat{f}-f)^2\\,{\\rm d}x]\\}$$\n",
    "\n",
    "The first exceptation term is the bias, and the sencond term is the variance. This is the case of bias-variance tradeoff, where the bias could be reduced at the expense of the variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *Finding optimal value of the bandwidth for KDE (Python)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the library for KDE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n",
      "/usr/lib/python3/dist-packages/sklearn/grid_search.py:43: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. This module will be removed in 0.20.\n",
      "  DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()\n",
    "import numpy as np\n",
    "from sklearn.neighbors import KernelDensity\n",
    "from sklearn.grid_search import GridSearchCV\n",
    "from sklearn.cross_validation import LeaveOneOut"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate and initialize the pseudo-random number as the observed dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_data(N, f=0.3, rseed=1):\n",
    "    rand = np.random.RandomState(rseed)\n",
    "    x = rand.randn(N)\n",
    "    x[int(f * N):] += 5\n",
    "    return x\n",
    "\n",
    "x = make_data(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameter optimization. *GridSearchCV* provide an automatically interative hyper-parameter tuning. *GridSearchCV* method traverses multi-combination of hyper-parameters, and determine a set of hyper-parameter with the best performance by cross-validation. <br>\n",
    "\n",
    "In this example, the esimator is a kernel density with guassian kernel; *cv* is the hyper-parameters for cross-validation, or it can also be the train/test generator; *LeaveOneOut* provides train/test indices to split data in train/test sets; *grid.fit()* is for running the grid search; *best_params_* describe the parameters which leads to the optimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'bandwidth': 0.35111917342151316}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bandwidths = 10 ** np.linspace(-1, 1, 100)\n",
    "grid = GridSearchCV(KernelDensity(kernel='gaussian'),\n",
    "                    {'bandwidth': bandwidths},\n",
    "                    cv=LeaveOneOut(len(x)))\n",
    "grid.fit(x[:, None]);\n",
    "grid.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the histogram of the observed dataset, and its KDE results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.02, 0.4)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAH1RJREFUeJzt3WuQXGd95/HvSKOLRxrZI6lHlpBsLrH/YGSyETigVdkgG8iSpaJ4C9auZRccm0SUs6YC8SIXtU6wy4CDjbPRCxO8QAEbHGWdxCsXRlBkfQGEgZXkjCWB/5Y0Y1kajWZac7/2TF/2RffIrdHM9Jme7jndp3+fN+pzebr/6suvn3n6OefUZTIZREQkWhaFXYCIiJSewl1EJIIU7iIiEaRwFxGJIIW7iEgEKdxFRCKoPshOZnY/cBOwHNjp7gem2efLwFZ3f1/QNiIiUh4Fe+5mth24zt23AZ8AHplmn2uAG+bSRkREyifIsMx2YC+Aux8BNphZw5R9HgY+P8c2IiJSJkGGZdYDLXnLcWAd0AZgZrcBzwIng7aZTjKZytTXLw5UtIiInFc33cog4T4+zR1lAMxsNfCfgX8HbAzSZia9vSMBSpleLNZIPD5YdPswqfZwqPZwqPbSi8Uap10fJNw7gOb8+wI6c7dvJNtL/xmwDHiLmf11gTYiIlJmQcbc9wE7AMxsC9Dq7qMA7v6P7v52d38PcDNwyN0/M1sbEREpv4Lh7u4HgRYzOwT8LfBZM7vNzG6eS5tSFSwiIoXVVcopf+PxwaILqdSxsCBUezhUezhUe+nFYo3T/qCqI1RFRCJI4S4iEkEKdxGRCFK4i4hEkMJdRCSCFO4iIhGkcBcRiSCFu4hIBCncRUQiSOEuIhJBCncRkQhSuIuIRJDCXUQkghTuIiIRpHAXEYkghbuISAQp3EVEIijIBbIxs/uBm4DlwE53P5C37U+A24EM8BLwKeC9wBPA0dxuh939rhLWLSIisygY7ma2HbjO3beZ2WbgUeCG3LYG4FbgenefMLNngK25ps+7+0fKVLeIiMwiyLDMdmAvgLsfATbkQh13H3H3G3PB3gA0AmfLVq2IiAQSZFhmPdCStxwH1gFtkyvM7B7gM8BD7t5qZlcA15jZPrKBf5+7/3i2B2lqaqC+fvFc6z8vFmssum3YVHs4VHs4VPvCCBLu41OW68iOr5/n7g+a2W7gaTP7FXAMeADYA1wJPGdmV7t7YqYH6e0dmVPh+Sr1quRBqPZwqPZwqPbSm+kLJ0i4dwDN+fcFdAKY2WrgHe7+nLuPmNkPgPe4+0+Ax3P7t5nZWWADeb19EREpnyBj7vuAHQBmtgVodffRvPbfNLMVueV3A25mt5rZF3Jt1pIdxmkvZeEiIjKzgj13dz9oZi1mdghIAneY2W1Av7s/aWb3Ac+aWZLs2PxTwArgFjPbT/YL4E53nzq8IyIiZRJonru77wJ25a06nLftu8B3pzQZAm6ed3UiIlIUHaEqIhJBCncRkQhSuIuIRJDCXUQkghTuIiIRpHAXEYkghbuISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIIU7iIiEaRwFxGJIIW7iEgEKdxFRCJI4S4iEkGBrsRkZvcDNwHLgZ3ufiBv258AtwMZ4CXgU+6ema2NiIiUV8Geu5ltB65z923AJ4BH8rY1ALcC17v7VuAqYOtsbUREpPyCDMtsB/YCuPsRYEMu1HH3EXe/0d0ncusagbOztRERkfILMiyzHmjJW44D64C2yRVmdg/wGeAhd281s4JtpmpqaqC+fvEcSr9QLNZYdNuwqfZwqPZwqPaFESTcx6cs15EdXz/P3R80s93A02b2qyBtpurtHQlQyvRisUbi8cGi24dJtYejVmu//cFnCu7zrXtuLOq+g6jV572cZvrCCTIs0wE0598X0AlgZqvN7H2QHaIBfgC8Z7Y2IiJSfkHCfR+wA8DMtgCt7j6a1/6bZrYit/xuwAu0ERGRMis4LOPuB82sxcwOAUngDjO7Deh39yfN7D7gWTNLkh1nfyo3FfKCNmX8P4iIyBSB5rm7+y5gV96qw3nbvgt8N0AbERFZIDpCVUQkghTuIiIRpHAXEYkghbuISAQp3EVEIkjhLiISQQp3EZEIUriLiESQwl1EJIIU7iIiEaRwFxGJIIW7iFzgxJl+MplZL78gVUDhLiIXOHyim18c7WQimQq7FJkHhbuIXKSzd4SfvdTB2Hgy7FKkSAp3EZlW//A4Pz9ylsSEevDVKND53EWkNg0Mj/PLX3fybzdfTv3iC/uCYV+PVWannruIzKpnYIyDHtePrFUmUM/dzO4HbgKWAzvd/UDetvcCXwYywHHgj4AbgCeAo7ndDrv7XSWsW0TmaCKZLrptR/cwL5/s5W1vXF3CiqScCoa7mW0HrnP3bWa2GXiUbHhPegy40d3bzex/A78PDAHPu/tHylG0iMzdybOD82rvp/q4rHEZ69esKFFFUk5BhmW2A3sB3P0IsMHMGvK2v9vd23O3zwGrSluiiMxXOpOh9Uz/vO/n0CtxhscmSlCRlFuQYZn1QEvechxYB7QBuHsfgJmtB94P3AtcC1xjZvuARuA+d//xbA/S1NRAff3iOf8HJsVijUW3DZtqD0ct1X6qc5C6gJ+vFSuWzbrd2wd4/+9eGei+pquzlp73MAUJ9/Epy3Vkx9fPM7Nm4PvAp92928yOAQ8Ae4ArgefM7Gp3T8z0IL29I3MqPF8s1kg8Pr8/OcOi2sNRa7W/+JuzDA/P+PG7QKH9hocT7D/4WqD7mlpnrT3vC2GmL5wg4d4BNOffF9A5uWBmq4AfAve6+w8BcsM0j+d2aTOzs8AGcr19EVk4o4kkXb2jJb3PY6fnP8Qj5RVkzH0fsAPAzLYAre6e/075KrDb3Z+eXGFmt5rZF3K315IdxmlHRBbc6fhQyacxpjUtsuIV7Lm7+0EzazGzQ0ASuMPMbgP6gR8BHweuyq2DbI/9ceAWM9tP9gvkTnefOrwjIiUU5KAiqR2B5rm7+y5gV96qw3m3Z/r15eZiixIRkfnREaoiIhGkcBcRiSCFu4hIBOmskCJygaf2vxp2CVIC6rmLiESQwl1EJIIU7iIiEaRwFxGJIIW7iEgEKdxFRCJI4S4iEkEKdxGRCFK4i4hEkMJdRCSCFO4iIhGkcBeRskmMp8IuoWYp3EWkbM50D4ddQs0KdFZIM7sfuAlYDux09wN5294LfBnIAMeBP3L39GxtRKQ2nDk3zJvWrwq7jJpUsOduZtuB69x9G/AJ4JEpuzwGfDS3/RLg9wO0EZEacK5/TEMzIQkyLLMd2Avg7keADWbWkLf93e7enrt9DlgVoI2IlFA6kwm7hGllMhnO9oyEXUZNCjIssx5oyVuOA+uANgB37wMws/XA+4F7yYb7jG2m09TUQH394rnUfoFYrLHotmFT7eGIUu3d/aMhVTK7FSuWMTSeuqDeKD3vlSxIuI9PWa4jO75+npk1A98HPu3u3WZWsM1Uvb3Ff7vHYo3E44NFtw+Tag9H1Gp/5VRfSNXMbng4QevoOFdvWMWS+kWRe94rwUxfOEGGZTqA5vz7AjonF8xsFfBD4C/c/YdB2ohIacX7KrPnDpBKZ+iq4PqiKki47wN2AJjZFqDV3fNfqa8Cu9396Tm0EZESSaXT9Awmwi5jVme7Ne6+0AoOy7j7QTNrMbNDQBK4w8xuA/qBHwEfB67KrQN43N0fm9qmLNWLCL2DCVKpdNhlTCvIxba/dc+N5S+kBgWa5+7uu4BdeasO591eFrCNiJTBub6xsEuQCqQjVEWq3Ll+hbtcTOEuUsWSqTQ9gwp3uZjCXaSK9Q4mSKcr8wAmCZfCXaSKaUhGZqJwF6li5yr0yFQJn8JdpEolU2l6K3x+u4Qn0FRIEQnX7Q8+E3YJUmXUcxcRiSCFu4hIBCncRUQiSOEuIhJBCncRkQhSuIuIRJDCXUQkghTuIiIRpHAXEYkghbuISAQFOv2Amd0P3AQsB3a6+4G8bcuBx4Br3P1duXXvA54AjuZ2O+zud5WwbhERmUXBcDez7cB17r7NzDYDjwI35O3yEPAicM2Ups+7+0dKVqmIiAQWZFhmO7AXwN2PABvMrCFv++eBJ8tQm4iIFCnIsMx6oCVvOQ6sA9oA3H3QzNZM0+4aM9sHNAL3ufuPZ3uQpqYG6usXB6t6GrFYY9Ftw6baw1HNtUfJ2rUrqaurC7uMQKrpPRMk3MenLNcBha7rdQx4ANgDXAk8Z2ZXu/uMJ5/u7R0JUMr0YrFG4vHBotuHSbWHo5prj5oTr3Zz6cplYZdRUKW+Z2b6wgkS7h1Ac/59AZ2zNXD3duDx3GKbmZ0FNpDr7YuITIr3j1VFuFebIGPu+4AdAGa2BWh191mv7WVmt5rZF3K315IdxmmfX6kiEkXxPl0qsBwK9tzd/aCZtZjZISAJ3GFmtwH97v6kmT0BbALMzJ4jOy3yKeAWM9tP9gvkTnefOrwjIkJ3/xjpdIZFi6pj3L1aBJrn7u67gF15qw7nbfvoDM1unkddIlIjJq8Fu+bS5WGXEik6QlVEQqehmdJTuItI6OL9CvdSU7iLVLhMptDM4+rXO5ggmUqHXUakBBpzF5nq9gefKbjPt+65cQEqib7B0YmwSyi7dDpDd/8Y61Y3FN5ZAlHPXaTCxXtrY8iiS+PuJaWeu1wkSK9cFk5XjYS7flQtLfXcRSpYKp3mXI382DgwPM5oIhl2GZGhcBepYN39Y6TS0f9BdZJ676WjcBepYJ01MiQzSeFeOgp3kQrW2VP82VKrUVffaE1M/VwICncpm9PxobBLqGpDoxMM1cA0yHyJ8RT9wzoNVSko3KVsDrzcxZG2bvXEilRrvfZJtTI7qNwU7lJWx0/3c+iVc6QV8HN2tkbDvXMeF+6R1yncpexOdQ3ScuycevBzMJFM0d0/FnYZoegdSDCR1KkI5kvhLgviZOcg/lpf2GVUjc6e0Zr9ayedyWjWTAko3GXBvPxaL+36kTWQWh2SmaShmflTuMuCevHYuZqbATJXqXS65sOtq1dTIucr0LllzOx+4CZgObDT3Q/kbVtO9tJ617j7u4K0kcpV7rHOZCrNgZe7uOG3N+iyajPo6h2t+THn0USSgeFxXTh7Hgr23M1sO3Cdu28DPgE8MmWXh4AX59hGKlAmk6Hl+LmyP07fUIKXX+st++NUqzPnarvXPqnWh6bmK0jPfTuwF8Ddj5jZBjNrcPfJZ/7zwBrgY3NoIxXoVNdQSQ88emr/qwX30TnfL5RKpznbMxx2GRWhs3cUu6Ip7DKqVpBwXw+05C3HgXVAG4C7D5rZmrm0mU5TUwP19YuD1DytWKyx6LZhq4Tah0bGae3sYMWKhf0zePWalSwOaXimEp73qU53DbJ02RKW1tBoxEzvufFUhpWrLuGSZZVzZvJKfM/MJMizNvVY4Dqg0C8dc27TO48fkGKxRuLxwaLbh6kSak9nMux/qYO+gYWfV/2rlnZ+a+OlC/64lfC8T+ewdzE8nAi7jAU12//36CtdXHl5ZQRqpb5nZvrCCTJbpgNozr8voLMMbSQkJ9r76Q4h2CE7PVLn8M6aSKbp0DjzBTo0RFW0ID33fcAXga+Z2Rag1d0LHWFQTBsJwcDwOL85Gd6Pm8lUml+/2sM7rbnwzhHX0T1MqgYvEl3ot5nHnvq1fpspQsFwd/eDZtZiZoeAJHCHmd0G9Lv7k2b2BLAJMDN7DnjM3R+f2qZ8/wUpVjqd4dArcdIhXwziVNcQb1q/itWrlodaR9he69QBXlI6gX6pcPddwK68VYfztn00YBupMMdO99E3VBnju0faerj+Heupq6vNue9DoxM1czk9WRg6QrVG9Q0lKupcLz0DY5zprt3x5pOdlfdDXSVJpWtvuGq+KmeOkSyYVDqdHY6pgMO788dbZxp7jfp4azqd4ZSGZGYV7xvj8tUNYZdRVRTuNejl1/oY0NVuKsYnv/Js2CVUvPb4sMJ9jhTuNeT2B58JuwSRopztGSaVXsPiRRpJDkrPlEiIekI6vqDaTCTTuvzeHCncRUJ0vL0/7BKqRntcBzTNhcK9RnR064NRaQZHxumo4RlCc9XRM0KyBg/yKpbCvQaMJpL867Hyn8pX5uaVU/26IMUcpFJpzpxTJyUohXvEZTIZXjwWJzGRCrsUyTM4Ml7S0yvXCj1nwSncI+54e79+iKpAL7/Wp157EeJ9YzrRXEAK9wjrHUyEelIwmV7vYEIXCi9SJpPhVJeeuyA0zz0iojyHvXcwQVNjNK5ekclkONzaHXYZVe1U1xBXbby0Zs9DFJR67lLxjrR1R2YI41TXkOa2z9PgyDg9A5VxwrtKpnCXitfdPxaJKYOJiRRHX+0Ju4xI0InWCtOwjFSFo209rFt9SVUffn6ktYfEuGYtFUMXW5+76v2kSE0ZHpvg+OnqPZqzo3uYU13qbcrCCdRzN7P7gZuA5cBOdz+Qt20r8NXctn929wfM7H3AE8DR3G6H3f2uUhYuteeV0/1sal5Jw/IlYZcyJ6OJJP96XAeRycIqGO5mth24zt23mdlm4FHghrxdvkM2+NuBF8zs73Prn3f3j5S6YKldqVSalhPdvOeadVUzUyKdyV7KUMMxstCCDMtsB/YCuPsRYIOZNQCY2ZuBHnc/5e5p4PvAB8tVrEyvlg7q6OwZ4XQVnUDKT/YS79NBZLLwgoT7eiCetxwH1s2wrQu4PHf7GjPbZ2Y/M7MPzLtSmVYm1zOsJS+dOFcVX2jt54bxU5VzKUOpLUHG3KdesqcOyBTYdgx4ANgDXAk8Z2ZXu/uMk1Obmhqor18cqOjpxGKNRbcN23xqP9razchEtM+Ut2LFxQcwvXJmkBvftYlFi4ofninne6a7f5R7v/HLst2/XGwhMqCaciZIuHcAzXnLMaBzhm2XA2fcvR14PLeuzczOAhuAtpkepLe3+HnMsVgj8Xh1zkSYT+29gwl+0XKmIq6FWk7Dwxf3CdqGE/ysLsPb3ri6qPss53tmaHSCn750piz3LTMrdwZUas7M9IUTJNz3AV8EvmZmW4BWdx8FcPfTZrbEzK4g+4Pqh4GbzexW4K3u/gUzW0t2GKe9BP8PyZlIpjnwclfkg302fqqPVSuW8obYygV/7Cif7kGioWC4u/tBM2sxs0NAErjDzG4D+t39SeDPyP7gmgH+zt1Pmdn3gVvMbD/Zcf073V1XZC6RTCZDy/FzDI9NhF3KgpjtAJb8bTqIpbaNT6RYuqT4od2oCTTP3d13AbvyVh3O2/YT4Hem7D8E3FyKAuVip7qGdF5rkSmOt/dzTZHDdFGkI1SrzMDIOC0ndFbB6XT2VP/5Z6R4rWcGdDxBHp1bpgpofDeYX/y6k7decRlXbbqMRVVykJOUTjKV5pXTfVz75jVhl1IR1HOXyMhkMvzmZC/7X+pgYFg/8dSito4BhkZr47eoQhTuEjndA2M8+2I7L74SL8sHfXBEXxyVKp3OcLRNp1UGDctIRGUyGU52DnKyc5DmpkvYGFvJ+jUNLJnHgXKJ8RTHTvfR2jFQwkql1Dq6h+nsHWFdU0PYpYRK4V7hRmpkumMpBDnndxD5UypT6TTd/WOcjg/THh8ila7d4wqqyUvHu9m+ZTn1i2t3cELhXsHGJ1K8cLSz8I5SUi+d6GZ8IsXw2AQDw+MK9Co0PDbByyd72VzDP64q3CtUMpXml7/u1PhuCFrPVO9FQeR1J84McPnqBtZedknYpYSidv9mqWDJVJpf/aaTbl1IWaRomUyGgzV8Ln2Fe4WZDPauXp0DXGS+RhNJDnhtnoNJ4V5BEuMpfn7krIJdpITifaMcPtFNpsYCXmPuFaJ/KMEvf9Ol2TEVoFSzbmRhBX3dauUEcwr3kGUyGdo6BjjS2q1ZGSILIJPJVM01eOdD4R6iodEJbr/7qbDLEKkpLce7ecdb1szrKl7VQOEeAh3pKBKeV88OMDAyzrssRsPyJWGXUzYK9wU0ODJOW8cAJzuHSKWifd1TkUrWMzDGM4faeduVTbxp/apI9uIV7mWkU/WKVK5kKs3h1m5aOwa4euNlbGxeweJF0ZlAGCjczex+4CZgObDT3Q/kbdsKfDW37Z/d/YFCbaIsMZ6iZ3BMByCJVInh0QlePBbn6Ks9vGHtCtavaWDNpcurPugLhruZbQeuc/dtZrYZeBS4IW+X75AN8XbgBTP7e+CKAm2qmnrkItUryJTJP7z+TVy2chmrViyl8ZIlrGxYQizWWP7iSihIz3072Qtg4+5HzGyDmTW4+4iZvRnocfdTALkLY38QWD9Tm/L8N0onk8mQTKVJTKRJjKcYm0gxNp5kLJFiNJFkJJEMu0QRKbP/89O2BXuscs27DxLu64GWvOU4sA5oy22L523rAjYUaDOtpqYG6os813Ymk2Ht2pV5y7l/cwuZ3LpMJkM6nSGd+zeVzgZ5MpVhIpliIplmfCJ1vj31i1lev5jllyy94PEeuitWVJ0iUp3qFy/itzZdFnYZcxIk3KeelrCOXG7Osm22NtPq7S2+Ux+LNRKPDxbdfhGwrA6WLVkESxZ2nG2+tYdJtYdDtYenEmufabgoSLh3AM359wV0zrDtcuAMkJyljYiIlFmQcN8HfBH4mpltAVrdfRTA3U+b2RIzu4LsD6ofBm4mG+zTtomCr3zlS3zuc58PvM9M++/Y8SF+/vOfnd8HuKjNjh0fAmDv3n0XrQM4deo1br31Y3z964+yefO1ALzwwn62bt12vs2ePd+jvz97jvITJ06zceNaPv3pz/Lwww+yaNEiVqxYyebN13LkyGF27ryTPXu+R1dXJ+Pj2T/Ali5dyvj4OEuXvj481dy8jsHBAfr7+89vX7RoEel0+vzy1q3bOHLkMIODA9x99z08/PCDbN26jRde2E9j4yoSiTHe+c7rzu8z2X7r1m1s23b9+cearGdy382br2XbtuvPPxdTb+/Z8z1uvfVjFz3fn/vc5/nKV77E7t2PkEgkLnjOt2x5O5s2XXH+Oct/LWZ6DbZtu579+3/K3r372LHjQ+dfi0m7dz/CsmXLOXHi9AWvGXD+/7d79yM0N6/j9OlTbNy4CYD+/v7z/8c9e77H6dOnWJQ3cyOTycx6EqzGxlUXPJ+l0Ni4CoDBwbkfeDf1/ZFvct3dd9/D7t2PnH+NE4mxC95DO3feyde//igAmzdfy8GD/+/8vpPP7969+wDYsuXt51+H6T6Dk/sG+RxPNVObYu6rFG1nUxfkTGlm9lfAB8j2yO8A3gn0u/uTZnYD8Ddkh13+zt0fma6Nux+e7THi8cGiT6yy0H/qNTevoqtr9jd5/j4z7d/cvIpMJkM8Pkhzc/bDM7VN/vqp6wqZbt+5tK900z0/QffPP7/ITPeT/1rkb5v6WIUeO0rPeSXKf12Ai16byXXTfbYKfY7zxWKN1NXVzfhZnst9laptrq5pj8AKNM/d3XcBu/JWHc7b9hPgdwK0ERGRBVLds/RFRGRaCncRkQhSuIuIRJDCXUQkgnRWyCLcffc9c9pnpv23bt024z6Ty9Ptk79utqmQk23yp0JCdmpaFKZCTj4XU2/PNBVy8vnYvfuRi57zjRs3sWnTFdO+FjO9BpNTISeXZ5oKOfU+QFMhJ5ViKiRc+Pxu3Ljpotc//zXN/2zM1UxtirmvUrSdTaCpkCIiUl00LCMiEkEKdxGRCFK4i4hEkMJdRCSCFO4iIhGkcBcRiSCFu4hIBEXqICYzWwe8DNzs7s+FXE4gZrYY+J/AVcAy4L+5+/PhVlWYmd1P9sLoy4Gd7n4g5JICM7Mvkb028BLgr9z9iZBLCszMLgGOAve7+7dDLmdOzOw/AX9O9sps97r70yGXFIiZrQT+F9BE9v1+n7vvC7eqwqLWc38IaA27iDn6GDDm7tcDtwFfDbecwsxsO3Cdu28DPgE8EnJJgeWuP/Bv3H0r2Yu5/3XIJc3Vfwe6wy5irnIB+efANrIX9fnDcCuak9sAd/f3AR8B/keo1QQUmXA3sxuBAfLONV8l9gB3526fA6rhqg7bgb0A7n4E2GBmDeGWFNjPgf+Yu90HLDWzqvgcmNlbgbcBVdHjneL3gKfdfczdz7j7H4dd0Byc4/XLhjYB8RBrCSwSwzJmthS4l2xv4G9CLmdO3H2c1y8o/mfA4yGWE9R6oCVvOQ6sA9rCKSc4d08CQ7nFTwI/cPfSnICl/B4G/ivZnmS12QTEzGwfsBL4S3d/JuSagvoH4A4zc2A18Ach1xNI1YW7mX2S7Icy3z7ga+7eb2YhVBXMDLX/pbv/yMz+FHgX2T9ZK934lOU6spdZrBpmtgP4Y7KXgqx4ZvZx4Cfu/molv8dnsYxsB+DDwJuB/2tmb6ySL9b/Apx09w+Y2W8D3yT7Wa1oVRfu7v4N4Bv568xsP/AhM/ss8Bbgd83so+5+NIwaZzJd7QBmdgfZvzr+INeTr3QdvP5nKkAM6Aypljkzs98D/gL4oLv3hV1PQP8eeJOZ/QdgI5Aws9Pu/i8h1xXUWeAFd08Bx8xsAFgLdIVbViBbyXYgcfcWM3uDmdXn/gqsWFUX7tPJ/bAHgJl9G/h2pQX7TMzszcCfAte7+2jY9QS0D/gi8DUz2wK0VkvtZnYp2R+Ab3T3qvlh0t1vmbxtZl8AXq2iYAf4F+AbZvYw2c5AI9mx7GpwgmxP/Z/M7A3AYKUHO0Qk3KvcJ4HLgKfz/tz+YCX34N39oJm1mNkhIAncEXZNc3AL2R/F/iHv+f64u78WXknR5+7tZvYk8CzZYL+rSoZkAP4W+I6ZPQ8sBT4Vcj2B6HzuIiIRVBVTwEREZG4U7iIiEaRwFxGJIIW7iEgEKdxFRCJI4S4iEkEKdxGRCPr/y+9W3Jz89egAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x6b0e1a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist = plt.hist(x, bins=30, normed=True)\n",
    "from sklearn.neighbors import KernelDensity\n",
    "x_d = np.linspace(-4, 8, 1000)\n",
    "# instantiate and fit the KDE model\n",
    "kde = KernelDensity(bandwidth=0.35, kernel='gaussian')\n",
    "kde.fit(x[:, None])\n",
    "\n",
    "# score_samples returns the log of the probability density\n",
    "logprob = kde.score_samples(x_d[:, None])\n",
    "\n",
    "plt.fill_between(x_d, np.exp(logprob), alpha=0.5)\n",
    "plt.plot(x, np.full_like(x, -0.01), '|k', markeredgewidth=1)\n",
    "plt.ylim(-0.02, 0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modelling: Model Specification and Types of Models\n",
    "Model specification is to choose the specific type of model to apply to the data [6]. It can be classified into two types of models:<br>\n",
    "* **Descriptive Models** presents the main feature of the data, and serves as a description characterization of data curve-fitting.<br>\n",
    "* **Explanatory Models** is used to explain the phenomenon of the data and represents the complex functions with respect to the combinations of simplified components.<br>\n",
    "\n",
    "The figure below shows the family of the model specification [copied from 7].<br> \n",
    "\n",
    "![1.JPG](1.JPG)\n",
    "\n",
    "In machine learning, we applied a variety of models to data, and the model could learn to perform multi-tasks, such as prediction, classification. In terms of the purposes, models can be divided into (1) the models that can be used to perform supervised learning, and (2) that for unsupervised learning.<br>\n",
    "\n",
    "Supervised learning aims mapping between the input variables and the output variables, while unsupervised learning is trying to find the features or characteristics on the given dataset.<br>\n",
    "\n",
    "### Linear Models\n",
    "Linear models can be expressed as the sum of the multiplications between weights and inputs. Suppose that the linear equation $y=ax+b$, so we can regards the parameters $a$ and $b$ as the weights and vecterized as: \n",
    "\n",
    "$$w=[w_1,w_2]=[a,b]$$\n",
    "\n",
    "Also the inputs can be rewritten as: \n",
    "\n",
    "$$x=[x1,x2]=[x,1]$$\n",
    "\n",
    "The model can be rewritten as:\n",
    "\n",
    "$$y=w_1x_1+w_2x_2=wx^T=\\sum_{i=1}^{n}{w_i}{x_i}$$\n",
    "\n",
    "Linear models are readily available, and easy to fit.\n",
    "\n",
    "### Non-linear Models\n",
    "Nonlinear models are models in which the prediction cannot be simply expressed as a weighted sum. In general, the nonlinear models consists of three types: Linearized models, Parametric non-linear models, and Non-parametric nonlinear models.\n",
    "\n",
    "#### Linearized Models\n",
    "Linearized models are similar to the linear models, but the input space differs from that of linear model. Theinput space has been transformed by using nonliear functions, such as Gaussians, Polynomials, Logistics, or Sinusoids.\n",
    "\n",
    "Suppose that a non-linear model $y=ax^2+bx+c$, where $a$, $b$, and $c$ are all parameters. Obviously, if $a\\neq0$, the outputs or predictions are not linear with respect to $x$. However, if we present the weights and input set as:\n",
    "\n",
    "$$ w=[w_1,w_2,w_3]=[a,b,c]$$\n",
    "\n",
    "$$ x=[x_1,x_2,x_3]=[x^2,x,1]$$\n",
    "\n",
    "then the model can be rewritten as:\n",
    "\n",
    "$$y=w_1x_1+w_2x_2+w_3x_3=wx^T=\\sum_{i=1}^{n}{w_i}{x_i}$$\n",
    "\n",
    "In this case, the model can be viewed as a linear model with respect to the the new input space $[x^2,x,1]$.\n",
    "\n",
    "The figure below shows the Linear and Linearized models [copied from 6].\n",
    "\n",
    "![4.JPG](4.JPG)\n",
    "\n",
    "#### Parametric Nonlinear Models\n",
    "Parametric nonlinear cannot be rewritten or transfered into any form of linearized model, but it can be expressed by the certain format of mathematical operations. In other words, the parametric nonliner models can not be linear with respect to the free parameters in the models.\n",
    "Consider the model $y=x^n$, where $n$ is a parameter, but the output of this model cannot be linear in terms of $n$, because the term $x^n$ has no way to be expressed as $ax^2+bx+c$, where $a$ $b$,and $c$ are weights.\n",
    "\n",
    "#### Nonparametric Nonlinear Models\n",
    "Nonparametric nonlinear models that can be generically applied to arbitrary datasets. In order to distinguish between parametric and nonparametric nonlinear model, there are some generalities concluded:<br>\n",
    "* Nonparametric nonlinear models are nonlinear in terms of the free parameters, so these cannot be categorized as linearized models.\n",
    "* Unlike nonparametric nonlinear models, parametric nonlinear models make more assumptions on the form of the nonlinearity of the model.<br>\n",
    "\n",
    "Nonparametric nonlinear models are powerful, but the trained or fitted models can be less intuitive, and these models are suffering from overfitting and local minima. With the data size increased, the computational complexity of nonparametric nonlinear models will increase as well.\n",
    "Examples: Nearest-Neighbor, Local Regression, and Neural Network.\n",
    "\n",
    "\n",
    "The figure below shows the Parametric and Nonparametric Nonlinear Models [copied from 6].\n",
    "\n",
    "![5.JPG](5.JPG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model examples\n",
    "This example illustrates the construction of a linear regression model for load cell data that relates a known load applied to a load cell to the deflection of the cell. The model is then used to calibrate future cell readings associated with loads of unknown magnitude [8]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data collected in the calibration experiment consisted of a known load, applied to the load cell, and the corresponding deflection of the cell from its nominal position. Forty measurements were made over a range of loads from 150,000 to 3,000,000 units. The data were collected in two sets in order of increasing load. The systematic run order makes it difficult to determine whether or not there was any drift in the load cell or measuring equipment over time. Assuming there is no drift, however, the experiment should provide a good description of the relationship between the load applied to the cell and its response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array([0.11019, 0.21956, 0.32949, 0.43899, 0.54803, 0.65694,\n",
    "               0.76562, 0.87487, 0.98292, 1.09146, 1.20001, 1.30822,\n",
    "               1.41599, 1.52399, 1.63194, 1.73947, 1.84646, 1.95392,\n",
    "               2.06128, 2.16844, 0.11052, 0.22018, 0.32939, 0.43886,\n",
    "               0.54798, 0.65739, 0.76596, 0.87474, 0.98300, 1.09150,\n",
    "               1.20004, 1.30818, 1.41613, 1.52408, 1.63159, 1.73965,\n",
    "               1.84696, 1.95445, 2.06177, 2.16829])\n",
    "x = np.array([150000, 300000, 450000, 600000, 750000, 900000, 1050000,\n",
    "         1200000, 1350000, 1500000, 1650000, 1800000, 1950000,\n",
    "         2100000, 2250000, 2400000, 2550000, 2700000, 2850000,\n",
    "         3000000, 150000, 300000, 450000, 600000, 750000, 900000,\n",
    "         1050000, 1200000, 1350000, 1500000, 1650000, 1800000,\n",
    "         1950000, 2100000, 2250000, 2400000, 2550000, 2700000,\n",
    "         2850000, 3000000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in analyzing the data is to select a candidate model. In the case of a measurement system like this one, a fairly simple function should describe the relationship between the load and the response of the load cell. One of the hallmarks of an effective measurement system is a straightforward link between the instrumental response and the property being quantified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEGCAYAAABlxeIAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFHdJREFUeJzt3XuQnXV9x/F3zKUQNsVYlw1UKSW031oYJyZQRBQCBIKgXKqSKmotl7HYOkNtC9pLVIp2rCKdArYKyEXxAiLgEIFJsAYERUiIA1S/lojK9ZBK0KxBk5DtH+dZ58yye87Z5Dzntu/XTIbzXPY53x9P2A/P73l+v2fayMgIkiS9qNMFSJK6g4EgSQIMBElSwUCQJAEGgiSpMKPTBTRrw4ZNffc41Ny5s9m4cXOny2ipfmwT9Ge7bFPv2Jl2DQ7Omdbsvl4hdNCMGdM7XULL9WOboD/bZZt6R7vaZSBIkgADQZJUMBAkSYCBIEkqGAiSJMBAkCQVDARJ6mKVSoXLL7+cSqVS+nf1zMA0SZpqKpUKCxfuz9atW5g5cxZr1z7E0NBQad/nFYIkdalVq25j69YtAGzduoVVq24r9fsMBEnqUkuWLGXmzFkAzJw5iyVLlpb6fXYZSVKXGhoaYu3ah7jnntUcfPDhpXYXgYEgSV1taGiI008/nQ0bNpX+XXYZSZIAA0GSVDAQJEmAgSBJKhgIklSi4eFh1qy5l+Hh4U6X0pBPGUlSSYaHhzn66MNYv/5h5s/fj5Ur72BgYKDTZU3IKwRJKsm6dWtZv/5hANavf5h169Z2uKL6DARJEmAgSFJpFixYyPz5+wEwf/5+LFiwsMMV1VfqPYSI+ChwBDAT+FhmXlez7RDgAmAX4KuZeX6ZtUhSuw0MDLBy5R1kfp+IV3T1/QMo8QohIg4DFmTmIcAxwIVjdrkKWAYcCLwxIuaXVYskdcrAwACLFh3U9WEA5XYZ3Q2cUnx+FpgVES8CiIh9gWcy89HM3A7cTDU0JEkdUlqXUWZuA0YfvD0D+Hrxyx9gT2BDze5PA3vVO97cubOZMWN6y+vstMHBOZ0uoeX6sU3Qn+2yTb2jHe0qfRxCRJwInAkcXbN6y5jdpgEj9Y6zcePmFlfWeYODc9oyg2E79WOboD/bZZt6x860azJBUvZN5aXAcuCYzHy2ZtOTwB41y/OAJ8qsRZJUX5k3lXcHPgkcl5k/q92WmY8BMyNi74iYDrwBuKWsWiRpR1UqFa655uq2vOS+08q8QlgGzAW+HBGj674BPJCZNwBnAzdR7Sr6fGY+WmItkjRp7X7JfaeVeVP5M8Bn6my/A3hVWd8vSTtrvJfcn3rqOztcVXkcqSxJE2j3S+47zdlOJWkCoy+5X7XqNpYsWdrX3UVgIEhSXUNDQ33dTVTLLiNJEmAgSJIKBoIkCTAQJEkFA0FSX5tKI413lk8ZSepbU22k8c7yCkFS3xpvpLEmZiBI6ltTbaTxzrLLSFLfmmojjXeWgSCpr02lkcY7yy4jSRJgIEiSCgaCJAkwECRJBQNBUldzpHH7+JSRpK7lSOP28gpBUtdypHF7GQiSupYjjdvLLiNJXcuRxu1lIEjqao40bh+7jCRJgIEgSSoYCJIkwECQJBUMBEkSYCBIKplTT/QOHzuVVBqnnugtXiFIKo1TT/QWA0FSaZx6orfYZSSpNE490VsMBEmlcuqJ3mGXkSQJMBAkSQUDQZIEGAiSpIKBIKkuRxpPHaU+ZRQRBwA3ARdm5sVjtv0YeBR4vlh1amY+XmY9kibHkcZTS2mBEBG7ARcBt9fZ7fWZOVxWDZJ2zngjjX2EtH+V2WX0a+A44IkSv0NSiRxpPLVMGxkZKfULIuJDwP9N0GX0beDlwF3A+zNzwmK2bXt+ZMaM6eUVKmlcTz31FCtWrOD4449n3rx5nS5Hkzet2R07OVJ5ObASeBq4HngLcO1EO2/cuLlNZbXP4OAcNmzY1OkyWqof2wT92a5m2zR9+m6ccMIpAF3/76AfzxPsXLsGB+c0vW/HAiEzrx79HBG3Aft3qhZJUoceO42IORGxOiJmF6teCzzYiVokSVVlPmW0CLgA2AfYGhFvBr4GPJKZN0TEtcC3ImIzcD/wlbJqkSQ1VlogZOYaYHGd7ZcAl5T1/ZKkyXGkstTnhoeHWbPmXoaHHfKj+nwfgtTHhoeHOfrow1i//mHmz9+PlSvvYGBgoNNlqUt5hSD1sXXr1rJ+/cMArF//MOvWre1wRepmBoIkCTAQpL62YMFC5s/fD4D58/djwYKFHa5I3cx7CFIfGxgYYOXKO8j8PhGv8P6B6jIQpD43MDDAokUHdboM9QC7jCRJgIEgSSoYCJIkoMl7CBFxFPDXwEuomVs7Mw8rqS5JUps1e1P5U8D5wE9LrEXSOCqVCqtW3caSJUt9n7FK1Wwg/CQzP1dqJZJewJfcq52aDYQVEfGXwH8DW0dXZuaPSqlKEuBL7tVezQbC3wLPA+fWrBsB9m15RZJ+Y/Ql96NXCL7kXmVqKhAyc++yC5H0QkNDQ6xd+5D3ENQWzT5ltBfwEeBAYDvwHeCfMnNDibVJohoKdhOpHZodh3AZsBZYBrwN+CFweVlFSZLar9l7CLtm5kU1yw9FxBvLKEiS1BnNXiHsGhG/O7pQfN6lnJIkSZ3Q7BXCvwD3RcQTVEcq7wGcUVpVkqS2a+oKITNXUH3E9LTizx9k5q1lFib1i0qlwjXXXE2lUul0KVJdda8QIuIvMvOKiDhvnG1k5vLySpN6nyON1UsaXSFsL/75/Dh/tpVYl9QXxhtpLHWrulcImXlV8fHnmfnvtdsi4oOlVSX1CUcaq5c06jI6AjgSeHtEvKRm00zgXcCHyytN6n2ONFYvafSU0Q+APYvPz9es3wq8tZSKpD7jSGP1ikZdRk8CX4iIu4FZmflDgIhYmJlr21GgJKk9mh2Ydhbw/prlcyPiYyXUI0nqkGYD4bDMPG10ITOXAa8rpyRJUic0GwgjETFrdCEidgOml1OSJKkTmp264lLgfyLifqohciDgY6eS1EeanbriCmAx8Hngc8ChmXlleWVJ3aNSqXD55Zc79YT6XlOBEBG7AG+hei/hRuD3ImKg1MqkLjA69cQZZ5zBwoX7Gwrqa83eQ7iU6niEQ4vlA4ArSqlI6iJOPaGppNlAeFlmngM8B5CZnwbmlVaV1CVGp54AnHpCfa/Zm8qjwTECv3nKyBfkqO+NTj1xzz2rOfjgw516Qn2t2SuEGyNiJbBvRPwHcD/VG8xS3xsaGuL00083DNT3mrpCyMwLI+JOqvcQRoC3ZuaaUiuTJLVVo9lOjxyz6oHin7tHxJGZ+Y0GP38AcBNwYWZePGbbIcAFVLuevpqZ50+qcklSSzW6QvjnOttGgAkDobjPcBFw+wS7XAUcBTwOfDsivpiZ6xvUI0kqSaNAuCozr4yIMzPz0kke+9fAccC5YzdExL7AM5n5aLF8M3AM8J+T/A5JUos0CoQPFnMY/U1EPD92Y2Z+dqIfzMxtwLaIGG/znsCGmuWngb3qFTJ37mxmzOi/6ZMGB+d0uoSW67Y2PfXUU6xYsYLjjz+eefN2/GnpbmtXK9im3tGOdjUKhL8DXg+8mPFnN50wEBrYMmZ5GsUjrRPZuHHzDn5V9xocnMOGDZs6XUZLdVubWvWS+25rVyvYpt6xM+2aTJA0ekHO9cD1EfGm4nOrPAnsUbM8D3iihceXgPFHGvv2Mml8jZ4y+kBm/itwUkScOHZ7Zu7Qf1mZ+VhEzIyIvaneVH4DcPKOHEuqx5fcS81r1GU0+prMu4EBql09G6l28dQVEYuoPla6D7A1It4MfA14JDNvAM6m+kjqCPD50RvMUiv5knupeY0CYXVE3AgsBNYB24FXA7cBp9f7wWLg2uI62+8AXjWZYqUd4UvupeY0mrrig8BjwL6ZeUJmngTsDQwDHyq5NklSGzUKhGOAs4tHSAHIzC1Uu3uOLbMwSVJ7NQqEzbVhMCozt1K9SpAk9YlGgbC1zrbtrSxEktRZjW4qvyYifjrO+mnAS0uoR3qBSqXiU0JSGzQKhHHnnZDapVUjjSU11mik8k/aVYg0HkcaS+3T7BvTpI7wncZS+zT7TmWpIxxpLLWPgaCu50hjqT3sMpIkAQaCJKlgIEiSAANBklQwECRJgIGgNqhUKlxzzdVUKpVOlyKpDh87VamcekLqHV4hqFTjTT0hqTsZCCqVU09IvcMuI5XKqSek3mEgqHROPSH1BruMJEmAgSBJKhgIkiTAQJAkFQwENeRIY2lq8Ckj1eVIY2nq8ApBdTnSWJo6DATV5Uhjaeqwy0h1OdJYmjoMBDXkSGNparDLSJIEGAiSpIKBIEkCDARJUsFAkCQBBkLfc9oJSc3ysdM+5rQTkibDK4Q+5rQTkiaj1CuEiDgPOArYBXh3Zt5Xs+3HwKPA88WqUzPz8TLrmWpGp50YvUJw2glJ9ZQWCBFxBHBQZh4aEQcAnwIOG7Pb6zNzuKwapjqnnZA0GWV2GR0B3ASQmQ8Ce0XE7BK/T+MYnXbCMJDUSJldRnsC36tZ3gAMAY/UrLs0Il4O3AW8PzNHJjrY3LmzmTFjeimFdtLg4JxOl9By/dgm6M922abe0Y52lRkIW8YsTwNqf+EvB1YCTwPXA28Brp3oYBs3bm51fR03ODiHDRs2dbqMlurHNkF/tss29Y6daddkgqTMQHgS2KNmeRD4zcPwmXn16OeIuA3Yv8RaJEkNlHkP4RbgRICIWAj8KDOfK5bnRMTqmnsKrwUeLLEWSVIDpQVCZq4BvhcRa4H/At4XEe+KiJMzcxPV7qFvRcS3gGeAr5RVSy9zpLGkdil1HEJmngucW7PqgZptlwCXlPn9vc6RxpLayZHKXcyRxpLayUDoYr7gXlI7ObldF3OksaR2MhC6nC+4l9QudhlJkgADQZJUMBAkSYCBIEkqGAiSJMBAKJ1TT0jqFT52WiKnnpDUS7xCKJFTT0jqJQZCiZx6QlIvscuoRE49IamXGAglc+oJSb3CLiNJEmAgSJIKBoIkCTAQJEkFA6EBRxpLmip8yqgORxpLmkq8QqjDkcaSphIDoQ5HGkuaSuwyqsORxpKmEgOhAUcaS5oq7DKSJAEGgiSpYCBIkgADQZJUMBAkScAUCASnnpCk5vT1Y6dOPSFJzevrKwSnnpCk5vV1IDj1hCQ1r6+7jJx6QpKa19eBAE49IUnN6usuI0lS8wwESRJgIEiSCgaCJAkwECRJBQNBkgTAtJGRkU7XIEnqAl4hSJIAA0GSVDAQJEmAgSBJKhgIkiTAQJAkFQwESRIwBaa/LktELAauAx4qVj0A/AtwNfBi4DHg1Mz8dUScDJwD7AJclJmfjYjpwKeAA4Bpxb6PRMQfApcCs4H7gPdk5khEnAW8o1j/D5n59Ra35wDgJuDCzLw4IvZoZ1siYjfgCuBlwC+BZZn5TAntuhJYBPys2OXjmbmil9oVER8FjgBmAh8DVtPj52qcNh1PD5+niJgNXAkMAbsB5wHfocvPk1cIO2d1Zi4u/rwX+DhwRWa+GvgxcGpEzAE+ARwLHAqcExEDwDuB7Zl5KPBR4MPFMS8FzsnMg4BB4IiImA+8GzgcWAp8IiKmtaoRxV+ci4Dba1a3uy3nAGsy8zVUf4GfXVK7AD5Qc95W9FK7IuIwYEFmHgIcA1xIj5+rCdoEPXyegBOA+zLzcOBNRd1df54MhNZaDHyt+HwT1ZNzEHBvZv48MzcDdwGvo/p/QzcV+94KLI6IWcB+mXnPmGMcDtyamVszswI8CfxRC+v+NXAc8EQH21J7jNF9y2jXeHqpXXcDpxSfnwVmAUfS2+dqvDaN97upZ9qUmV/KzH8rFl9G9YpgMV1+nuwy2jl/HBG3AHOoJviczHyu2PY0MA/YE9hQ8zMvWJ+Z24pLxEHgmXH2/cUEx/h+KxqRmduAbRFRu7rdbak99ui6MtoF8N6I+HvgKeCveqldRZuGi8UzgK8DJ/TyuZqgTdDD52lURNxTHOs44M5uP09eIey4/wXOp9rX+Q7gMqp9faOmASPAljE/N9F6gG1N7ju6vky139mOtoxdX1b7Pgf8Y3Epfy/Vvt2ea1dEnAicSbUboC/O1Zg29cV5ysyDgZOBL42pqSvPk4GwgzLz8cz8QmZuz8xHqP5fzEBxMwmqafwE1cu3PWp+9AXri0vBrVTT/MX19h2zvkyb2tyW2vV7UlL7MvP2zLyvWFwB7N9MTd3UrohYCiwHjs3MZ+mDczW2Tb1+niLiwIjYGyAz11L9XfvLbj9PBsIOiog/i4gPFZ9fSvVpgsuAE4td/pTqX+TvAq+MiN2Lm0UHA3cCtwAnFfu+AViZmduBdRFxSLH+5OIYK4GlETEzIvYCXpKZPyy5ibe2uS21xxj9vpaLiGsj4pXF4muBB3upXRGxO/BJ4LjMHH0Cp6fP1Xht6vXzBLyG4iZuRAxR7Va+mS4/T95D2HE3A8si4i6qwfoe4H7gixHxPiCBLxf9f8uBO4DtwHmZ+VxE3AicEBH3AZuBtxXHPRe4IiJmAN/MzLsAIuKzVC+dt9OCJ3BqRcQi4AJgH2BrRLwZOBW4po1t+XTx7+4+qv2dy0pq13LgsojYDGwCTsvMLT3UrmXAXODLNfdG/hy4qofP1Xht6vXz9Oniu+8Efovq74c1tPf3w6Tb5PsQJEmAXUaSpIKBIEkCDARJUsFAkCQBBoIkqWAgSOOIiH0i4rEWH3NJRHyzlceUWslAkCQBDkyTJiUiTgPOAp6jOtjn9Mz8eUScSXVitl8Vf5Zl5rMRcRLwEarTBmSHypaa4hWC1KRibprzgKMy8zDgp8D7is27AicWk7H9BHh7sf5i4E2ZeTTwfJtLlibFQJCat5Dq3PW/KJZvBw4sPv8KuC4iVlOdd/6lEfE7wK6Z+YOa/aWuZSBIzRv7lrppwPaI+H2qb7U6pbhCuGGC/f3vTV3Nv6BS8+4FDoyI3y6Wj6X6nty5wIbMfLKY+fZoqhOa/YzqpHqjM7a14i1wUmmc3E4aR0TsQ/Um8LdrVn+X6lvqzqL6es5Hqb7Q5Tngq1TnnF8PXEf13sEyqq9PXA78CPgB8CfF/Qep6xgIkiTALiNJUsFAkCQBBoIkqWAgSJIAA0GSVDAQJEmAgSBJKvw/TdB/Lz3GZjkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x72071470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, y, \"k.\")\n",
    "plt.xlabel(\"Load\")\n",
    "plt.ylabel(\"Diflection\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using software for computing least squares parameter estimates, the straight-line model:\n",
    "\n",
    "$$ D=\\beta_0+\\beta_1L+\\epsilon$$\n",
    "\n",
    "is easily fit to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Least squares estimate values: [  7.22102581e-07   6.14968421e-03]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEVCAYAAADgh5I1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFXWx/FvQsIeIEAHUUBF8KiDKCIiLqi4gMy4j+M2uOEuvoIKiDLuihI0qLiigvvu6DjugqLI4iAygMp1ZFVAuoGEsEOW94+qxO6QkISk053O7/M8Pnbdqq46tzvU6VtV996kwsJCREREkmMdgIiIxAclBBERAZQQRETEp4QgIiKAEoKIiPiUEEREBICUWAcgicPMCoH2zrnfonycy51z48tYdw4wHGgENADmAtc451ZEM6bqYGZ7AYsBV8rq451zy3fy3uLPxMwWAMc451ZVQ0xtgJ7OuX9VdV8S/5QQpFYxs3pAJrBDQjCzA4BHgMOcc0vNLAXIAp4D+kU7LudcfjXsKt85t18lj70bMAz/M6ns+8txHHACoIRQByghSNSZWX1gDNAXyAdeBe5xzhWa2WHAY0AzIA+43jn3uX8yfwLojfd3+l/gYuBdoLn/K/hk59zisEPtD4Scc0sBnHN5ZjYMaOjH0co/dmdgPpANLHPOjSzZuglfNrOL8VodKcBKYICfcC4GTgNaAN8DN5jZ5cAQvMux3wFXOefWm9kxeMmpkR/rbc65Nyv5OTYFXgT2w2v9TAKuAaYB7fzPpCuwFWgPGHAv8A1wJrAGuBovoe4HTHDOjfD3fTNwKVAP+An4O9ARGAekmFlT59y5ZnYqcA9QH1gIDHTO/V6Zekj80j0EqQnX4Z2cugA9gDOA/v66p4GxzjkDHgCe9Mv74p2Q9gM64Z3Ae+GdtPKdc/uVSAYAU4E9zOx9MzvTzFo65zY757L99cPxEsbewPV4J/O8nQVuZq2Bp4C+zrnOwC/AP8I2ORG40jl3g5l1B+4A+vi/0tcDt/nbjQGGOOf2B/7ifwaVdRGQ4+9jX2A78Ce8z2SZ/5lsC9s+H+gGvI/3WebjJdk/A8cAQ82siZkdDNyE9910xkugg5xzs/ESwlt+MtgDmAic7ddvKvDoLtRD4pQSgtSE04BnnXPbnXMbgRf444TYA3jNf/013okL4He8X/xnAI2dc7c55z7Z2UH8a+aHAMuBh4GQmX1uZl39TXoDb/nbLvKPt1POudVAc+fcslJiBPjZOfdzWD3fC/vF/GRYPVcCF5rZ/s65hc6588s4ZD0zW1Div9fC9tHLzE4C6jnnrnHOzSmnCjnOuS+dc4XAj8AU59xG55wDCoEMfx/tnHPrnHMFeC2KjqXsqy8ww38veMnlNDNLKicGqSV0yUhqQgYwxszu8pcbADP91wOAQWbWBO9yRRKAc+47M7sKGAw8b2bvAYPKO5DfargKwMz2B24GPjSzDkBLvMtERdaUtz//0tW9ZtYX799LM+DnsE3WlqjnX82sj7+c7NcVvMtdtwOfmdlmYLhz7p1SDlnmPQTn3Dtm1gK4C9jfzF4AbiynCuvD913Kcj0zSwPGmlkvP+ZWwAel7CsDONK/NFVknb/96nLikFpACUFqwu/AQyWfVPGfqnkS7ymW782sE/C/ovX+9v/yT4LP4J38ni3rIGZ2ELCl6Besc+4nMxsE5ALpQA7QPOwtbfCe6gEowE9G/gmyyPl4l7d6O+dCZnYZ3vX1sur5gnPuhpIrnHNr8S5TXe8njHfN7BO/xVRhzrnngOf8G8lv4SXUhZXZRymGAPvg3YzfYGb3AO1K2e53YJJz7swqHk/ilC4ZSU34F3C5maWaWZKZjfR/cTcHNgA/mVkqfgvAzJqa2aVm9g8A51wOXqJIwrtunlzipF3kZOAV/2SJfynjPOBH59waYDrezVXMrDNwRNh7l+PdrwD4G97lFPASyWI/GWQAFwBNd1LPM/3tMLPTzGyYX+8pZtbW324u3r2Lggp8dsXM7DYzu9RfXAUs44/PpKnfmtkV6Xif0QY/KZ/CH3XcjnfTHOBT4Gh/G8zsMDPL2sVjShxSC0Gq25dmFn6j9jK8G5N7AvPw/uZmAmOBjcCHeJdgVgP/h3dP4VO8G68Tzex/eCdPB1yCd4liKrDYzE51zk0LO9Zof/9f+MmgAPgB7wQHcD/wmpkt9svDL4vcAjxuZkvwbsLm4F3CehU438wWAovwHu9838xGUaK/gHNutv/r+gv/5BzCewpnu5k9A0zy48rDu2m7uZTPr16JSzJFRuA9YTTBzIbjJaxpePdjGuBduvrNf2qrsp4C3vY/67l4LZm3zexavO/iBjOb7pzrZWYD/XUN8L6/63bheBKnkjQfgtRV/kn6N+fcHbGORSQe6JKRiIgASggiIuLTJSMREQHUQhAREV+tecooFFqfcE2Z9PTGZGdvinUY1SoR6wSJWS/VqfaoSr0CgbQK9yRXCyGGUlLqxTqEapeIdYLErJfqVHvUVL2UEEREBFBCEBERnxKCiIgASggiIuJTQhAREUAJQUREfEoIIiJxKi8vj8zMUZxwwglkZo4iL2+nM75WWa3pmCYiUpfk5eVx8Z/24ZPsbK4HMidN4vnnn2XOnAWkpETn1K0WQhUtX/4bw4YN4bLLLuSyyy5kzJhRbN684zD3M2ZM45//fKvM/bz44kTmz59bqWPPnj2LkSOH7VB+zDE9GTToCq666lIuv/wiJk36dKf72bhxA99+O6NSxxaR6Aq0D/BJtjfj68N+WTAYJCsrM2rHrFMthIKCAiZP/oy5c/9L164H0afPiSQn73pOLCgoYOTIYQwaNITu3XsA8MorLzJq1F3cddeoiO0OP/yIsnYDwIABF+9yHCU1bdqUceOeBiA7ey3Dhg2hVavWHHzwIaVu79wCvv12Bocddni1xSAiu6b+Zx/T/IK/RZSFT9E3c+b0qB27ziSEgoICLr74fD7++MPisn79+jNx4iu7nBS+/XYGHTrsWZwMAM477++cd95ZrFmzmiefHEdKSgo5OTkcffQxLFq0kEGDBjN27Bjmz5/LAQfsh3P/47bb7mbixGc49tjjadCgAe+++xZJScksXryI448/kUsuuZxZs75l/PgnSE1NJS0tjbvuur9CMaant+Taa6/nlVde5OCDD+G1117iiy8mUVBQQK9eR3LppVfw0EOj2bRpI+3bd6BLl66MGXMfKSmpJCcnc/fd99OsWfPyDyQiVVNQQGC3FhFFNwAl5yjt2bNX1EKoM5eMJk/+LCIZAHz88YdMnvzZLu9z2bIldOq0b0RZUlISHTvuw9KlSwBo3rwFo0aNKV6/ePEi5s37L08/PZHLLruMBQt+jEhIycnJ/PjjD9x66x089dRzvP32GwBs2LCeW265nXHjnqZp07RK/Uro3Hlfli715pIvLIRHHnmCp5+eyEcffcDGjRs4//wB9OlzIqeddibZ2Wu59tohPProU3TtejCffvrRrn48IlJBDV96fodksHLpKl7OyIgoy8jIYMiQoVGLo860EObO/W+p5fPmzeWEE/ru0j7z8vLJz8/fobywsICkJG+Awf33PyBi3ZIlizjggC4kJyezzz77kJHRZof3m+1Hw4YNi/cFkJbWjAcfvJ+CggJWrFhOt27dady4cYXiTEpKpqDAGyw2JSWFIUMGkZycTE5ONrm5uRHbNm/enCeeeJRt27axenWIE0/sV6FjiMgu2LKFQIfIk37u+IlsPe1MUoA5cxaQlZXJ99//h27dejBkyNCo3VCGOpQQunY9qNTyAw/susv73GefTrz7buSN4oKCAhYvXsRee3UEICUlNWJ9YSEkhQ1Gm5S048i09ertOLLhfffdSWbmw3TsuA9jxozaYf3OzJ8/l06dOrN8+W+89dZrTJjwMo0bN+GCC/66w7YPP/wgF1xwIb16HcVLL01k69atlTqWiFRM49H30WRM5KXf0Kp1ESeIlJQUhg4dQSCQRii0Puox1ZlLRn36nEi/fv0jyvr160+fPifu8j67d+/BihXLmTZtanHZa6+9zAEHdCE9Pb3U97Rr146ff3YUFhbyyy+/EAyuqtCxNm/eTNu2u7NuXQ6zZ89i+/btFXrfunU5PPnko5x99rmsX7+e9PSWNG7chPnz5xEMBtm+fTtJSUnFLZ3c3HXssUd7tm7dyrRpU6P+3LNIXZOUk00go1lEMsh590NCwdzIX4sxUGdaCMnJyUyc+AqTJ3/GvHlzOfDArlV+yiglJYXRo8cydmwmEyaMp7CwkP33/xPDh99a5nv23Xc/2rfvwOWXX0T37t3Ye++OpbYISjrrrL9xzTUDad9+Ty66aCATJoxn4MArS912w4YNDBp0BYWFhWzdupULLriIQw45lPz8fJo0aco111xGly4HcuaZZzN2bCaDBg3mySfH0aZNG84++zxuvXUou+++B+ed93fGjh3D8cefRKdOnXf5cxIRT9p1V9Hw9VeKl/Pb7MbaeT/HMKJItWZO5USZMW3btm1MmvQpJ5/8F5o0qUffvv144433onpdsCbVVNO2piVivVSn6MvLyyMrK5NFX07mrf/MjFi3dsoM8kvcYyxLVepVmRnTEuMsVIvUr18f537irbdeJzW1HpdffnXCJAMR+UNeXh4HH7wfbweDHBVWvvWIo8h998My3xdLOhPFwODB3mNj8fZrRkSqz2sjbmJVMBhR1g44/8ijid6Do1WjhCAiUs1a7bcXQ9auLV5+Brjcfx3NnsZVpYQgIlJNUr+eQouzTokoaw6E9/aJZk/jqqozj52KiERNYSGBjGYRyWD98Ftpk5ERkQyi3dO4qqLaQjCz+4DjgFTgAefcm2HregEPAg2Bd5xz90QzFhGRaGjw9hs0u/qyiLLQryFo0IA5199IVlYmM2dOp2fPXlHvaVxVUYvMzHoDBzvneplZS2Au8GbYJs8DxwPLgelm9qpzbmG04omGlStXcOGF52K2X3FZ585G//5/4auvvmTgwCuZOnUKPXseQWpqKl988TnHHXdChfb99tuvk5OTE9HXYPbsWbzzzhvcc8/oiG1vv30Et9xyOw0aNKyeiolI+bZtI9CudURR7iNPsPXcC4qXi3oa1xbRTFXTgKIxXHOA+maW7JwrMLOOwFrn3K8AZvZv4CTgiSjGExUdOuxZPNR0uM6dDfB6Lh9ySA9SU1N56aXnK5wQKuPOOys3lIWIVE2jR8fS9O7bIspCv+dAFTq6xoOoJQTnXB6wwV+8DPjQOVfgL7cFQmGbB4Hdd7a/9PTGpKSU36O3Jm3d2oSUlGQCgbSI8pkzZ/Lyyy/Tp08ffvxxPiNGDOHYY49l4cL/ceedIxg3bhxZWVl899135OXlceGFF9K/f3+mT5/OfffdRyAQoHXr1rRv3z5i3y1aNKZBg9QdjtenTx/ef/99xo8fz6ZNm1i0aBHLli1j5MiR9O7dm08//ZQJEyaQnJxM165dGT58eFQ/l5LxJYpErJfqVEm5udC8xHDwH38MffsSiN5RgZr5rqJ+McvMTsN74ip80KBtJTZLAnbaEzk7e9NOj9PkjpE0eP/dXQmxTFtPOZ2Nd5R9a2Pt2o3k5RXs0JcgJ2cTW7du58gjj6dly7GMGpVF48aNeeaZZ7j99lF8/vnXLF68lJdeeonfflvNwIED6Nr1MB54YDQjR95Fx46dGDr0elq2zIjYd9F+Sx4vP7+A1as3sGVLHkuWLGPUqCymT/+GF154mb322o9x4x7nySefIzU1lZEjh/PFF9/QpcuuD+q3M4natyIR66U6la+op/GMGdO4f/16jpgzu3hdYeMmrF6y0luI8udYxZ7KFd422jeV+wK3ASc553LCVq0Ewsd83Q1YEc1YomXZsqUMGnRF8XKPHj058MDSR1YtsmDBD8yfP48BAwawbVsehYUFrFmzmt9//52OHTsBcPDBh7Bly5ZKx3PQQd0AaNOmDRs2rOfXX5fy++8rGDLkWsCbLnPlyhVRSwgiiaKop3FyMMjKEuuyP/+KvK4HxySuaIrmTeXmwENAH+fcmvB1zrnfzCzVzDrg3VT+C3BGVY638Y57dvprPlpKu4cwe/asct/Xv/8p3HTT4DKzftH8BZUVPlBe0ThVnTsbY8c+vkv7E6mrsrIymRgMcnJY2SzgjaEjGJqAyQCi2w/hHCAdeN3MvvT/u83Mik78g4H3gP8ArxTdYE40SUnJxUNLF53kDzigC9OmTSU/P5+tW7cydqw3aXYgkMGyZUsA+P7776rl+B067MWyZUvJzvZ6TT777FOsXh0q510idVu9//3M6MxREclgb6AH8d3TuKqieVP5aWDHx2/+WP8V0C1ax48X3bodwqBBV/Dww0+w777GlVdewlNPTeCQQw7l3HPPZfv2PE4/3Zuo5oorruHWW4fRtu3utGnThtJGop0zZ3bEJaqRI+/a6fEbNmzI9dffyE03XU9KSgpm+9GqVeudvkekLmvZoyv1/ClwAV4DzgtbH889jatKw1/HkG7q1R6JWC/VKVLKzBmkn3JSRNn+rVuzYPXq4uWMjAzmzFlQ453LNPy1iEhNKCwk0CbyUdJN19/IxltvZ7L/lFFt6WlcVYlbMxGRctT/4H2aX3JBRFlo6Spo1AiofT2Nq0oJQUTqnrw8Aru3jCha/8BDbLnksjLeUDcoIYhIndLw2adIGxE54mhoZTZUYG7zRKeEICIJK7yn8dHde3Dv2DER69c9/yrbTv5zjKKLP0oIIpKQinoaB4NB7geGfz0lYn1o1TpIqvADOHVC7R6aT0SkDFlZmRQEgxQC4cM5jjt/AKFgrpJBKZQQRCQhnTxhfMSQygvwRtF887eEHBShWighiEhCSV6ymEBGM44L61C2L7C//zqRexpXlRKCiCSOgw+m1WF/jDb8af36JAH/85fjfU7jWNNNZRGp9VLmzCb9pGMjylbP/4UDW7ZkaB3qaVxV+mREpFZr3a41Sdv+mHNr0xVXs/GeBwDvBFeXehpXlRKCiNRK9T//hObnnx1ZuH49Gzcn3DiYNUb3EESkdikoIJDRLCIZbLjzPu9R0qZNYxhY7acWgojUGg1ffoG0IYMiykLL10BqaowiSixKCCISt4qGnpg9bSqff/N1xLrc8RPZetqZMYosMSkhiEhcKhp64qpgkM9LrNOwE9GhhCAicenJUXexKhiMKDsG6Dl0BEOVDKJCCUFE4k7adVdx5+uvFC+vAPYoWkjgSe5jTQlBROJG8m+/0uqQP0WUdQF+CFvW0BPRo8dORSQutPjLSRHJYOsRR9EmIyMiGWjoiehSC0FEYqreD/NpedwREWVr5vxEwe57MKeOTXIfa/pkRSRmWu2/N8lr1hQvbx5wMRsefKR4ua5Nch9rSggiUuNSv55Ci7NOiShb/cuvFDZrHqOIBJQQRKQmFRYSaBN50t844h9s0n2BuKCEICJREz7J/VXNmnP+B/+KWB/6NQQNGsQoOilJCUFEoqKop3F2MMi2EutyH3mCredeEJO4pGx67FREoiIrK5OLSkkGw2+6WckgTqmFICLVLml9LqMzR0WU9QU+BXp/OyMmMUn51EIQkWrVdPgNtN6nXfFyLpCElwxAPY3jmVoIIlItklf9TqsD940oOz49ncnZ2cXL6mkc39RCEJEqa37umRHJYHu3QwgFc3nph4UMHTqC3r2PZejQEcyZs0A9jeOYvhkR2WX1fna0PKpHRNmaWfMo6LAnoJ7GtY1aCCKyS1oeemBEMthy5l8JBXOLk4HUPmohiEilpMyYTvqpfSPKVi9YTGHLVjGKSKqLWggiUqa8vDwyM0dx1lmnkDn6PgIZzSKSwcbBNxEK5ioZJAi1EESkVEU9jYPBIGcAo7+eErE+tHQVNGoUm+AkKtRCEJFSZWVlsiYYpBB4J6z87RP7EQrmKhkkoKi2EMysC/AekOWcG1di3RLgVyDfL7rAObc8mvGISMXt/s6b5JUoqwcctXULvWMRkERd1BKCmTUBHgUm7WSzk51zG6IVg4jsgo0bCezdlsFhRacC7/uv1dM4cUXzktFWoD+wIorHEJFq1OSOkQT2bhtRlsQfyUA9jRNbUmFhYVQPYGZ3AKvLuGQ0HWgPfAPc7JwrM5i8vPzClJR60QtUpC4LhSAjI7Js+nTyDj2Ue++9l6lTp3LUUUdx6623qqdx7ZNU4Q1jmBAuBD4DgsDbwCvOuTfK2k8otD66gcZAIJBGKLQ+1mFUq0SsEyRmvYrq1OySv9MgbOKavM77kv3NrBhGtusS8XuCqtUrEEircEKIWap3zr1Q9NrMPgH+FKtYROqkhQsJdOoUUbR2+nfk79M5RgFJrMXksVMzSzOzKWbW2C86Cpgfi1hE6qL0Y3pBWDLY2u/PhIK5SgZ1XDSfMuoOPAjsBWw3s78C/wIWO+f+aWZvAFPNbBPwPfBWtGIREU/K7Fmk9+sTUbZ6/i8Ulrx/IHVS1BKCc+474NidrH8MeCxaxxeRsEnup3/DF998Tb3we4aDBxO65a7YBSdxR48LiCSooqEnDg0G+arEutCiFd7jpQl4A1Z2nRKCSIIa++ADrAoGI8oGA/WHjmBo06axCUrimhKCSAJq+NLzPPDgAxFlqUAe0Hvm9JjEJPFPCUEkkWzZQqBD5A3ivwFvhi1r6Akpi0Y7FUkQje+/Z4dk0CYQiEgGGnpCdkYtBJFaLil7La1tr4iynPc+YnuvI5njP2U0c+Z0evbsxZAhQzX0hJRJfxkitVjatVfQ8M3Xipfzd9+DtXN+Kl7WJPdSGUoIIrVQ8q/LaNW9S0TZ2q9mkr/f/jGKSBKB7iGI1DIt+p8QkQy2HX2MN+yEkoFUUYVaCGZ2PDAIaEnYUKrOOU2cJBJFxT2NZ0zjjI77MOT55yLWr/nvAgra7h6j6CTRVPSS0ePAPcCyKMYiImHCJ7lfC6SHTXK/ecDFbHjwkdgFJwmpoglhqXPuxahGIiIRsrIyOTAY5PMS5bf93xCuG3lnTGKSxFbRhPCBmV0FfAFsLyp0zi2KSlQidV1BAaMzR0UU3QKMAnrP+Z7rYhKUJLqKJoQbgXxgeFhZIdCx2iMSqeMavPkaza69IrIM2Oa/Vk9jiZYKJQTnXIdoByJS523bRqBd64iiQWlpPLb+jxFJ1dNYoqmiTxntDtwLHAoUADOAkc65UBRjE6kzGj3yEE3vuSOiLPR7DrcWFNBYPY2lhlT0L+sZ4CMgE++x037As8CpUYpLpE5IWp9L633aRZTlvP5Pth93PAApycnqaSw1pqIJoZFz7tGw5R/M7JRoBCRSVzQdOoRGzz9bvFyQ1ow1C3+LYURS11W0p3IjM9ujaMF/3TA6IYkktuTfVxLIaBaRDLInfa1kIDFX0RbC3cAsM1uBd8koA7gsalGJJIjwnsaHH34Ed8+YRoOwDmbbux1Czidfxi5AkTAVfcroAzPrCOyLlxCcc25zVCMTqeXCexrvD4wOSwQAa2bNo6DDnrEJTqQUO00IZnaJc26Cmd1Vyjqcc7dFLzSR2i0rK5NgMMgSIPy0//3+B9BuyowYRSVStvJaCAX+//NLWVdYzbGIJJRNn3y4wz+SVkCXQAZvxSIgkXLsNCE45573X65zzo0NX2dmt0ctKpHarLCQQJvmjAsruhcY6b9WT2OJV+VdMjoO6AP83cxahq1KBS4GNMKWSJj6779H84EDIsoaAVv81+ppLPGsvEtGC4C2/uvwy0bbgfOiEpFIbZSXR2D3lhFF6zPHsuGCC7lOPY2llijvktFK4BUzmwbUd879DGBmhzjnZtdEgCLxruEzT5J2y7CIstDKbKhXjxRQT2OpNSraMe1q4Oaw5eFm9kAU4hGpPTZuJJDRLCIZrHvhNULBXKhXL4aBieyaiiaE3s65S4sWnHPnAEdHJySR+NfkjpEE9m4bURZatY5t/frHKCKRqqvoxcxCM6vvnNsGYGZNAP0EkjonafVqWh8QOQ1I9oefk3foYTGKSKT6VDQhjAd+NLPv8VoVhwJ67FQSXtHQE7Nnf0vWsl858H8//7Gu875kfzMrhtGJVK+KDl0xwcw+A7rjDV1xvXNOI3FJQisaeiItGOSXEuvWTv+O/H06xyQukWip0D0EM2sInI13L+FdYE8zaxrVyERiLCsrk89LJIN3gWFDRygZSEKq6E3l8Xj9EY70l7sAE6ISkUgcSJk9i9GZozgwrKwNcAYwc+b0GEUlEl0VTQjtnHPDgM0AzrmngN2iFpVIrBQW0rptOun9+hQXZeFdJw36yxp6QhJVRRNC0XaFUPyUkSbIkYRS/7OPCbRpTlL+H53y927dmhvCttHQE5LIKvqU0bv+TeWOZvYI3pzKj0UvLJEalJ9PoG16RNGGu0ex+cprmeY/ZfT99/+hW7ceGnpCElpFnzLKMrOv8e4hFALnOee+i2pkIjWg4UvPk3bDdRFloeVrIDUVgJSUFIYOHUEgkEYotD4WIYrUmPJGO+1Tomie///mZtbHOTe5nPd3Ad4Dspxz40qs6wU8iHfp6R3n3D2VilykKrZsIdAhI6Iod/xEtp52ZowCEom98loI/9jJukKgzITg32d4FJhUxibPA8cDy4HpZvaqc25hOfGIVFnj+++hyUOjI8pCq9ZBUlKMIhKJD+UlhOedcxPN7HLn3PhK7nsr0B8YXnKFPz/zWufcr/7yv4GTgCcqeQyRnQqf5P64g7px57iIeZ7Iee8jtvc6sox3i9Qt5SWE282sPjDEzHaYRtM591xZb3TO5QF5Zlba6rZAKGw5COxefrgiFRc+yf3zwIVhk9zn774Ha+f8FLvgROJQeQnhJuBkoAWlj25aZkIox7YSy0mUM0dzenpjUlISbzy9QCAt1iFUu3ip05133knDYHCHP6zHr7mGax57jEAl9xcv9apOqlPtURP1Km+CnLeBt83sLP91dVkJhN/R2w1YsbM3ZGdvqsbDx4dEfHIlnup01oMPRYzAOAk4Aeg9/0fOrmSM8VSv6qI61R5VqVdlEkl5TxmNcM6NAk43s9NKrnfOXVj58MA595uZpZpZB7ybyn/BGxVApMrqzZ9Hyz5HRrQA9uCPXxzqaSxSuvIuGRVNkzkNaIp3qScb7xLPTplZd7wbXEquAAAQS0lEQVTHSvcCtpvZX4F/AYudc/8EBuM9kloIvFR0g1mkKlp17kDyupzi5RcbNuTCLVuKl9XTWKRs5SWEKWb2LnAIMAcoAA4HPgEG7uyNfse1Y3ey/iugW2WCFSlL6pQvaHF2ZCN29S+/ckLjJgzVJPciFVLuU0bAb8Bf/aeG8J86ygLuAEZGNTqR8hQUENitRUTRxltuY9PgmwA0yb1IJZQ3uN1JwOCiZADgT6M5GG88I5GYafDmazskg9Bvq4uTgYhUTnkJYVN4MijinNsObIhOSCLl2LaNQEYzml17RXFR7rinCAVzoX79GAYmUruVlxC272RdQXUGIlIRjR7JItCudURZ6Pcctv7tvBhFJJI4yruHcISZLSulPAloXUq5SLUqGnpi7tSv+Hj6NxHrct54l+3Hlhx/UUR2VXkJodRxJ0RqQtHQE7cHg4QPRVeQlsaahctjFpdIoiqvp/LSmgpEpKRn776NVcFgRFk34MSrBqGeBCLVTw9kS1xqfvZpjJzyRfHyt0BP/3UzTXIvEhVKCBJX6v3saHlUj4iyPYHwG1kaekIkOsp7ykikxrTs3iUiGWw+82zaZGREJAMNPSESPWohSMylzphGi1Mj+zmudksoTG/JHP8pIw09IRJ9+pclsVNYSKBN84iijTcMZdPNf8zcWjTJvYhEnxKCxET999+j+cABEWWhpaugUaMYRSQiSghSs/LyCOzeMqJo/ZiH2XLhJTEKSESKKCFIVIVPcn9jaiqnTf48Yn1oZTbUS7ypUUVqIyUEiZqinsYbgkE2lli37sXX2db35JjEJSKl02OnEjVZWZkMLZEM8oBhN92sZCASh9RCkKhICoUYnTkqoqwnXo/j3t/OiElMIrJzaiFItWt28QW0/tM+xcs/4A2P+62/rJ7GIvFJLQSpNsmLFkJGNxqElR3esiUz164tXlZPY5H4pRaCVIv0Yw6n1eHdipe39vszoWAu/5z/C0OHjqB372MZOnQEc+YsUE9jkTilf5lSJSmzZ5HeL3KSmtU/LKQwEPDWq6exSK2hFoLsmsJCWrdpHpEMNl15LRQWFicDEald1EKQSqv/yUc0H3BORFlo8Upo0oTGMYpJRKpOCUF2qqCggMmTP2Pu3P9yUJcunPv3cyPWb7jnfjZfcU2MohOR6qSEIGUqKCjgoovO45NPPuJy4NwS60Mr1oJuEIskDN1DkDJ9/vmnTPnkIwqBp8PKv/q/GwgFc5UMRBKMEoKUqeG9d7C5RFkS8NjvK2IRjohEmRKC7CApey2BjGac9dOPxWVH4yUDgMLCmIQlIlGmhCAR0q6+jNa2V/HyUrxEMDVsm9NOO6OGoxKRmqCEIAAkL1tKIKMZDd9+o7hs9ZQZXFliVNK+fU/mhBP61nR4IlIDdFdQaHFyH1K/m1W8vO2Y41j35nsAPP/8q0ye/Bnz5s3lwAO70qfPiSQn63eESCJSQqjD6s2bS8vjj4ooWzPXUbBb2+Ll5ORkTjihr1oFInWAfurVUa06tY9IBpsvHkgomBuRDESkblELoY5J/XIyLf52ekTZ6oW/UZjWLEYRiUi8UEJIYOET3Pfq2YsHxtwfsX7DyDvZ/H9DYhSdiMQbJYQEVTTBfTAYZADwwNdTItaHflsN9evHJjgRiUu6h5CgsrIyyQ4GKQReCCt/tf8p3rATSgYiUoISQoLa583X2FaiLBl4csP6WIQjIrVAVC8ZmdldwPFAQ+BK59yssHVLgF+BfL/oAufc8mjGUxck5a6jdaf2XBlWdgIwyX+tCe5FpCxRSwhmdhzQwzl3pJl1AR4HepfY7GTn3IZoxVDXNL3xehq9OKF4OScpifSwgYc0wb2I7Ew0LxkdB7wH4JybD+xuZppQKwqSV64gkNEsIhmsnTSVzcvXaIJ7EamwaJ4d2gL/DVsOAW2AxWFl482sPfANcLNzTuNoVlLzs06l/tdfFi9v734oOR9NBrwvVxPci0hFRTMhlLynmQSEn/BvAz4DgsDbwNnAG5QhPb0xKSn1qjvGmAsE0nbtjT/8AF26RJYtWULqnnsS6ynud7lOcS4R66U61R41Ua9oJoSVQEbYcgBYVbTgnCt+GtLMPgH+tLOdZWdvqu74Yi4QSCMUqvxTPy0P3p96K/64/77lr+ew/vHx3sIu7K867Wqd4l0i1kt1qj2qUq/KJJJoJoSPgHuBJ8zsEGCRc24zgJmlAf/Gu6m8CTgKeDeKsdRK4T2NDz/8CIb3OpJWZ/4lYpvVbgmF6S1jFKGIJJKoJQTn3Hdm9l8zmw3kAQPN7GJgnXPun2b2BjDVzDYB3wNvRSuW2ii8pzHAVyV6Gm+8YRibbh4Zi9BEJEFF9ZET59xwYHhY0bywdY8Bj0Xz+LVZVlYmwWCQs9gxU4aWBaFhw1iEJSIJTM8gxqn/TP+Gko9cXQ783PtY3lIyEJEoUEKIQ42efpwvpn4VUVYPKACGqqexiESJEkI82bCBQMfdI4pOwbv7DuppLCLRpcHt4kSTf4yISAaF9eqxcvka9ldPYxGpITq7xFhSMEjrLp0iyrI/mkRe9x7qaSwiNUothFg6/fSIZJBn+xEK5pLXvUcMgxKRukothBhIXrSQVod3iyhbO2M2+R07lfEOEZHoUwuhhqX37hmRDLb2+zOhYK6SgYjEnFoIURQ+9MTZHfbkupdfiNxg1SpykxrFJjgRkRKUEKIkfOiJkh3MNl15LRvvHuUNOpWAA3GJSO2kS0ZRkpWVyWGlJINbr7+RjXePiklMIiI7o4QQDfn5jM4cxfthRdfjTQgx9fvvYhSUiMjOKSFUs4YvTCDQNj2iLAV4xH+tSe5FJF7pHkJ12byZwJ5tIooGNmvGc7m5xcsaekJE4pkSQjVofP/dNHkoM6IstGod9+Tn0york5kzp9OzZy+GDBmqoSdEJG7p7FQFSdlraW17RZTlvPcR23sdCUBKSoqGnhCRWkP3EHZR2tWXRSSD/D3aEQrmFicDEZHaRi2ESkpetpRWhx4YUbb262/Jt/1iFJGISPVQQtiJkpPc3zPpU+rP/uOx0W29j2PdW+/FMEIRkeqjhFCG8J7GBwGjS0xyv2auo2C3trEJTkQkCnQPoQxFk9xnA3PCymcc1I1QMFfJQEQSjloIZcj76N87DDvRDOjWvDlvxSIgEZEoUwuhpIICAhnNeGj+vOKim/GGnViPehqLSOJSQgjT4PVXCOzWIqKsPvCA/1o9jUUkkemSEcC2bQTatY4oyh33FBvPPJvB6mksInVEnT+7NXr4QZree2dEWej3HEhO1iT3IlKn1NmEkJS7jtad2keU5bz5HtuPOS5GEYmIxFadvIfQ9MbrI5JBQfMW3rATSgYiUofVqRZC8soVtDoocoiJtZOmkn9g1xhFJCISPxI6IYQPPfHs0qV0XrakeN327j3I+WhS7IITEYkzCZsQioae2B4MsrbEujXfzaegfYeYxCUiEq8S9h5C0dATr4aVvQgMGzpCyUBEpBQJmxBmzJgGwDnAc0BL4EJg5szpMYxKRCR+JWxCOPzwIwBYBwwEsv1yDT0hIlK6hE0IQ4YMJSMjI6JMQ0+IiJQtYW8qp6SkMGfOArI09ISISIUk9NlRk9yLiFRcwl4yEhGRylFCEBERQAlBRER8SggiIgIoIYiIiC+psLDkVPIiIlIXqYUgIiKAEoKIiPiUEEREBFBCEBERnxKCiIgASggiIuJTQhARESDBRzuNJjM7FngT+MEvmgfcDbwAtAB+Ay5wzm01szOAYUBD4FHn3HNmVg94HOgCJPnbLjazfYHxQGNgFnCNc67QzK4GBvjltzjnPqzm+nQB3gOynHPjzCyjJutiZk2ACUA7YCNwjnOu5HTY1VGviUB3YI2/SaZz7oPaVC8zuw84DkgFHgCmUMu/q1Lq9Gdq8fdkZo2BiUAboAlwFzCDOP+e1EKominOuWP9/64DMoEJzrnDgSXABWaWBowB+gFHAsPMrCnejJ4FzrkjgfuAO/19jgeGOed6AAHgODPbB7gSOAboC4wxs6TqqoT/h/MoMCmsuKbrMgz4zjl3BN4JfHCU6gUwIux7+6A21cvMegMHO+d6AScBWdTy76qMOkEt/p6AU4FZzrljgLP8uOP+e1JCqF7HAv/yX7+H9+X0AP7jnFvnnNsEfAMcjfdr6D1/24+BY82sPtDJOTezxD6OAT52zm13zq0CVgL7VWPcW4H+wIoY1iV8H0XbRqNepalN9ZoG/M1/nQPUB/pQu7+r0upU2rmp1tTJOfeac260v9gOr0VwLHH+PemSUdUcYGYfAWl4GTzNObfZXxcEdgPaAqGw9+xQ7pzL85uIAWBtKdvmlrGPn6qjEs65PCDPzMKLa7ou4fsuKotGvQCuM7OhwO/AtbWpXn6dNviLlwEfAqfW5u+qjDpBLf6eipjZTH9f/YGv4/17Ugth1/0PuAfvWucA4Bm8a31FkoBCYFuJ95VVDpBXwW2LyqMp/Jg1UZeS5dGq34vArX5T/j9413ZrXb3M7DTgcrzLAAnxXZWoU0J8T865nsAZwGslYorL70kJYRc555Y7515xzhU45xbj/Ypp6t9MAi8br8BrvmWEvXWHcr8puB0vm7fY2bYlyqNpfQ3XJby8LVGqn3NuknNulr/4AfCnisQUT/Uys77AbUA/51wOCfBdlaxTbf+ezOxQM+sA4JybjXeu3Rjv35MSwi4ys3PN7A7/dWu8pwmeAU7zNzkT7w/5W6CrmTX3bxb1BL4GPgJO97f9C/CZc64AmGNmvfzyM/x9fAb0NbNUM9sdaOmc+znKVfy4husSvo+i41U7M3vDzLr6i0cB82tTvcysOfAQ0N85V/QETq3+rkqrU23/noAj8G/imlkbvMvK/ybOvyfdQ9h1/wbOMbNv8BLrNcD3wKtmdgPggNf963+3AV8BBcBdzrnNZvYucKqZzQI2Aef7+x0OTDCzFOBL59w3AGb2HF7TuYBqeAInnJl1Bx4E9gK2m9lfgQuAl2uwLk/5n90svOud50SpXrcBz5jZJmA9cKlzblstqtc5QDrweti9kYuA52vxd1VanWr79/SUf+yvgQZ454fvqNnzQ6XrpPkQREQE0CUjERHxKSGIiAighCAiIj4lBBERAZQQRETEp4QgUgoz28vMfqvmfZ5gZl9W5z5FqpMSgoiIAOqYJlIpZnYpcDWwGa+zz0Dn3DozuxxvYLYt/n/nOOdyzOx04F68YQNcjMIWqRC1EEQqyB+b5i7geOdcb2AZcIO/uhFwmj8Y21Lg7375OOAs59yJQH4NhyxSKUoIIhV3CN7Y9bn+8iTgUP/1FuBNM5uCN+58azNrBTRyzi0I214kbikhiFRcyVnqkoACM9sbb1arv/kthH+Wsb3+vUlc0x+oSMX9BzjUzJr5y/3w5slNB0LOuZX+yLcn4g1otgZvUL2iEduqYxY4kajR4HYipTCzvfBuAk8PK/4Wb5a6q/Gm5/wVb0KXzcA7eGPOLwTexLt3cA7e9Im3AYuABcBh/v0HkbijhCAiIoAuGYmIiE8JQUREACUEERHxKSGIiAighCAiIj4lBBERAZQQRETE9/9LVJHePuxDPgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x6adce2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X=np.vstack([x,np.ones(len(x))]).T\n",
    "a=np.linalg.lstsq(A,y)[0]\n",
    "\n",
    "print('Least squares estimate values:', a)\n",
    "\n",
    "plt.plot(x, y, \"k.\",label='Original Data',markersize=12)\n",
    "plt.plot(x,np.dot(X,a),'r',label='Fitted Line')\n",
    "plt.xlabel(\"Load\")\n",
    "plt.ylabel(\"Diflection\");\n",
    "plt.legend()\n",
    "plt.title('Least Squares Estimate for Linear Regression')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After fitting a straight line to the data, many people like to check the quality of the fit with a plot of the data overlaid with the estimated regression function. Graphical analysis of the residuals is the single most important technique for determining the need for model refinement or for verifying that the underlying assumptions of the analysis are met."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x6aca20f0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEVCAYAAAAo63jjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+cXFV9//HXsgFikiUEsglJU0iV9s2volAgxhRIpAoaDD8spYICFvzm66+qfL/fqNViTSm1gtoKpWJTQKWCxsqPGsFSlB+lafghWiXyUSuBxMSwQAghCeTH7vePezcOk5ndmd07d+fefT8fDx7MnHvPzPnsZOez55x7zu3o6+vDzMwsC3uMdAPMzKw8nFTMzCwzTipmZpYZJxUzM8uMk4qZmWXGScXMzDIzZqQbYFaPpD7gf4AdadEY4F7gTyPihSG+5mPAiRGxvqr87cBFETF3iK87F1gSEQcPpX47kTQT+HlEtOT7Ybg/a2tv7qlYu5sbEYdExCHA4cB+wJ8P9cXS11o/+JntTVLnSLfBrBb3VKwwIuIlSXcApwJI2gu4AjgZ2AncCFwaEX2S3ge8F+gAngfeGRGPpr2f3wTWAp8HFgC/Au7pfx9Jd5P0Om6ofi7pOODvgX1IelAfiIh/r2ynpCOAfwQmAnsCfxcRV1Wd8yDw1xHxzfT5GcCHI+K1khYAlwJ7kfTULoyIX0n6i7TtrwZulHQj8BVgGrA3cFNEfCw9b0ZEXJS+9q7nkk4EPge8Im3KJRGxtNHPQNJ+wBeA1wC9wLUR8en02JuBT6VteSFt9w8k7VHvZ23l456KFYakfYG3Af+ZFr0fEHAEcCxwBvBmSV0kX8rHpT2cTwPzq17uFOCNwKHAicAJDTbji8DfRoSAvyH5gq32CeALEXEYMBs4SdLeVed8Azit4vnpwNcl/QZwPXBW2vb/AK6savebIuIzwIeAe9P3+V3gtyRNG6T9VwAfiohDSZLzGYMFXOUy4NmI+J00tvdIOj7tOX0FeHf6s7k1fa/+Ng/lZ20F5J6Ktbu7Je0g+at9P+DvSJIEJF/KV0XEdmC7pC+TfEl+B9gOXCjppoj4Ro3XPQG4PSI2A0j6OvCWBtpzLMlf6AD3Aa+scc464K2Sfgx8PyJqfXEvBVakX8YdJEnv4yS9rv+KiEjP+wfgKUkd6fMVEdFT8T5/KOkuYHlEnJPGMlD71wHnSXoqIn4CnDNoxC83nzQRRcQGSf8KvDEi7pM0PSJeSs+7D7ggfTzUn7UVkHsq1u7mpn+xH0fyZf7ViOifuJ8CXCHpsXQC/gPAhPT460n+kn5M0n2SDq963f2ADRXPn22wPe8AHpQUwJ0kCaHaIuC/gZuANZL+d/UJEfELYA3wOpIv3YiI1WlMcypi+i9gI7B/jXb+LfBN4GpgvaRPNND+C0iGpu6U9DNJZzZQp9IU4JmK588C3enjD0n67/Rncz2//n4Z6s/aCsg9FSuEiHha0udJein9w0a/Aj4bEbfVOP9HwNmS9gQuJhmmOr7ilA0kcx79plQ83snLk0UX7Loq6gvArIh4RNLBwM9qvPeLJL2Oj0s6Cvg3Sf8eET+vOvUbJPMMY4GvV8R0V0Ts9mVf3QOJiF6SIaYrJL0K+Lak/6jX/rTOsyTJ9wOSXg/cIuk7/b2IBjwFTAYeT59PJkloJwD/h2TI8XFJfwAsSc8Z6GdtJeOeihXJZ4DXpZPNALcB75K0p6QOSR+XdLKkIyUtlbRXOjT2A3bvUSwHTpH0CknjgLMqjv0SOARA0qHAYWn5RJK/8n+SJqv3pedMqHxhSd+q6Bk9RnKhQC1LgT8A3kySYAD+DTg+TVhIOk7S52pVlvRFSW9In64i+cLvqGr/eOBN6eM9Jd1TMe/y3yQXG/TSuH8F3pm+3mSSpLgMmASsB55Ifx7vAsanw3YD/aytZNxTscKIiE2SPkXyl/lxwFXAQcCPSP4tryAZEtoCPAGslLQNeI7kSrBK/0oyPxAkX8bLSIbMIEle/yxpDvAT4N+BTpIv4W8DPwWeBv6UZI7l34A/q3jtK4Gvplen9ZHM+1T3UoiIn6ZXRq2NiF+mZWslXQj8Szq5v5nkgoRargH+QVL/RP5twF3ABOAdkh4gGWJbBkyOiO2SlgB3pV/2O4D3RcTWGq/dmQ6/VXoT8DHgi5J+SpKMPhURD0gaS/Iz/gXJvM37gKOBrwJvp/7P2kqmw/dTMTOzrHj4y8zMMuOkYmZmmXFSMTOzzDipmJlZZkbV1V89PZtKd1XCpEnj2LBhy0g3I1NljAnKGVcZY4JyxjWcmLq7u2ot8q3JPZWCGzOmfJvVljEmKGdcZYwJyhlXXjE5qZiZWWacVMzMLDNOKmZmlhknFTMzy4yTipmZZWZUXVJsVkQrVq5n2fJVrH16C9Mnj2P+7JnMOmzqSDfLrCYnFbM2tmLleq657dFdz9f0bN713InF2pGHv8za2NLv7bZjflJ+d+1ys5HmpGLWxp7d9FLt8udrl5uNNA9/mZWc52QsT04qZm1sv669a/ZW9uvau6H6npOxvHn4y6yNnTXv4KbKq9Wdk6lTbjZc7qlY6Y308M+u939mC9P3b+79+89btvwJ1j2zmWn7j2f+7IMarl93TqZOudlwOalYqY308E8W7z/rsKkeqrLC8PCXldqy5avqlD+Ry/uP9CXB9eZe9tunsTkZs2Y5qViprX269k2J1j2zOZf3H+lLguvOycxtbE7GrFke/rJSmz55HGt6dk8g0/YfPwKtyd9w52TMmuWkYqWmAyfVTCo6cN+GX2M4E/3DvSQ4C56TsTzlnlQkLQZOAsYCCyPioYpjs4HPpMe+GRGX1qsj6TeA64C9gZ3A2yNiba7BWNuLJzfUKX+uofrDnWg/a97BL6tfWW5WRrnOqUiaBxwbEXOA84HPVp3yJeBs4BjgLZJeNUCdvwSWRMSJwDeAi/OIwYpluHMqw53on3XYVBYuOJwZ3RPo3KODGd0TWLjgcPccrLTy7qnMA24FiIgfS5ouaVxEbJH0SuDZiFgNIOlbwBuBabXqAH8KvJi+7tPAa3KOxQpguHMqWUz09w8/dXd30dOzqeF6ZkWUd1KZBvyw4nkPMBV4PD3WU3HsKWB6vToR8TiApE7gvcBfDPbmkyaNY8yYzmE0vz11d3eNdBMyl1VMbzv5EC6/4eEa5WroPfabOJann9u6e/k+Y4fUxiJ+Vvc+soald/2MJ9dv4sCpXZx10m9zwlEzdh0vYkyNKGNcecSUd1LZVvW8A+gb5FjdOmlC+Qpwd0R8b7A337Ch9l+dRVbGv36zjOnQGRNZuODw3a5+OnTGxIbeo3dnb83ynTt7m25jET+r6jmlVeue5/IbHub5518sde+rjHENJ6ZmklHeSWUdMKXieTewvs6xA4C1wI4B6lwHrIqIT7SktVYKw7n66bkXqv+mSWzcXLu8bAaaU/K8kNWS9+LH24HTACQdDfwiIrYCRMQaYE9JB6Y9kFPT82vWkXQu0BsRf5ZzDDaKTJ88rmb5aFnnMtKLR614cu2pRMTDkn4o6fskPZALJV0AbIyIm4EPkkzK9wE3pJP2q6vrpC/3XmCspLvT5ysj4j05hmOjwPzZM2teEjx/9kEj0Jr8jfbFo9a83NepRMSHgQ9XFP2o4ti9wFEN1CEiXteqNpr1G+0r0rNYPGqji1fUmw1iNK9IH+7iURt9vKGkmdXlORVrlnsq1vZG+iZbo5nnVKxZ7qlYW+tfJ7GmZzO9fX279t5asXL94JVt2ObPnlmnfHRcqGDNc0/F2prXSYys0X6hgjXPScXamsf0R95wL1Tw8OXo4qRibc1j+sU23FsHWPF4TsXamsf0i224tw6w4nFPxdqax/SLzcOXo4+TirW90bz4sOj2nbBXzdspT5yw1wi0xvLgpGIt54la203f4KdYMTmpWEt5onZ0G+23DhiNnFSspbJYZ+KeTnH56r3Rx1d/WUsNd6LWK+qLzVfvjT5OKtZSw73JlS9JLbZZh01l4YLDmdE9gc49OpjRPYGFCw53T7PEPPxlLTXcm1z5ktTi89V7o4uTirXUcNeZeEzerFicVKzlhvOX6mi/na9Z0TipWFvzinrz1X/F4qRibc9j8qOX1zkVT+5JRdJi4CRgLLAwIh6qODYb+Ex67JsRcelAdSS9H/gsMCkiXsg1EDNrOd9Pp3hyvaRY0jzg2IiYA5xPkhAqfQk4GzgGeIukV9WrI+k84ABgbV7tN7N8+eq/4sl7nco84FaAiPgxMF3SOABJrwSejYjVEdELfAt44wB1bo6Ij+FdhMxKa7jrnCx/eQ9/TQN+WPG8B5gKPJ4e66k49hQwvV6diHi82TefNGkcY8Z0Nlut7XV3d410EwZ07yNrWHrXz3hy/SYOnNrFWSf9NiccNWPAOu0e01CVMa5WxvS2kw/h8hserlGulv8s/VkNTd5JpXoXuQ5+3dOod2ygOk3ZsKF2V7rIuru76OnZNNLNqKt6onXVuue5/IaHef75F+uOibd7TENVxrhaHdOhMyaycMHhu139d+iMiS19X39Wu9dtVN5JZR0wpeJ5N7C+zrH++ZIdA9SxNueJVhsuX/1XLHnPqdwOnAYg6WjgFxGxFSAi1gB7SjpQUidwanp+3TrW/jzRaja65NpTiYiHJf1Q0vdJeiAXSroA2BgRNwMfJJmU7wNuiIjVwOrqOgCSPga8gaRHc7uk5RGxKM94bHDeZsVsdOno6xs9F0/19GwqXbDtPvZbPafSb6Cdats9pqEqY1xljAnKGdcw51Q6Gj3XK+qtpbzNig2Xt2kpFicVazlPtNpQeZuW4vFNusysbfkmbcXjpGJmbctXDxaPk4qZtS1v01I8Tipm1rbmz55Zp9w3aWtXnqg3s7blqweLx0nFzNracK8e9CXJ+XJSMbPS8iXJ+fOcipmVli9Jzp+TipmVli9Jzp+TipmVli9Jzp/nVGxQnui0opo/e2bNDU19SXLrOKnYgDzRaUU267Cp/PyXG7n3B2vZvrOXPTv34ITXTPe/3Rby8JcNaOn3ft5UuVk7WbFyPXc9vIbtO3sB2L6zl7seXsOKlb55bKs4qdiAnt30UlPlZu3EV3/lz0nFzErLV3/lz0nFBrRf1961y/epXW7WTnz1V/6cVGxAZ807uHb53NrlZu3EG1Lmz1d/2YC8oZ8Vmf/95s9JxQbl2wFbkfnfb75yTyqSFgMnAWOBhRHxUMWx2cBn0mPfjIhL69WRNAX4MrAvsAY4NyJ8SZKZ2QjKdU5F0jzg2IiYA5wPfLbqlC8BZwPHAG+R9KoB6lwOXBcRrwVWAefmEIKZmQ0g74n6ecCtABHxY2C6pHEAkl4JPBsRqyOiF/gW8MYB6swFbktf91bg5BzjMDOzGvIe/poG/LDieQ8wFXg8PdZTcewpYPoAdboiYmvFuQcM9uaTJo1jzJjOITe+XXV3d410EzJXxpignHGVMSYoZ1x5xJR3UtlW9bwD6BvkWCPlla9T14YNtRdCFVl3dxc9PZtGuhmZKmNMUM64yhgTlDOu4cTUTDLKe/hrHTCl4nk3sL7OsQOAtQPU2dQ/dFZxrpmZjaC8k8rtwGkAko4GftE/hBURa4A9JR0oqRM4NT2/Xp07+suBM4FleQZiZma7y3X4KyIelvRDSd8HdgAXSroA2BgRNwMfJJl07wNuiIjVwOrqOunLXQbcKOliIICv5RmLmZntrqOvb9CpiNLo6dlUumA99lscZYyrjDFBOeMa5pxKR6PnNtRTkXQ6sB9wHUlP4veBD0TEV4bUQjMzK6VG51QWAUtJ5jl6gd8F3tOqRpmZWTE1mlS2RsQm4BTg+oj4JbB1kDpmZjbKNDpRP1bSImAB8GFJBwO1b1RgbWfFyvUsW76KtU9vYfrkccyfPdMb7JlZSzTaU/lfJKvbz4+IF0h6LB9tWassMytWruea2x5lTc9mevv6WNOzmWtue9T36Dazlhiwp5LuxwXJUNfnK8q+3eJ2WUYGuke3eytmlrXBhr/uIlkzUutysj7glTXKrY34Ht1mlqcBk0pE/Fa9Y5LmZN8cy9r0yeNY07N7AvE9us0a4znJ5jS6TmUf4O3A5LRob+A84Ddb1C7LiA6cVDOp6MB9R6A1ZsXSPyfZr39OEnBiqaPRifqvAa8G3glMItlz632tapRlJ57cUKf8uZxbYlY8A81JWm2NJpWxEbEQeCIiPkRyg6yzWtYqy4znVMyGzr8/zWs0qbwivSf8HpKmRMTTwKta2C7LyPTJtZcTeU7FbHD+/Wleo0llCcm9478CPCrpUeBXLWuVZWb+7Jl1yg/KtyFmBeTfn+Y1NFEfEUv6H0u6E5gUEY+0rFWWmf7JxGXLn2DdM5uZtv945s8+yJOMZg3w70/zGr36a3GNsrdGxMezb5JlbdZhU/1LYDZE/v1pTqPDXzsr/usAjuPlt/g1MzNrePjrk9Vlkq7KvjlmZlZkQ7pHvaQxwGEZt8XMzAqu0TmV1SR7fUEy/LUPyV0gzczMdmn0fiq/X/G4D9iQ3rTLzMxsl8G2vj9vgGN7RMT1zbyZpE7gauAIkh7PuRHxeNU5Z5DcvngscGVEXFuvnqQ9gL8CLoqI7mbaYmZm2Rusp/KG9P+TSe5Lv4JkHmY2cD9wfZPvdx7QGxFzJJ0KfDItA0BSF3AFcDSwHfi+pK+TbAlTq95HgNXU3prfzMxyNtjW9+8AkPQNQBGxOX3eRbLKvlnzgK+mj+8g6X1UOhZ4MCI2pu9zP3D8APWujIhNtdbRmJllwVvfN6fROZUD+xMKQPpFPpR9CqYBPelr7JDUKakzInZWH089BRwwQL2m5nUmTRrHmDGdQ2h2e+vu7hrpJmSujDFBOeMqY0yQxHXvI2tqbn2/zz5jOeGoGSPYuqHJ47NqNKk8Kuk/gOXp81nAzweqIOki4KKq4lfXOLWv4vG2qmMd6fHq8up6DdmwofaOo0XW3d1FT0+5rpkoY0xQzrjKGBP8Oq4bv/NYzeM3fic4dMbEnFs1PMP5rJpJRo0ufnynpDeQzKt0AH8NfGeQOkuoGiKTtIR0Jb6kvYDtEdFbcco6Xr5S/wDgu5XldeqZmWXOW983b8DFj5KOSv//epItWn4APAK8RHJPlWbdDpyePj4VuLPq+APAkZImSppA0iO6r4F6ZmaZ89b3zRusp/IOkiTy5zWO9ZH0IppxC7BA0kPAFuAcAEkfAe6JiOWSLgHuBXqBxRGxVVK9eleS9J4mSrobuC0iPttkm8zMapo/e+bL5lR+Xe6t7+vp6OtrbmoiXRsyISKeb02TWqenZ1PT8zDtroxj2mWMCcoZVxljgpfHlVz9Vfyt74c5p9Lwso1Gt2lZCOwJfJFkrcrBkhZHxOVDaqGZWUF46/vmNLqh5AXAF4AzSIbD9gfObFGbzMysoBq9pPjFdH3IKcBNEbFN0tZWNszMrAxG2+LJRpNKn6SrgdcDCyW9Ftirdc0yMyu+FSvX11w8CZQ2sTQ6/HUu8FPgTRGxDZgJvLtVjTIzK4Nly1fVKX8i13bkqaGkEhHrgCeBk9OiR4CVrWqUmVkZjMbFkw0lFUl/C/wxSY8F4A+Bv29Vo8zMymA0Lp5sdPjr6Ij4I2ATQET8FbX38TIzs9T82TPrlJd38WSjE/X9uwj3wa6bbQ3p/vZmZqNF/2R8GRZPNqrRpPKgpH8Cpku6mGS9yvda1yyrNNouSTQrk9G2eLLRpPI3JDfKeh6YAXyOZLLeWmw0XpJoZsU12C7FJ0j6FRDAXwJXRcTFwHSS3YOtxUbjJYlmVlyD9VT+GjgpIh6VNB/4Yrqh5AbguJa3zkblJYlmVlyDTbbviIhHASJiGXAQ8PmIODMi1ra8dca+E2pvXDBxvDc0MLP2M1hS2Vn1fHVE3NyqxlgTGt6I2swsP81eFly6+5G0u+de2FazfGOdcjOzkTTYnMrrJD1Z8XxK+rwD6IuIA1vXNINkRe6ant3nT8q8ItfMimuwpKJcWmF1+XamZlYkAyaViPB1qyNsNK7INSuT0bZ4udHFjzaCRtuKXLOyGI2Ll71/l5lZi4zGxcu59lTSjSivBo4gmew/NyIerzrnDGARMBa4MiKurVdP0hFpOcBG4G0R8UI+0ZiZDWw0Ll7Ou6dyHtAbEXOAy4BPVh6U1AVcAZwCzAEWSZowQL3PA/83Ik4AfgJckEcQZmaN8P1UWm8ecGv6+A5gbtXxY4EHI2JjRGwB7geOH6DemRHxQPr4aWCf1jTbzKx5vp9K600DegAiYoekTkmdEbGz+njqKeCAAeo9ByBpPElv5q0DvfmkSeMYM6Yz04DaQXd310g3IXNljAnKGVcZY4Js4jr1xC722WcsS+/6GavXb+I3p3Zx1km/zQlHzcighc3L47NqWVKRdBFwUVVxrbtFVq7Sr14m3pEer7V8vP+GYeOB24DPRUQM1KYNG2qPbxZZd3cXPT2bRroZmSpjTFDOuMoYE2Qb16EzJnLJ+ce8rGwkfmbDiamZZNSypBIRS4AllWWSlgBT0sd7AdsjorfilHX9x1MHAN+tLK+sJ2kMSUL5ekT8U6tiMTOzxuQ9p3I7cHr6+FTgzqrjDwBHSpqYTtDPIrlvS716Hwbui4hrWtpqMzNrSN5zKrcACyQ9BGwBzgGQ9BHgnohYLukS4F6gF1gcEVsl1awHvBdYJWlu+vy7EbE4v3DMzKxSR1/f6Nl4uKdnU+mCLeOYdhljgnLGVcaYoJxxDXNOpeGbbXhFvZmZZcZJxczMMuMNJc3M2ljRdjl2UjEza1NF3OXYw19mZm2qiLscO6mYmbWpIu5y7KRiZtamirjLsZOKmVmbKuIux56oz0HRrt4ws/bQ/z2xbPkTrHtmM9P2H8/82Qe19feHk0qLFfHqDTNrH7MOm1qo7woPf7VYEa/eMDMbKieVFivi1RtmZkPlpNJiRbx6w8xsqJxUWqyIV2+YmQ2VJ+pbrIhXb5iZDZWTSg6KdvWGmdlQOamYmbWxoq1zc1IxM2tTRVzn5ol6M7M2VcR1bk4qZmZtqojr3HId/pLUCVwNHAF0AOdGxONV55wBLALGAldGxLX16klaAHwU2AH0AOdExIu5BWRm1kLTJ49jTc/uCaSd17nl3VM5D+iNiDnAZcAnKw9K6gKuAE4B5gCLJE0YoN6HgFMi4nhgI3BmLlGYmeWgiOvc8p6onwd8NX18B0nvo9KxwIMRsRFA0v3A8fXqRcS89Lw9gSnAmlY23swsT0Vc55Z3UplGMkxFROyQ1CmpMyJ2Vh9PPQUcMFA9SRcAlwL/EhH3DvTmkyaNY8yYzmwjagPd3V0j3YTMlTEmKGdcZYwJ2ieuU0/s4tQTD87ktfKIqWVJRdJFwEVVxa+ucWpfxeNtVcc60uPV5bvqRcT1km4Avizp7RFxQ702bdhQe9KryLq7u+jp2TTSzchUGWOCcsZVxpignHENJ6ZmklHLkkpELAGWVJZJWkIyTIWkvYDtEdFbccq6/uOpA4DvVpb31wP2kvT6iPh22nu5mWSYrG5SMTOz1sp7ov524PT08anAnVXHHwCOlDQxnaCfBdxXp9424B8kTUvLXwtEC9tuZmaDyHtO5RZggaSHgC3AOQCSPgLcExHLJV0C3Av0AosjYquk3epFRK+khcAtkraR9GY+nnM8DSnaNgtmZkPV0dfXN/hZJdHTsyn3YKu3Wei3cMHhmSQWj/0WRxnjKmNMUM64hjmn0tHouV5R32JF3GbBzGyonFRarIjbLJiZDZWTSov5dsJmNpo4qbRYEbdZMDMbKt9PpcWKuM2CmdlQOankwLcTNrPRwsNfZmaWGScVMzPLjJOKmZllxknFzMwy46RiZmaZcVIxM7PMOKmYmVlmnFTMzCwzTipmZpYZJxUzM8uMk4qZmWXGScXMzDLjpGJmZpnxLsUNWLFyPcuWr2Lt01uYPnkc82fP9K7DZmY1OKkMYsXK9Vxz26O7nq/p2bzruROLmdnL5ZpUJHUCVwNHAB3AuRHxeNU5ZwCLgLHAlRFx7WD1JC0EPhoRM7Nu87Llq+qUP+GkYmZtb9dIyzNbmL5/60da8p5TOQ/ojYg5wGXAJysPSuoCrgBOAeYAiyRNGKiepCnAma1q8Nqnt9QsX/fM5la9pZlZJvpHWtb0bKa3t2/XSMuKletb9p55J5V5wK3p4zuAuVXHjwUejIiNEbEFuB84fpB6nwYuaVF7mT55XM3yafuPb9VbmpllYqCRllbJe05lGtADEBE7JHVK6oyIndXHU08BB9SrR5JwXoiIFZIGffNJk8YxZkxnUw1+28mHcPkND9coF93dXU29Vqu0SzuyVMaYoJxxlTEmKEdca5+pP9LSqvhallQkXQRcVFX86hqn9lU83lZ1rCM9Xl0O0An8OU0MfW3YUPsHPJBDZ0xk4YLDWbb8CdY9s5lp+49n/uyDOHTGRHp6NjX9elnr7u5qi3ZkqYwxQTnjKmNMUJ64pu8/jjU9uw/VT9t/fFPxNZOAWpZUImIJsKSyTNISYEr6eC9ge0T0Vpyyrv946gDgu5Xl/fWAo4DfAL6T9lKmSbopIv4461hmHTbVk/JmVjjzZ8982dWrvy4/qGXvmffw1+3A6en/TwXurDr+AHCkpInATmAW8G6gq7peRKwADumvKGlVKxKKmZk1Lu+J+luAsZIeAj5IehWXpI9Imh0R20gm3e8F7gMWR8TWevXMzKy+0k/UpxPy59co/1TF46XA0kbqVZ0zM5tW7s4r6s2siEZiSYRX1A/CK+rNrKimT64/Ud8q3lByECPRfTQzy8L82TPrlJdnor5wvKLezIqqfzSleklEK0dZnFQGMRLdRzOzrPQvichr7Y2HvwYxEt1HM7Oick9lECPRfTQzKyonlQZ4Rb2ZWWM8/GVmZplxUjEzs8w4qZiZWWacVMzMLDNOKmZmlpmOvr6+wc8yMzNrgHsqZmaWGScVMzPLjJOKmZllxknFzMwy46RiZmaZcVIxM7PMOKmYmVlmvEvxCJI0F1gKPJoW/Qj4S+DLwL7AGuDciHhJ0hnAImAscGVEXCupE7gaOALoSM99XNLvAP8IjAMeAt4TEX2S3g28Iy3/s4j4dsbxHAHcCnwuIq6SNCXPWCS4ymExAAAFA0lEQVSNB64DZgCbgbMj4tkWxHU98HvAM+kpl0fEsiLFJekyYB6wJ/A3wD0U/LOqEdN8iv85jQOuB6YC44HFwH/Rxp+Veyoj756ImJv+937gcuC6iHgtsAo4V1IXcAVwCjAHWCRpAnAe0BsRc4DLgE+mr/mPwKKIOBboBuZJehWwEDgROBm4QlJHVkGk//CuBO6qKM47lkXAwxHxOpIk8MEWxQXw0YrPbVmR4pJ0AvCaiJgNvBH4HAX/rOrEBAX+nFILgIci4kTgrWnb2/qzclJpP3OB29LHt5J8uMcCD0bExojYAtwPHE/yV9mt6bl3AHMl7QUcHBErql7jROCOiNgeEeuBdcAhGbb7JeDNwNoRjKXyNfrPbUVctRQprv8E/ih9/BywF/B6iv1Z1Yqp1vdbkWIiIm6KiE+nT2eQ9Ezm0saflYe/Rt5hkm4Hukj+iuiKiK3psaeAA4BpQE9Fnd3KI2JH2tXtBp6tce7zdV7jJ1kEERE7gB2SKovzjqXytfvLWhEXwPsl/T/gV8B7ixRXGtML6dOLgG8DC4r8WdWJCQr8OVWStCJ9vTcD97XzZ+Weysj6GXApydjvO4AlJOOe/TqAPmBbVb165QA7Gjy3v7yVKt8zj1iqy1sV31eAj6VDEg+SjHMXLi5JpwHvIhnOKMVnVRVTKT4ngIiYBZwB3FTVrrb7rJxURlBE/DIivhoRvRHxOMlfUxPSyTlI/iJYS9INnVJRdbfytEu7neQvin0HOreqvJU25RxLZfk0WhRfRNwVEQ+lT5cBhzfSpnaKS9LJwCXAKRHxHCX4rKpjKsnndIykAwEi4vsk39mb2/mzclIZQZL+WNJfpI8nk1zhsQQ4LT3lTJJfhgeAIyVNTCffZgH3AbcDp6fnngrcGRG9wA8kzU7Lz0hf407gZEl7SpoO7BcRP21xiHfkHEvla/S/X+YkfV3SkenT3wd+XKS4JE0EPgu8OSL6r4wq9GdVK6aif06p15FOjEuaSjJM/i3a+LPynMrI+hZwtqT7SRL8e4BHgBslXQwE8LV0LPQS4F6gF1gcEVsl3QIskPQQsAU4J33dDwPXSRoD3B0R9wNIupZkGKCXbK5M2UXS7wGfAWYC2yX9IXAu8M85xnJN+rN7iGTs9+wWxXUJsETSFmAT8CcRsa1AcZ0NTAK+VjFXdD7wpQJ/VrViKvrn1P+a10m6D9ib5DviYfL9jmgqLt9PxczMMuPhLzMzy4yTipmZZcZJxczMMuOkYmZmmXFSMTOzzDipmLWIpJmS1mT8mn8g6e4sX9MsS04qZmaWGS9+NMuZpD8B3g1sJVlMdmFEbJT0LpLNEF9M/zs7Ip6TdDrwVyTbY8QINdusIe6pmOUo3cdpMXBSRJwAPAlcnB5+BXBaugHiE8Db0/KrgLdGxBuAnTk32awpTipm+Tqa5L4Xz6fP7wKOSR+/CCyVdA/JPSsmS9ofeEVEPFZxvlnbclIxy1f13TY7gF5Jv0VyZ74/SnsqN9c537+z1tb8D9QsXw8Cx0jaJ31+Csk9xycBPRGxLt2x+g0kGwg+Q7KRZf8uiVnczdKsZbyhpFmLSJpJMrG+vKL4AZK7bb6b5FbFq0luKrUV+CbJ/Sr+B1hKMpdyNsltZC8BfgE8BhyXzseYtR0nFTMzy4yHv8zMLDNOKmZmlhknFTMzy4yTipmZZcZJxczMMuOkYmZmmXFSMTOzzPx/MP9vZ+iqGKgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x6acae790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "residuals=y-np.dot(X,a)\n",
    "plt.plot(x, residuals, \"o\")\n",
    "plt.xlabel(\"Load\")\n",
    "plt.ylabel(\"Residuals\");\n",
    "plt.title('Residuals versus Load')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Reviewing the plots of the residuals versus all potential predictor variables can offer insight into selection of a new model, just as a plot of the data can aid in selection of an initial model. Iterating through a series of models selected in this way will often lead to a function that describes the data well.<br>\n",
    "The horseshoe-shaped structure in the plot of the residuals versus load suggests that a quadratic polynomial might fit the data well. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the residual plots, the function used to describe the data should be the quadratic polynomial:\n",
    "\n",
    "$$ D=\\beta_0+\\beta_1L+\\beta_2L^2+\\epsilon$$\n",
    "\n",
    "The data with a quadratic estimated regression function and the residual plots are shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Least squares estimate values: -3.16081908282e-15 7.32059161768e-07 0.000673564732418\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEVCAYAAADgh5I1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FOX2wPFvQkINPQmKgIjIQUSkqIAFEVBQr2LHXsEG9wd4jcgV0asoKmAAsYAFFBsKKnZUVFSkiEhTORaQDgmkQCBAkt3fHzOJuyGVZLObzfk8Dw8778zOnndmM2ffd2beifB6vRhjjDGRwQ7AGGNMaLCEYIwxBrCEYIwxxmUJwRhjDGAJwRhjjMsSgjHGGACigh1AuBMRL9BcVTcH+HMGqeoLhcwbAIwAagE1gFXAXaq6NZAxlQcRaQmsB7SA2b1VdUsR783bJiKyFjhLVXeUQ0xNgK6q+kEp3zcT6APcqKqfH+ZnPwQ0U9WB+cqPAuapavvDWW9ZiMgM4AJgl1sUCewE7lbVxRUdT3FE5FTgEVXtG+xYQo0lhDAgItWAccAhCUFE2gGTgVNVdYOIRAGJwMtAv0DHpao55bCqHFVtW8rPPgK4F3eblPb9xTgb58BeqoQAXAO0UdW/SrJwabafmxgDngyKiGmSqo7xWW4A8J6INFXVMt3sJCIRAGVdTy5VXQpYMiiAJYQgEZHqwHicL2YO8CYwRlW97i+YZ4B6QDYwVFW/dA/mzwE9cPbdSuAm4H2gvvsr+DxVXe/zUccDyaq6AUBVs0XkXqCmG0dj97OPA9YAqcBGVR2Vv3XjOy0iN+G0OqKAbcD1bsK5CegPNAB+Bu4WkUHAcJxfjj8Bd6jqHhE5Cyc51XJjHa2q75RyO8YAM4G2OK2f+cBdwA9AM3ebdAAOAM0BAR4FFgKX4vyqvRMnobYFpqvqSHfd9wG3ANWA34DrgFbAFCBKRGJU9SoRuQgYA1QH/gJuVdXt+eL8xq3/PBH5P3dbvwAcg7OPx6rqTLdFtBh4CzgFOL2E26El8KeqRonIrTjfq3TgTMADDFDV1SJSD3jWXXcWzoE8txV1PvC4ux0z3HqsEJGebvlmnO/qgBKENNetQ2NgZ2HbSERigddxvn8/AruB7ar6gIj8DbwCXAv0EZGDwDR32f3AKFX9sIi/i32FlHcGXlTV1iJSE5gI9HK300fACFXNEZFNwMPAQOAo4B1VHV6Culdadg4heP6Nc3Bqj/PHeQlwvjtvGjBRVQV4AnjeLe+Lc0BqC7TGOah0xzlo5ahq23zJAOB74CgR+VBELhWRRqqaqaqp7vwROAnjGGAozsE8u6jA3T/iqUBfVT0O+BN4wGeRc4DbVfVuEekCPAT0cn+l7wFGu8uNB4ar6vHAv9xtUFo3AmnuOtrgHOROwNkmG91tctBn+RygE/AhzrbMwTloXACcBSSISB0R6Qjcg7NvjsNJoENUdTlOQpjtJoOjgBnAFW79vgeezh+kqvZ0X/ZU1U9w9vHXqtoGp6U2WUSOcZdpBKxQ1RIlgwLkuPV5zo3pa2CYO+9RnOTYFqelM1pE2rutzJnAne73bi7O/sl1EvCsqhabDNxf9HcCv6vqzmK20Qic/dcKeAy4yo0/15Gq2lpV/8b5zi1147sMmO5+Fwv7uyis3NcwnB8Kx7t1PBOnJZe7HXu67+kKDBaR5sXVvzKzhBA8/YGXVDVLVfcCr/LPAfEUnF9XAN/hfKkBtuN8cS8BaqvqaFWdV9SHuH3mnYEtwCQgWUS+FJEO7iI9gNnusuvczyuSqu4E6qvqxgJiBOdA8LtPPef6/GJ+3qee24AbROR4Vf1LVa+hYNVEZG2+f2/5rKO7iJwLVFPVu1R1RTFVSFPVb9wuiF+BBaq6V1UV8ALx7jqaqWq6qnpwWhStClhXX2Cx+15wkkv/3G6OgohINE7SfBHA3Y7f4/xKBYgG3i2mDsX51U1e4LTKmrmv+wNTVdWrqsnAHOBitxuoqaoudJfLv08zVfWrIj5vaO6+wWld9MZJSlD0NjoTmAWgqitxWna+PgJwWwHn4f44UtU/cbZZPwr/uyjJ38sFOK3CHFU9ALwNnOsz/w1V9bhdctt9tmNYsi6j4IkHxovIw+50DWCJ+/p6YIiI1MHprsjtQ/1JRO7A+VXziojMBYYU90Fuq+EOABE5HrgP+EREWuD8Gk31WXzXoWvw5/5xPioifXG+Q/WA330WSclXz8tFJPdgF+nWFZzm+4PAFyKSidNUL+hAWOg5BFV9V0Qa4DTtjxeRV4H/FFOFPb7rLmC6mojUBSaKSHc35sbAxwWsKx443T0Q5kp3l99ZyOfHuuv03dYpQFxuDKq6u5g6FCfd57UH53uUG++bIpLlTtfCPSADw0XkGpz9UxMnOfrGV5S8cwgi8iZO99WfPp9Z2DZqlG/d+U/6586LdevwrYjkzosB5hf2d1HCv5d4Ct8PuXHm8t2OYckSQvBsB57Kf6WK2xf8PM5VLD+LSGvgj9z57vIfuAfBF3EOfi8V9iEichKwP/fXmar+JiJDcPpqGwJpQH2ftzTBuaoHnD+ACHc9dX2WuQane6uHqiaLyECc/vXC6vmqqt6df4aqpuB0Uw11E8b7IjLPbTGVmKq+DLzsnkiejZNQS3TitgjDgWNxTsZniMgYCv51uB3noHRpKda9E2fbxgLJblkshx4MA2E7cImqrvItFJEeON+lU1V1vYj0wW3BHIYHgJ9E5HmfX9YFbiMR2Y3zgyJXC2BdAevciZOsTy0oWRbyd/FAIeXzfd6ahLPtc1XUfghJ1mUUPB8Ag0QkWkQiRGSU+4u7Pk6T+ze3a2EIOCdPReQWEXkAQFXTcBJFBE6/eWS+g3au84A33INlbv/u1ThdCruARTgnVxGR44DTfN67Baf/FeBK/vnF2BBY7yaDeJyTfjFF1PNSdzlEpL+I3OvWe4GIHOkutwrn3IWnBNsuj4iMFpFb3MkdwEb+2SYxbmvmcDTE2UYZblK+kH/qmIVz0hzgc+BMdxlE5FQRSSxqxaqaBcwDbnbf0wqnn/rLw4y1ND4A7nS/c1EikiginXDquwPY4J6oHwTUKarrqzBuy2AWzklkKHob/Qhc7JZ3wOnHL2id2cAnOBcMICK1ReRlEWlW2N9FEX8vvj7E6baMFJHaOH8bBbUEqwRrIVSMb0TE90TtQJwTk0cDq3H2wxKcqx324nzxf8f5VfR/OOcUPsc58TpDRP7AOXgqzkElHac/db2IXKSqvv2wT7rr/9r94/YAv+Ac4MC5euQtEVnvlvv+MfwXeNa92uNDnNZENZyrkq4Rkb9wfs3dC3woImPJd7+Aqi53f11/7R6ck3GuMMkSkReB+W5c2TjN/MwCtl+1fN0NuUbinAidLiIjcBLWDzjnY2rgNP83i3PVVmlNBea423oVTktmjogMxtkXd4vIIlXtLs5VPXNEpAbO/vt3CdZ/B/Ci27o6CAxS1U1uC7E4l4vIGT7TK3C6AUtiFM4VbGtx9uUnON/B34DBOPtzG84Pkc7AGzjborQeBlREJqrqyiK20VjgbRFZh/Md/gj/ripfdwDT3O9CNZyW52YReZ+C/y4iCin3TTqTcVqCa93PnQWU6kq3cBJhz0MwvtyD9GZVfSjYsZiqQUQi3BP8iMg7wHeqOjnIYVVJ1mVkjAkat8X1gdtl0wTnUthFQQ6ryrIuI2NMMM3Auf/jD5zuzCdU9cegRlSFWZeRMcYYwLqMjDHGuCpNl1Fy8p6wa8o0bFib1NR9wQ6jXIVjnSA862V1qjzKUq+4uLolvnTYWghBFBUVfjc9hmOdIDzrZXWqPCqqXpYQjDHGAJYQjDHGuCwhGGOMASwhGGOMcVlCMMYYA1hCMMaYkJWdnc24cWPp06cP48aNJTu7yIcZllmluQ/BGGOqkuzsbDp2bEtSUhIA8+fP55VXXmLFirVERQXm0G0tBGOMCUGJieNolpSEF2cMfICkpCQSE8cF7DOthVBGW7ZsZtKkCaSkOE/ha9v2eAYPHkatWrX8llu8+Ae2bdvKJZdcXuB6Zs6cQadOnWnfvkOB8wuyfPky3n33bcaMedKv/KyzunLiiSeRnZ1NTk4OV111Lb17n1vIWmDv3gx++WUNp57arcSfbYwJIK+Xx8Y/nneAnoTz4AaAJUsCNxhslUoIHo+Hr776glWrVtKhw0n06nUOkZGH30jyeDyMGnUvQ4YMp0uXUwB4442ZjB37MA8/PNZvuW7dTitsNQBcf/1Nhx1HfjExMUyZMg2A1NQU7r13OI0bx9KxY+cCl1ddy9Kliy0hGBMCqs/7lPrXD/Ar830UYteu3QP22VUmIXg8Hm666Ro+++yTvLJ+/c5nxow3DjspLF26mBYtjs5LBgBXX30dV199Gbt27eT556cQFRVFWloaZ555FuvW/cWQIcOYOHE8a9asol27tqj+wejRjzBjxov07NmbGjVq8P77s4mIiGT9+nX07n0ON988iGXLlvLCC88RHR1N3bp1efjhx0sUY8OGjRg8eChvvDGTjh0789Zbr/H11/PxeDx07346t9xyG0899ST79u2lefMWtG/fgfHjHyMqKprIyEgeeeRx6tWrX/wHGWPKJjubuKaN/IoeiIlhTEZG3nR8fDzDhycELIQqcw7hq6++8EsGAJ999glfffXFYa9z48a/ad26jV9ZREQErVody4YNfwNQv34Dxo4dnzd//fp1rF69kmnTZjBw4EDWrv3VLyFFRkby66+/cP/9DzF16svMmfM2ABkZe/jvfx9kypRpxMTULVWz8bjj2rBhw3oAvF6YPPk5pk2bwaeffszevRlcc8319Op1Dv37X0pqagqDBw/n6aen0qFDRz7//NPD3TzGmBKq+eLzhySD5K0p3PX7RhISRtKnTx8SEkYG9IQyVKEWwqpVKwssX716FX369D2sdWZn55CTk3NIudfrISLCGWDw+OPb+c37++91tGvXnsjISI499lji45sc8n6RttSsWTNvXQB169ZjwoTH8Xg8bN26hU6dulC7du0SxRkREYnH4wwWGxUVxfDhQ4iMjCQtLZXdu3f7LVu/fn2ee+5pDh48yM6dyZxzTr8SfYYxpvQiMvYQ2+oov7L0GW9w8Px/Ac4BOiFhJHFxdUlO3hPweKpMQujQ4aQCy088seQncfM79tjWvP/+bL8yj8fD+vXraNmyFQBRUdF+871eiPAZjDYi4tCRaatVO3Rkw8ce+x/jxk2iVatjGT9+7CHzi7JmzSpatz6OLVs2M3v2W0yf/jq1a9fh2msPPcE9adIErr32Brp3P4PXXpvBgQMHSvVZxpiSibnvP9R6+QW/suQd6f4HiApWZbqMevU6h379zvcr69fvfHr1Ouew19mlyyls3bqFH374Pq/srbdep1279jRs2LDA9zRr1ozff1e8Xi9//vknSUk7SvRZmZmZHHlkU9LT01i+fBlZWVklel96ehrPP/80V1xxFXv27KFhw0bUrl2HNWtWk5SURFZWFhEREXktnd270znqqOYcOHCAH374PuA3whhT1URu20pcfD2/ZJD62VckJ+0OajKAKtRCiIyMZMaMN/jqqy9YvXoVJ57YocxXGUVFRfHkkxOZOHEc06e/gNfr5fjjT2DEiPsLfU+bNm1p3rwFgwbdSJcunTjmmFYFtgjyu+yyK7nrrltp3vxobrzxVqZPf4Fbb729wGUzMjIYMuQ2vF4vBw4c4Nprb6Rz55PJycmhTp0Y7rprIO3bn8ill17BxInjGDJkGM8/P4UmTZpwxRVXc//9CTRtehRXX30dEyeOp3fvc2nd+rjD3k7GVFX5r2y8bMokqvv8gMw68STS5n8XxAj9VZpnKofLE9MOHjzI/Pmfc955/6JOnWr07duPt9+eG9ATRRWpovo6K1o41svqFFgej4cbb7yaefM+pQOQ/yzmrh9X4Tm6ZYnWVZZ6leaJaeFxFKpEqlevjupvzJ49i+joagwadGfYJANjzD++/PJz5s37lBTAtwN53Rk9qPvuR8EKq0h2JAqCYcOc64hD6deMMaZ8/T31GfJ3azQCzmnalCnBCKgELCEYY0x5cm8w8z2T+AAwxn0dyr30VeYqI2OMCbSar04/5AazmvyTDAD697+kQmMqjYC2EETkMeBsIBp4QlXf8ZnXHZiAs73eVdUxBa/FGGNCXEYGca2a+hWlT5nKgA/f58C8f+7279v3vMO+EbYiBCwhiEgPoKOqdheRRsAq4B2fRV4BegNbgEUi8qaq/hWoeIwxJhDqPDya2lMm+pUlb02BqCheuXxAuV7qHmiBjOwH4Er3dRpQXUQiAUSkFZCiqptU1QN8BBQ+PnOI2rZtK+ec04MhQ27L+zdp0gT++EN56aWpAHz//YK8m8i+/vrLEq97zpxZeevItXz5MkaNuveQZR98cCQHDuwvQ02MMaUVuWM7cfH1/JJB2rsfOTeYuVcORkZG0qdPX4YPT6BPn74hnQwggC0EVc0GcofpGwh84h78AY4Ekn0WTwL821v5NGxYm6io4m/gqkgHDtShVatjmDXrzUPmnXbayQC8++4szj33bOrUqcNbb83kyiv9+w/j4uoe8l6AmJiaZGXV8JvfoEFtatSIPuQ9zz4bWtcsFFanyi4c62V1OkwXXQQffvjPdIsWsH49DQJ4wK+IegX8KiMR6Q8MAnzHiDiYb7EIOOQKLT+pqfvKObKyS0nZS3a255BLR3MfXHPGGWexYsUKbrrpFk477QzWrl3LoEF38Nhj45g69RnWrl1DZuYBrrjiKnr3Ppdly5YyefIEGjeOpVGjxjRtepTfutPS9nHgQNYhn3f55Rfy6quzeP31V8jM3MeGDRvYsmUTw4Yl0K3baSxY8BVvvfU6kZGRtGvXnsGDhwZsm4TrpbThWC+rU+lVW/sbjXp09StLWbCYnOPbwa69AfvcMt6YVuJlA31SuS8wGjhXVdN8Zm0D4n2mjwC2luWz6jw0ihofvl+WVRziwIUXs/ehwz/X3a/fBbz44vOMHz+Z2rVr8+abM3nssXGsXLmCpKTtvPbaa2zevJNbb72eM844i6lTp/DQQ4/SqlVrEhKG0rTpUcV/iI/IyEh27NjOhAmTWbRoIXPnvstJJ3Xi1Ven8/zzLxMdHc2oUSNYs2ZVqZ7MZkxVlZ2dTWLiOBYvWsj7q1fSKD09b96B3uew+43ZQR9/qDwF8qRyfeApoJeq7vKdp6qbRSRaRFrgnFT+FxC612IVYePGDQwZclve9CmnOI+vLMratb+wZs1qrr/+eg4ezMbr9bBr1062b99Oq1atAejYsTP795f+vMBJJ3UCoEmTJmRk7GHTpg1s376V4cMHA87jMrdt22oJwZhi5D7k/oSkJL7NN2/XyrV4jiyyl7tSCmQLYQDOHduzRCS37Ctgtaq+BwwD5uJ0Fb2mqpvK8mF7HxpTpl/zh6tFi6PzHleZa/nyZcW+7/zzL+See4YV2gzMfX5BafkOlJc7TtVxxwkTJz57WOszpqqaNP5xdiQl+ZVNBLYmjCQhDJMBBPak8jRgWhHzvwU6BerzQ0VERGTe0NK5B/l27dozZcpEcnL+zYEDB3juuckMG5ZAXFw8Gzf+TYsWLfn555844YT2Zf78Fi1asnHjBlJTU2jYsBEvvTSV/v0vJTY2rszrNiZc1Zj1Bo8/9aRfWQMgHegRwIfcB5sNXRFgnTp1di9HfY42bYTbb7+ZqVOn07nzyVx11VVkZWVz8cXOg2puu+0u7r//Xo48silNmjShoJFoV6xY7tdFNWrUw0V+fs2aNRk69D/cc89QoqKiEGlL48ax5VtJY8JFATeY/R/wtM90IB9yH2w2/HUQ2VUelUc41svq5K/2+Mep8+RjfmUt4+LYkPzPFfLx8fEBf65xQWz4a2OMqQARO3YQe6L/A6DSX5rJwQv7s8i9ymjJkkV07dqd4cMTwnq4+vCtmTHGFKPunQOpOeftvGlv9ersXLcVqlcHnKciJiSMDFZ4Fc4SgjGmyqn25x80Oq2LX1nqh5+T3bVbkCIKDZYQjDFVh9dLg3N7Er3y57yi7BNOJPXLb6EEzzYPd5YQjDFhK+9O48U/MKBZCwa/OdNvfsp3S8mRtkGKLvRYQjDGhKXcO413JSWxF6jhM2//RZew54UZYTXsRHkI7bFYjTHmMCUmjuPMpCSy8U8Gj94+mD0vvmLJoACWEIwx4WfvXp4cN5bZPkVjcYZV/vy3X4IUVOizhGCMCSu1np9C3DFH+pXFAv91X4fzncZlZecQjDFhIWLnToivR4xP2YiYGJ7MyMibjo+PZ/jwhIoPrpKwhGCMqfRiRt5DrZf8x9Lc+ddmhtWqjbcK3WlcVrZljDGVVuS6v2jczX/Q5PTpr3PwggsB5wBXle40LitLCMaYysfrpf6VF1N9wdd5RZ6GDYncsoWDGdlBDKxys5PKxphKJernn4hrUt8vGaR+MI9dugFq1QpiZJWftRCMMSHL907j7qd247E3X6Pa1i1587NOPpW09z/JG4zOlI0lBGNMSMq90zgpKYnzgSe+W+A3P2XBYnKObxec4MKUJQRjTEhKTBzH7qQk8j8Za/nxJ9D8q+9tMLoAsHMIxpiQ1Pj9OWTmKzsGGB4XZ8kgQCwhGGNCSkRaKnHx9Rjxx+95ZY/iDDvxN3ancSBZQjDGhIzaj48hts3RfmVNgVHua7vTOLDsHIIxJugiN2+icecT/Mr2PD6BjOtv4rqJ4+1O4wpiW9YYE1R1B99GzXfe8ivbuXY93kaN7U7jCmYJwRgTFNV+/YVGPf3PB6S/NJODF/YPUkTGEoIxpmJ5vTQ4rxfRy3/KK8qJb0Lq4uV4Y+oGMTBjJ5WNMRUmeuF3xDWp75cM0t77mJQ1f1gyCAHWQjDGBEzu0BM//vA9H//0I7X378+bl9W1O2nvzIWaNYMYofFlCcEYExC5Q0+clZTE1/nmpX75LdkdOgYlLlM4SwjGmIB45olH2ZGU5Ff2GrDm7nv5jyWDkGQJwRhT7mpNe5Yxkyb4lbUFFOixbGlQYjLFs4RgjCk3EUlJxLZv7Vc2Frgf8gaps6EnQpddZWSMKRd1Hrz/kGTQpXFj/ss/ycCGnght1kIwxpRJQc81zhjzOJmD7uSjnBwS7SH3lYbtGWPM4fF6qXfzddT45EO/4p2r/8DbpAkAUVFRNvREJWJdRsaYUotasZy4JvX9ksHu518iOWl3XjIwlY+1EIwxJefx0LD3mUT9sjqvKCe+CanfL8XboGEQAzPlwRKCMaZQvg+5v6HJEQycPctvftpbc8jqdU6QojPlzRKCMaZAuXcapyUlkQLU8ZmX1eUUZ9iJmJhghWcCIKAJQUTaA3OBRFWdkm/e38AmIMctulZVtwQyHmNMySUmjuPcpCRm5it/5prruXLiM0GJyQRWwBKCiNQBngbmF7HYeaqaEagYjDGHJ2J3Ok+OG+tXNgu4Aei2eRNXBiUqE2iBvMroAHA+sDWAn2GMKWe1pkwitnVzv7KTgKuAg9idxuEswuv1Fr9UGYjIQ8DOQrqMFgHNgYXAfapaaDDZ2TneqKhqgQvUmKpu2zZo2tSvaFrt2ty1b19ev+4RRxzBpk2b7OayyiWipAsGc6+OBr4AkoA5wBXA24UtnJq6r4LCqjhxcXVJTt4T7DDKVTjWCcKzXr51qjNqBLWnPec3P2XRT1x49DH8me9O49TUzGCEWyLhuJ+gbPWKiyv5g4eClhBU9dXc1yIyDzghWLEYU1VV++sPGnXv4leWMep/ZA4ZCpGR9pD7KiYoCUFE6gIf4ZxU3gecAbwfjFiMqZK8Xujfn0YffOBXvOvnX/Ec1SxIQZlgC+RVRl2ACUBLIEtELgc+ANar6nsi8jbwvYjsA34GZgcqFmPMP6KWL6Nhv15+ZXsmPsP+q6+DiBJ3N5swFLCEoKo/AT2LmP8MYBczGxNAHo+Hr776glWrVnJS+/Zc9tgjRP/6yz/z6zcg5Yef8MbFBTFKEyrsUgFjwpTH4+HGG69m3rxP6YfzoBo/777LrjP6BCEyE6osIRgTpr788nO+nvcpuwHf60xSjm6J58tvnXsNwvCKHHP4bPhrY8JU6qQJHMA/GfQC7uzaDW/9BkGKyoQySwjGhJmI3enExddjyI9L8so+xBmc7mucC4yMKYglBGPCSEHDTnQFLgJyb+3s3/+Sig7LVBJ2DsGYMBCxYwexJx7nV7bv+psYsH0rS7/4PK+sb9/z6NOnb0WHZyoJSwjGVHJ1HriP2lOf9StL+WYROe1O4EX3stPVq1dx4okd6NXrHCIjrWPAFMwSgjGVVLV1f9KoW2e/sr1338u+hJFQzRkIMjIykj59+lqrwJSIJQRjKhuvl3oDb6TGh/6jvexa/DOeVscGKSgTDqztaEwlErViOXFN6vslgz2PPUny9jRLBqbMStRCEJHewBCgET5ja6tqjwDFZYzhn4fcL1m0kNf+/J3m27fnzfNGRZHy0xo8RzYtYg3GlFxJu4yeBcYAGwMYizHGR+5D7jskJbEg37zdz77AgcuutMHoTLkqaULYoKr5n7VtjAmgyeMfZ01SEr7Dzv0JvDZ4KIMvHxCssEwYK2lC+FhE7sC50TErt1BV1wUkKmOquBpz3mbsU0/6lf0L+BjosXolg4MSlQl3JU0I/wFygBE+ZV6gVblHZEwVFpGxh9hWR/mVfQ1cDOx2p+0h9yZQSpQQVLVFoAMxpqqrNe1ZYkbd51d2cYMGzE1Ly5uOj49n+PCEig7NVBElvcqoKfAocDLgARYDo1Q1OYCxGVMlROzcSWw7/8b2/ksvZ89TU3iuenXa5HvIfVSU3T5kAqOk36wXgU+BcTiXnfYDXsIZM8sYc5jqjHmI2pOf8itL/fwbsjs6dyDbQ+5NRSppQqilqk/7TP8iIhcGIiBjqoLIDX/T+JQOfmX7bh/M3gf+B9WrBykqU9WV9E7lWiKSd6bLfV0zMCEZE97qDr7tkGSQ8t1S9j4y1pKBCaqSthAeAZaJyFacLqN4YGDAojImTOTeabx48Q9c0upYhr/yst/8jFEPkTl4aN5gdMYEU0mvMvpYRFoBbXASgqpqZkAjM6aSy73TOCkpia+Bnt/7GRjwAAAViklEQVT532+8a+lKPC2PCU5wxhSgyC4jEbnZ/f9hYCRwGXApMNItM8YUIjFxHJKUhBfo6VM+59x+JO9It2RgQk5x5xA87v85BfzLDmBcxlRu2dkMfTqRb32KdgFNgcn799sYRCYkFdllpKqvuC/TVXWi7zwReTBgURlTiVX/5CPq33SNX9mVwDvu6+vsTmMToopMCCJyNtALuE5EGvnMigZuAv4XuNCMqWQyM4lt3YyIrLzhvlgdFUXP7GxS3Gm709iEsuK6jNYCv7mvfbuL9gFXBzAuYyqVmm/MJO7oJn7JIP3N2cRuTOLWhJH06NGThISRrFix1u40NiGruC6jbcAbIvIDUF1VfwcQkc6qurwiAjQmlEWkpxF7nP9QXwf6nMueadPxxtS1O41NpVLSG9PuBHxH3RohIk8EIB5jKo1akxMPSQZpH3zG7jdm442pG6SojDl8JW279lDVvDNhqjrAbTUYU+VE7thO4xPb+JVlXnM9GY9PgJp2A7+pvEraQvCKSN499SJSB7BbK02VU+eB+w5JBqlffkvGxGcsGZhKr6QthBeAX0XkZ5wkcjJgl52asJc79MSO77/hlUWL/Obt+7+72TvifoiODlJ0xpSvkg5dMV1EvgC64AxdMVRVNwc0MmOCLDs7m44nCZOSk7kq37yUhcvIOa5Nge8zprIqUZeRiNQErsA5l/A+cLSIxAQ0MmOC7O37/sOOfMngHmDEPfdZMjBhqTRdRtuA093p9sAwnCRhTHjxeGhwUT+GLl3sV9wS2AD0yFduTLgo6UnlZqp6L5AJoKpTgSMCFpUxQRK98DvijmhAtM9B/xacftIN7rQ95N6Eq5K2EHIThxfyrjKySypM+MjKolG3TlTbtDGvKOeIIzkp6yC/7NqVV2ZDT5hwVtIWwvvuSeVWIjIZ+Bl4LXBhGVNxqn/8IXFHNfZLBukvzSRllfLF6j9ISBhJnz59bOgJE/ZKepVRooh8h3MOwQtcrao/BTQyYwItM5PYY44kwuPJK8rq1Jn0We/hbdAQgKioKBISRhIXV5fk5D3BitSYClHcaKe98hWtdv+vLyK9VPWrYt7fHpgLJKrqlHzzugMTcLqe3lXVMaWK3JgyqPn6q9QdPsSvLG3We2Sd3TtIERkTfMW1EB4oYp4XKDQhuOcZngbmF7LIK0BvYAuwSETeVNW/ionHmDKJSEslts3RfmUH+p3P7mdfhBi7ktpUbcUlhFdUdYaIDFLVF0q57gPA+cCI/DPc5zOnqOomd/oj4FzguVJ+hjFF8n3I/X0REZz/7Td+81M/+oLsU7sGJzhjQkxxCeFBdwyj4SKSk3+mqr5c2BtVNRvIFpGCZh8JJPtMJ+E8XbBQDRvWJioq/IZPiosLv1ExQ6VO2dnZNG/enIjt29maf+agQTB5Mg1LMf5QqNSrPFmdKo+KqFdxCeEe4DygAXBmAfMLTQjFOJhvOgL3ktbCpKbuO8yPCl3heKIylOo0btxYErZv5+585Yk33MJ1j06APVnOvxIIpXqVF6tT5VGWepUmkRT3gJw5wBwRucx9XV62AfE+00fAoT/ijDlckev+4slxY/3KHsUZkfH0v9dxXVCiMia0FXeV0UhVHQtcLCL9889X1RsO50NVdbOIRItIC5yTyv8CLjmcdRmTX907bqHmu7P9ytoAf7iv7U5jYwpWXJdR7mMyfwBicLp6UnG6eIokIl1wLittCWSJyOXAB8B6VX0PZyykuThdRa/lnmA25nBVW72KRr3P8CsbHRPDmIyMvP5Iu9PYmMIVlxAWiMj7QGdgBeABugHzgFuLeqN741rPIuZ/C3QqTbDGFMjrpf7F51N90UK/4l3Lf+HOI45kf+I4lixZRNeu3Rk+PMHuNDamEMVeZQRsBi53rxrCveooEXgIGBXQ6IwpRvTiH2hwUT+/st2Tn+PAgGsgIsIecm9MKRSXEM4FuuYmAwBVPSgiw4BFWEIwwZKdTcMzTiFq3T/3MubENyH16x/wxsUFMTBjKq/iBrfb55sMcqlqFpARmJCMKVr1eZ8S17SRXzJIn/EGKWv+sGRgTBkU10Io6iJtTxHzjCkXvncan3HyqTz63NNE7N+fNz+rYyfS35mLt36DIEZpTHgoLiGcJiIbCyiPAGIDEI8xebKzs+nYsS1JSUlcDzz23QK/+WmzPyCrR8+gxGZMOCouIRQ47oQxFSExcRyZSUmH3ML+67GtifvyO6hTJyhxGROuirtTeUNR840JpOZz3mZ3vrJuQI2jmjHbkoEx5a6kT0wzpsJE7NxJXHw9/u1z0vgFoAawBLvT2JhAsYRgQkrtxx8htl0rv7ITgdtwbpO3O42NCRy7ZdOEhMgtm2ncqZ1f2d5h97D77ns57+lEGtmdxsYEnP1lmaCLuWcYtV71H0l91+Kf8bQ61u40NqYCWUIwQVPtj99pdPrJfmUZDz9G5m13QaT1ZhpT0SwhmIrn9VLvxmuo8dnHfsW7VvyGp+lRQQrKGGM/w0yFilr5M3FN6vslg92TnyN5R7olA2OCzFoIJqByh55Ysmghr/35O823b8+b52nYkJSFP+GNtZvejQkFlhBMwOQOPdEuKYkF+ealv/ImB8+7IChxGWMKZl1GJmAmTXiCZUlJfO1TthIY/X/DLRkYE4IsIZiAqP7ZJzw+4Qma+5T1BDoCC1b8HJygjDFFsi4jU7727ycuvp5f0WfApUCmO21DTxgTmqyFYMpNjVlvQK1afmV9GzbkPP5JBjb0hDGhy1oIpswi9uwm9thmfmX7B1zDngmTmREZ6VxlZENPGBPy7C/TlEmtac8SM+o+v7KUb5eQ0/Z4ABt6wphKxBKCOSwRu3YRe/wxfmX7Bg+l9qQJ5KTsC1JUxpiysIRgSq32E49SZ8ITfmW7lq7E0/IYalerFqSojDFlZQnBFMn3Ifd9T2jP/c8/4zc/Y/QjZA7+P4iICFKExpjyYgnBFMr3IffPAnfme8j9rlWK54gjgxOcMabc2WWnplCJieOo7z7k/k6f8rf7XeAMRmfJwJiwYgnBFOpfL0/jd5/pbKAx8My+vdZFZEwYsoRgDlFtzWri4uvRY9euvLKLgWggBbvT2JhwZecQzD+8XhpccA7Ry5bmFW2MjKS9x8Med9ruNDYmfFkLwQAQtWQxcU3q+yWDtPc+JnrzTu5IGEmPHj1JSBjJihVr7U5jY8KU/WVXdTk5NOrakWobN+QVZZ18KmlzPoRatexOY2OqEEsIVVj0V1/S4KpL/cpSP/+G7I6dgxSRMSaYLCFURQcPEtvyCCKys/OKDlxwEbunTYfo6CAGZowJJksIVUz1D96j/sAb/cpSvv+RnDYSpIiMMaHCEkIY8x124swup/DoxPF+8zNvHkjG2PEQadcWGGMsIYQt32EnbgUezT/sxLLVeFocHZzgjDEhyRJCmEpMHMd+d9gJX5+fdgad3vvY7jQ2xhzC+grC1DHvvEV6vrIjgEejoiwZGGMKFNAWgog8DPQGagK3q+oyn3l/A5uAHLfoWlXdEsh4qoKInTuJbdfKbzC6u4Dn3Nc32LATxphCBCwhiMjZwCmqerqItAeeBXrkW+w8Vc0IVAxVTZ2HR1N7ykS/skZAqvvahp0wxhQlkF1GZwNzAVR1DdBURGoH8POqrMgtm4mLr+eXDHZPm862rSkMtGEnjDElFOH15j/tWD5E5AXgM1Wd404vAq5R1fXu9N/AIqA5sBC4T1ULDSY7O8cbFWWPZzzEzTfDjBn/TEdHQ0oKxMQELSRjTEgp8UnDQP5cPJhvOgL8LnoZDXwBJAFzgCuAtwtbWWpq+D24PS6uLsnJe4pfsACR6/6icbdOfmVps94j6+zekOmFzMNbb1mVpU6hLBzrZXWqPMpSr7i4uiVeNpAJYRsQ7zMdB+zInVDVV3Nfi8g84IQAxhJW6l19GTXmf5E3ndOiJSnfL4WaNYMYlTGmsgtkQvgUeBR4TkQ6A+tUNRNAROoCH+GcVN4HnAG8H8BYKiXfO427dTuNhH4XENf7DL9lUj+dT3aXU4IUoTEmnAQsIajqTyKyUkSW4zx98VYRuQlIV9X3RORt4HsR2Qf8DMwOVCyVke+dxgBPf7eAuHFj8+Znde1O2nsfg50kNsaUk4AeTVR1BDDCp2i1z7xngGcC+fmVWWLiOJKSkugKLM43L+W7peRI22CEZYwJY3ancohasmgh2/FPBm8DPc88y5KBMSYgrL8hBEUv+JoF33/rV3Y0sBFI6HZaUGIyxoQ/SwihJDubuKaN/IomA0Pd13ansTEmkKzLKERU/+iDQ5LBjpVr2Wx3GhtjKogdXYJt/37iWsT7Fe0d+QD7hicQiT3g3hhTcSwhBNNLLxE3cKBf0c4/NuKt3yBIARljqjJLCMGQkUFcq6Z+RXvGTWT/jbcEKSBjjLFzCBWu1pRJhySD5PXbLBkYY4LOWggB5Dv0xNkdO/G/p/2fVcA775B8Vt/gBGeMMflYQggQ36EnHgdG5HvIffKmZOKaxUIYjsxojKmcrMsoQBITx+F1H3LvO3bHtCuvJjlpN9SoEazQjDGmQJYQAuSsV6eT5DOdhNMce337tiBFZIwxRbOEUM5yH2d5wY7teWUnA02AHKCrPeTeGBOiLCGUo3rXXkHjTu3ypldFRREJ/ORO29ATxphQZieVy0G1v/6gUfcufmUpC5cRd0wr7kkcx5Ili+jatTvDhyfY0BPGmJBlR6cyatCnB9GrVuRNH+hzLrtffwciIojChp4wxlQelhAOU9TqlTTsfaZf2a6ff8VzVLMgRWSMMWVj5xBKy+ulUfvj/JJB5jXXk5y025KBMaZSsxZCEfI/5P7e084g9pIL/JbZ+dt6vI0bBylCY4wpP5YQCpH/Ifff5rvTeN/Q/7D3/geDEZoxxgSEJYRC5D7k/nzg43zzktdthZiYYIRljDEBYwmhEEsXLcSbr2w4sLxHT2ZbMjDGhCFLCAWoMedtvsn3kPsawEEgwe40NsaEKUsIvrKyiDvK/wTxDcBM97XdaWyMCWd22amr5kvTDkkG2zYlc4Q95N4YU0XY0S0zk7ijm/gVpc+cxcG+59mdxsaYKqVKtxBqT3jikGSQvCOdg33PC1JExhgTPFWyhRCxZzexx/rfVZw291Oyup8epIiMMSb4qlwLoc5/E/ySgSc2luQd6ZYMjDFVXli3EHyHnujdoSMPPjPJb37K/O/JObFDkKIzxpjQErYJwXfoiTeBq3yGnsg+4URSv14YvOCMMSYEhW2XUe7QE2OAq3zKnxx4uyUDY4wpQNgmhMWLfwDgZ3f6EyAC+Ph3DVZIxhgT0sI2IXTrdhoAc3ASQe6g1faQe2OMKVjYJoThwxOIj4/3K7OhJ4wxpnBhe1I5KiqKFSvWkmgPuTfGmBIJ66NjVFSUDT1hjDElFLZdRsYYY0rHEoIxxhjAEoIxxhiXJQRjjDGAJQRjjDGuCK83/6PkjTHGVEXWQjDGGANYQjDGGOOyhGCMMQawhGCMMcZlCcEYYwxgCcEYY4zLEoIxxhggzEc7DSQR6Qm8A/ziFq0GHgFeBRoAm4FrVfWAiFwC3AvUBJ5W1ZdFpBrwLNAe5xk+16rqehFpA7wA1AaWAXepqldE7gSud8v/q6qflHN92gNzgURVnSIi8RVZFxGpA0wHmgF7gQGqmhKAes0AugC73EXGqerHlaleIvIYcDYQDTwBLKCS76sC6nQBlXg/iUhtYAbQBKgDPAwsJsT3k7UQymaBqvZ0//0bGAdMV9VuwN/AtSJSFxgP9ANOB+4VkRjgBsCjqqcDjwH/c9f5AnCvqp4CxAFni8ixwO3AWUBfYLyIRJRXJdwvztPAfJ/iiq7LvcBPqnoazgF8WIDqBTDSZ799XJnqJSI9gI6q2h04F0ikku+rQuoElXg/ARcBy1T1LOAyN+6Q30+WEMpXT+AD9/VcnJ1zCvCjqqar6j5gIXAmzq+hue6ynwE9RaQ60FpVl+Rbx1nAZ6qapao7gG1A23KM+wBwPrA1iHXxXUfusoGoV0EqU71+AK50X6cB1YFeVO59VVCdCjo2VZo6qepbqvqkO9kMp0XQkxDfT9ZlVDbtRORToC5OBq+rqpnuvCTgCOBIINnnPYeUq2q220SMA1IKWHZ3Iev4rTwqoarZQLaI+BZXdF18151bFoh6AfxbRBKA7cDgylQvt04Z7uRA4BPgosq8rwqpE1Ti/ZRLRJa46zof+C7U95O1EA7fH8AYnL7O64EXcfr6ckUAXuBgvvcVVg6QXcJlc8sDyfczK6Iu+csDVb+ZwP1uU/5HnL7dSlcvEekPDMLpBgiLfZWvTmGxn1S1K3AJ8Fa+mEJyP1lCOEyqukVV31BVj6qux/kVE+OeTAInG2/Fab7F+7z1kHK3KZiFk80bFLVsvvJA2lPBdfEtP5IA1U9V56vqMnfyY+CEksQUSvUSkb7AaKCfqqYRBvsqf50q+34SkZNFpAWAqi7HOdbuDfX9ZAnhMInIVSLykPs6FudqgheB/u4il+J8kZcCHUSkvnuyqCvwHfApcLG77L+AL1TVA6wQke5u+SXuOr4A+opItIg0BRqp6u8BruJnFVwX33Xkfl65E5G3RaSDO3kGsKYy1UtE6gNPAeerau4VOJV6XxVUp8q+n4DTcE/iikgTnG7ljwjx/WTnEA7fR8AAEVmIk1jvAn4G3hSRuwEFZrn9f6OBbwEP8LCqZorI+8BFIrIM2Adc4653BDBdRKKAb1R1IYCIvIzTdPZQDlfg+BKRLsAEoCWQJSKXA9cCr1dgXaa6224ZTn/ngADVazTwoojsA/YAt6jqwUpUrwFAQ2CWz7mRG4FXKvG+KqhOlX0/TXU/+zugBs7x4Scq9vhQ6jrZ8xCMMcYA1mVkjDHGZQnBGGMMYAnBGGOMyxKCMcYYwBKCMcYYlyUEYwogIi1FZHM5r7OPiHxTnus0pjxZQjDGGAPYjWnGlIqI3ALcCWTi3Oxzq6qmi8ggnIHZ9rv/BqhqmohcDDyKM2yABilsY0rEWgjGlJA7Ns3DQG9V7QFsBO52Z9cC+ruDsW0ArnPLpwCXqeo5QE4Fh2xMqVhCMKbkOuOMXb/bnZ4PnOy+3g+8IyILcMadjxWRxkAtVV3rs7wxIcsSgjEll/8pdRGAR0SOwXmq1ZVuC+G9Qpa3vzcT0uwLakzJ/QicLCL13Ol+OM/JbQgkq+o2d+Tbc3AGNNuFM6he7oht5fEUOGMCxga3M6YAItIS5yTwIp/ipThPqbsT5/Gcm3Ae6JIJvIsz5vxfwDs45w4G4Dw+cTSwDlgLnOqefzAm5FhCMMYYA1iXkTHGGJclBGOMMYAlBGOMMS5LCMYYYwBLCMYYY1yWEIwxxgCWEIwxxrj+Hy06Bosty3VLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x6adee730>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEVCAYAAAA2IkhQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucXHV9//HXkoAxZIVINrdGiIK+NaGo1BAiBRJRgiQGo7W0IoqCv1SKrfL7GbXaqJRSEbwVqiLholIB8aeGGqFFVEAaF4LghehHLRcJiWEhkcQkXMJu/zhnwjjM7s7O7JkzZ+b9fDzyyJzrfD9zZudzvpdzTtfAwABmZmZZ2iPvApiZWftzsjEzs8w52ZiZWeacbMzMLHNONmZmljknGzMzy9zYvAtgNlKSBoD/AXals8YCNwN/FxF/qHOfvwSOjohNFfPfApwWEfPr3O98YGVEHFTP9q1E0kzgNxGRye9Go5+1tTbXbKyo5kfEiyPixcBs4LnAP9a7s3Rfm4Zfs7VJGpN3Gcyqcc3GCi8iHpd0PbAYQNJewPnAQuAp4Erg7IgYkHQG8LdAF7AVeHtE3J3Wlp4HbAD+FVgC/A64qfQ+kn5AUku5onJa0mHAvwHPIalx/X1EfLe8nJIOBi4G9gH2BD4bERdWrHM78C8R8Y10einw/og4XNIS4GxgL5Ka3akR8TtJH03L/lLgSklXAl8BpgHPAq6KiA+l682IiNPSfe+elnQ08Gng2WlRVkTENbUeA0nPBb4AvAzoBy6NiE+ky44HPp6W5Q9pue+StMdgn7W1H9dsrPAk7Qv8NfDf6ax3AwIOBuYAS4HjJXWT/FgfltaIPgEsqtjdccCxwEuAo4GjaizGF4HPRISAc0l+eCt9BPhCRMwC5gHHSHpWxTpfB04om3498DVJfwJcDrwpLfsPgQsqyv3aiPgk8F7g5vR9/hR4vqRpw5T/fOC9EfESkqS9dLiAK5wDbI6IF6WxnS7pyLSm9RXgXelnsyp9r1KZ6/msrYBcs7Gi+oGkXSRn+c8FPkuSPCD5sb4wIp4EnpT0ZZIfz/8EngROlXRVRHy9yn6PAq6LiO0Akr4GvK6G8swhOaMHuAV4QZV1NgJvlPRz4McRUe0H/RqgN/2R7iJJhh8mqaX9KCIiXe/zwEOSutLp3ojoK3ufv5B0I7AmIt6cxjJU+TcCb5X0UET8AnjzsBH/sUWkCSoitkj6D+DYiLhF0vSIeDxd7xbglPR1vZ+1FZBrNlZU89Mz/MNIfuS/GhGlAQOTgfMl/TLt+P97YEK6/FUkZ96/lHSLpNkV+30usKVsenON5TkZuF1SADeQJIpKy4GfAlcB6yX9TeUKEXEPsB54JcmPcUTEA2lMR5TF9CPgUWC/KuX8DPAN4HPAJkkfqaH8p5A0cd0g6deS3lDDNuUmA4+UTW8GetLX75X00/SzuZynf3fq/aytgFyzsUKLiIcl/StJrabU/PQ74FMRcW2V9X8GnChpT+BMkuauI8tW2ULSp1Iyuez1U/xxEumG3aO0vgDMjYg7JR0E/LrKez9GUkv5sKSXA/8l6bsR8ZuKVb9O0o8xDvhaWUw3RsQzkkBljSUi+kmaqs6XdCDwHUk/HKz86TabSZLy30t6FfAtSf9ZqnXU4CFgEnBvOj2JJNEdBfxfkqbLeyW9GliZrjPUZ21txjUbawefBF6ZdnIDXAu8U9KekrokfVjSQkmHSLpG0l5pE9tdPLMGsgY4TtKzJY0H3lS27EHgxQCSXgLMSufvQ1Ir+EWaxM5I15lQvmNJ3y6rSf2SZIBCNdcArwaOJ0k8AP8FHJkmMiQdJunT1TaW9EVJr0kn7yNJBF0V5d8beG36ek9JN5X16/yUZJBDP7X7D+Dt6f4mkSTL1cBEYBNwf/p5vBPYO23+G+qztjbjmo0VXkRsk/RxkjP5w4ALgQOAn5F8x3tJmpZ2APcD6yQ9AfyeZGRauf8g6X8Ikh/p1SRNb5AktX+XdATwC+C7wBiSH+fvAL8CHgb+jqQP57+Afyjb9wXAV9PRcgMk/UqVtRoi4lfpSK0NEfFgOm+DpFOB/58OKthOMhCimouAz0sqDSC4FrgRmACcLOk2kqa61cCkiHhS0krgxjQJ7ALOiIidVfY9Jm3GK/da4EPAFyX9iiRJfTwibpM0juQzvoekX+gM4FDgq8BbGPyztjbT5efZmJlZ1tyMZmZmmXOyMTOzzDnZmJlZ5pxszMwscx6NBvT1bWu7URITJ45ny5YdeRdjVDmm4mjHuNoxJmgsrp6e7moXL1flmk2bGju2/W7+65iKox3jaseYoHlxOdmYmVnmnGzMzCxzTjZmZpY5JxszM8uck42ZmWXOQ58tF73rNrF6zX1seHgH0yeNZ9G8mcydNSXvYplZRpxsrOl6123iomvv3j29vm/77mknHLP25GY0a7rVa+4bZP79TS2HmTWPk4013YaHq1+tvPGRWh8KaWZF42RjTTd90viq86ftt3eTS2JmzdL0PhtJZwHHkDxffVlErC1bNo/kaYjjgG9ExNmDbSNpMvBlYF+Spw6eFBGPl+3rSuDxiDilKYFZzRbNm/lHfTZPzz8gh9KYFVPRBtk0tWYjaQEwJyKOAN4GfKpilS8BJwKvAF4n6cAhtjkPuCwiDid5zvpJZe/zGuDALGOx+s2dNYVlS2Yzo2cCY/boYkbPBJYtmd3SfyhmraQ0yGZ933b6BwZ2D7LpXbcp76INqtk1mwXAKoCI+Lmk6ZLGR8QOSS8ANkfEAwCSvg0cC0yrtg0wH/ibdL+rSJ5tfmn6fPYPAWcDb2heaDYSc2dNcXIxq9NQg2xa9e+q2X0204C+suk+YMogyx4Cpg6xTXdE7KxYF+CDwL8BW0e15GZmLaKIg2yaXbN5omK6CxgYZlkt87uAAUkvBA6JiI9Kml9roSZOHN+Wtw/v6enOuwijzjEVRzvG1Sox7T+1m/s2PvN8+nlTuusqYzPianay2QhMLpvuATYNsmwqsAHYNcg220pNcGXrLgIOlPQj4DlAj6TlEfGJoQrVjg9E6unppq9vW97FGFWOqTjaMa5WimnhnOdVHWSzcM7zRlzGRuIaSZJqdjPadcAJAJIOBe4pNYVFxHpgT0n7SxoDLE7XH2yb60vzSfpmVkfEZyLipemggdPTeUMmGjOzoiniIJum1mwi4g5JP5H0Y5Iay6mSTgEejYhvAu8h6ewfAK5IBws8ULlNurtzgCslnQkEcHUzYzEzy1PRBtl0DQwMDL9Wm+vr29Z2H0IrVflHi2MqjnaMqx1jgoab0bpqXdd3EDAzs8w52ZiZWeacbMzMLHNONmZmljknGzMzy5yTjZmZZc7JxszMMudkY2ZmmXOyMTOzzDnZmJlZ5pxszMwsc042ZmaWOScbMzPLnJONmZllzsnGzMwy52RjZmaZc7IxM7PMOdmYmVnmnGzMzCxzTjZmZpa5sXkXwCwPves2sXrNfWx4eAfTJ41n0byZzJ01Je9imbUtJxvrOL3rNnHRtXfvnl7ft333tBOOWTbcjGYdZ/Wa+waZf39Ty2HWSZxsrONseHhH1fkbH9ne5JKYdQ43o1nHmT5pPOv7nplYpu23dw6lsU7UiX2GrtlYx1k0b+Yg8w9obkGsI5X6DNf3bad/YGB3n2Hvuk15Fy1TrtlYxymdQa5ecz8bH9nOtP32ZtG8A9r+zNJaw1B9hu38HXSysY40d9aUtv7DttbVqX2GTjZm1nHy7DPp1D5D99mYWUfJu8+kU/sMXbMxs46Sd59Jp/YZOtmYWUdphT6TTuwzdDOamXWU6ZPGV53f7n0meXPNpkN14kVlZpD0mZTfG+/p+e3dZ5I3J5sO5BtRWifr1D6TvDnZdKC8O0jN8taJfSZ5c59NB2qFDlIz6yxNr9lIOgs4BhgHLIuItWXL5gGfTJd9IyLOHmwbSZOBLwP7AuuBkyLicUlvAv5fusvvR8QHmhRaYXTqRWVmlp+m1mwkLQDmRMQRwNuAT1Ws8iXgROAVwOskHTjENucBl0XE4cB9wEmSxgGfIElMhwPzJR2ccViF06kXlZlZfprdjLYAWAUQET8HpksaDyDpBcDmiHggIvqBbwPHDrHNfODadL+rgIUR8Rjwsoj4Q0QMAJuB5zQruKKYO2sKy5bMZkbPBMbs0cWMngksWzLbbdhmlplmN6NNA35SNt0HTAHuTZf1lS17CJg+xDbdEbGzbN2pABHxKEBao3kesJZhTJw4nrFjx9QRTmvr6ekedNnio7tZfPRBTSzN6BgqpqJqx5igPeNqx5igOXE1O9k8UTHdBQwMs6yW+eX7QdILgauAkyOicvtn2LKleod5kfX0dNPXty3vYowqx1Qc7RhXO8YEjcU1kiTV7Ga0jcDksukeYNMgy6YCG4bYZlupCa5sXSTNIGlee3tE3DnaAZiZ2cg1O9lcB5wAIOlQ4J5SU1hErAf2lLS/pDHA4nT9wba5vjQfeAOwOn19KXB6RNzenJDq07tuEysu6eW0c7/Pikt62/4pfWbW2ZrajBYRd0j6iaQfA7uAUyWdAjwaEd8E3kPS2T8AXBERDwAPVG6T7u4c4EpJZwIBXC3pRcCfAx+R9JF0vU9FRGkgQUvwFfxm1mm6BgYGhl+rzfX1bWvqh7Dikt6q17nM6JnAWaceNirv0Y7ty46pONoxrnaMCRrus+mqdV3fQSAHvoLfzDqNk00OfItzM+s0TjY58BX8ZtZpfNfnHPgW5/nz83zMmsvJJie+xXl+PBrQrPncjGYdZ6jn+ZhZNlyzsY7j0YDF52bQ4nGysY7j5/kUm5tBi8nNaNZxPBqw2NwMWkyu2VjH8WjAYnMzaDE52VhH8mjA4nIzaDG5Gc3MCsXNoMXkmk2dPBomX/78O5ebQYvJyaYOHg2TL3/+5mbQ4nEzWh08GiZf/vzNisfJpg4eDZMvf/5mxeNmtDq0w2iYRvs88uwzaYfP34rNfYYj55pNHYo+GqbU57G+bzv9AwO7+zx6121qyvaNKvrnb8WW9/e/qFyzqUPRR8MM1edRSwyNbt+oon/+rcBn5vXL+/tfVE42dSryaJhG+zxaoc+kyJ9/3jyarzGt8P0vIiebDtRon4f7TPLXSM3EZ+aN8fe/PjX12Uh6vaR3SOqSdK2kzZJOzrpwlo1G+zzcZ5KvRvsMfGbeGH//61NrzWY5sBBYDPQDfwp8HfhKRuWyDDXa5+E+k3w1WjPxmXlj/P2vT63JZmdEbJN0HHB5RDwoaWeWBbNsNdrn4T6T/DRaM1k0b+Yf9dk8Pb9zzswbHSDh7//I1ZpsxklaDiwB3i/pIGB8dsUys8E0WjPp9DNzD5DIR63J5v8A7wTeFhF/SGs4H8yuWGY2mNGomXTymbkHSORjyGQj6QXpy53Av5bN+07G5TKzQcydNYXfPPgoN9+1gSef6mfPMXtw1Mum+4eyRh4gkY/hajY3AgNAV5VlA8ALqsy3DuCLAvPTu24TN96xfvf0k0/1c+Md6znoT/ap+Rh08vHzAIl8DJlsIuL5gy2TdMToF8eKwG3e+Wq0GajTj58HSCR2n3A8soPp+2V/wlFTn42k5wBvASals54FvBV4XkblshbmNu98NdoM1OnHr9MHSEA+Jxy1DhC4GvgtcCzwLeA1wBmZlMhqkmcziNu889VoM5CPX2cPkIB8TjhqvevzuIhYBtwfEe8F5gNvyqRENqy87zo7fVL1Ue9u826ORq9g9/GzPE44ak02z5Y0GdhD0uSIeBg4MLNS2ZDyflKlb9eRr7mzprBsyWxm9ExgzB5dzOiZwLIls2s+I22F49e7bhMrLunltHO/z4pLen17/ibL44Sj1ma0lcCJJLenuVvSQ8CvMiuVDSnvZhC3eeevkWagvI9fpw9QaAV5DJKoKdlExMrSa0k3ABMj4s7MSmVDaoWhm53e5l10eR6/Th+g0AryOOGodTTaWVXmvTEiPjz6RbLheOimFVneNXNLlE44enq66evblvn71dqM9lTFNnNIRqeNWJq4jgHGAcsiYm3ZsnnAJ9Nl34iIswfbJu1D+jKwL7AeOCkiHpe0lOQu1eOACyLi0nrK2crybgYxa0Qr1Myt+WptRvtY5TxJF470zSQtAOZExBGSDgY+BxxVtsqXSJLKg8AaSVcC+w+yzXnAZRFxtaTzgZMkXQOcDxwKPAn8WNLXIuIPIy1rq3MzlhWVa+adqdbRaH9E0lhgVh2bLgBWAUTEz4Hpksan+3wBsDkiHoiIfuDbJNf1DLbNfODadL+rSJ63Mwe4PSIejYgdwK3AkfXEaGbZaHQ0nRVTrX02D5DcCw2S+6Q9B7isjvebBvykbLoPmALcmy7rK1v2EDB9iG26I2Jn2bpTB9nH1OEKNXHieMaOHTOiQIqgp6c77yKMOsdUHEPFtfjobhYffVATSzM6OvFYjZZa+2z+vOz1ALAlIurpUXqiYrqLp5PYYMtqmT/cukPasqV6h2WRNavTr5kcU3G0Y1ztGBM0FtdIktRwjxh46xDL9oiIy2svFgAbgcll0z3ApkGWTQU2ALsG2WabpPFpc1lp3Wr7+N4Iy2hmZqNsuJrNa9L/JwF/CvSS9PPMI+kPuXyE73cd8M/A5yUdCtxTagqLiPWS9pS0P8kAgcXAUpLk8YxtJF0PnABcCbwBWA3cBhwiaR+SEXRzgXeNsIzWBJ18i3uzTjTcIwZOBpD0dUARsT2d7ia5q8CIRMQdkn4i6cckNZZTJZ0CPBoR3wTeQ9LZPwBcEREPAA9UbpPu7hzgSklnAgFcHRG7JK0Abgb6gbPK+nWsRfgKcrPO0zUwMGyXBpJui4jDKub9KCIOz6xkTdTXt234D6FgWrl9ecUlvVWvs5jRM4GzTj2syhaJVo6pXu0YE7RnXO0YEzTcZ1PtwZpV1TpA4G5JPwTWpNNzgd+MtGBm4CvI24GbQW2karrOJiLeDnyMpAP+d8C/kDw8zWzEfIv7Ysv7ERdWTEMmG0kvT/9/FUmH+13AncDjJBdVmo1YK9zi3uqX9yMurJiGa0Y7mSS5/GOVZQN4WLHVwfd2KzY3g1o9hhuNdmb6/4LSPEl7ABMiYmvGZbM25nu7FZdvpGn1qKnPRtIySWdI2gu4A3hQ0vuyLZqZtSI3g1o9ah2NdgrJDS3fSNKsNhe4ieTOy2bWQdwMavWoNdk8ll4weRxwVUQ8IckXS+bIQ08tT402g/r723lqTTYDkj4HvApYJulwYK/simVD8RX4VmT+/namWp9ncxLwK+C1EfEEMBPfcyw3HnpqRebvb2eq9aLOjSSPgV6YzroTWJdVoWxoHnpqRebvb2eqdTTaZ4C/IqnhAPwF8G9ZFcqG5ivwrcj8/e1MtTajHRoRfwlsA4iIfwZemlmpbEgeempF5u/v6Ohdt4kVl/Ry2rnfZ8UlvS1/u6BaBwg8lf4/ACBpDLUnKhtlHnpqRebvb+OKOMii1mRzu6RLgOnp82OWAt/Prlg2HF+Bb0Xm729jhhpk0aqfa63J5lxgAbAVmAF8mmSQgJmZNVkRB1kMd9fnoyT9juRJmP8EXJjeL206cEsTymdmZhWKOMhiuJrNvwDHRMTdkhYBX0xvxLkFGPyRimZmlplF82b+UZ/N0/NrH2Sx+y4Oj+xg+n7Z38VhuGSzKyLuBoiI1ZI+C7wvIr6ZWYnMzGxIjQ6yyGOAwXDJ5qmK6QecaMzM8tfIIIs8BhiMdPjyQCalMDOzpsljgMFwNZtXSvpt2fTkdLoLGIiI/TMrmZmZZSKPB+ANl2yU2TubmVkuRmOAwUgN91ho34bVzKzN5HEXh1ov6jQzszZSGmDQ09NNX9+2zN/P9zczM7PMuWZjloNmX1BnljcnG7MmK+Ide80a5WY0sybzY5GtEznZmDVZEe/Ya9YoJxuzJiviHXvNGuVkY9ZkfiyydSIPEDBrsnZ4LPLu0XQP72D6JI+ms+E52ZjloNkX1I0mj6azergZzcxGxKPprB6u2ZjVodFmpCJf1OnRdFYPJxuzEWq0GanozVB53J7eiq+pyUbSGOBzwMEkz8Q5KSLurVhnKbAcGAdcEBGXDradpBcBFwPjgbXA6RExIOkM4GRgDPCliLigORFaJ2j0KYd5PCVxNOVxe3orvmb32bwV6I+II4BzgI+VL5TUDZwPHAccASyXNGGI7S4GlkfEHKAHWCDpAOCd6favBM6UtG/mkVnHaLQZqejNUHNnTWHZktnM6JnAmD26mNEzgWVLZhciUVp+mt2MtgD4avr6epLaSrk5wO0R8SiApFuBI6ttJ2kv4KCI6E3nrwIWAh8EjoyIXek+dpLUfH6fSUTWcRptRmqHZqjSaDqzWjU72UwD+gAiYpekMZLGRMRTlctTDwFTq21HUpPZXLluRPQDWwEkHQs8FBEbhirUxInjGTt2TMPBtZqenu68izDqWiGmv174Ys674o4q81VT+RrdviiyjOXmO9dzzY2/5rebtrH/lG7edMwLOerlMzJ7v5J2Oj7lmhFXZslG0mnAaRWzX1pl1YGy109ULOtKl1fOB9g1yLql9z8c+BRw/HBl3bKlerNGkRXx+o3htEpML5mxD8uWzH7GRZkvmbFPTeVrdPsiyPJYVQ6wuG/jVs674g62bn0s09pWq3z/RlsjcY0kSWWWbCJiJbCyfJ6klcDk9PVewJNpTaRkY2l5airwvfL5pe1Iajr7Vqy7IV3npcClwOKI+O3oRfU0X0Hd2RptRiryRZ15K/oAi07V7AEC1wGvT18vBm6oWH4bcIikfdKBAXOBW6ptlyapuyTNS+cvBVanTWyXA2+MiHuyCKJ0ZrW+bzv9AwO7h672rtuUxduZWZmiD7DoVM3us/kWsETSWmAH8GYASR8AboqINZJWADcD/cBZEbFTUtXtgPcDl0kaC/wgIm5N+2kOAD4vqfS+yyPittEKwmdW+XPNsnO1wwCLTtQ1MDAw/Fptrq9v24g+hNPO/T79VT63MXt0cfHyBaNWrka0Y/NMKabKNvuSIg6/bcfjBM3tsynJ+vj7WFXdtqvWdX1vtDr4eST58r25Opuv8ykm366mDr6COl9uszdf51M8TjZ1aIfnkRSZ2+zNisfJpk4+s8qPa5ZmxeNkY4XjmqVZ8TjZWCG5ZmlWLB6NZmZmmXPNxqwOvqjUGtGJ3x8nG7MRKvqTNi1fnfr9cTOa2Qj5olJrRKd+f5xszEbIF5VaIzr1++NkYzZCvl2RNaJTvz9ONmYjtGjezEHm+6JSG16nfn88QMBshHxRqTWiU78/TjZmdfBFpdaITvz+uBnNzMwy52RjZmaZc7IxM7PMOdmYmVnmnGzMzCxzTjZmZpY5JxszM8uck42ZmWXOycbMzDLnZGNmZplzsjEzs8w52ZiZWeacbMzMLHNONmZmljknGzMzy5yTjZmZZc7JxszMMudkY2ZmmXOyMTOzzDnZmJlZ5pxszMwsc2Ob+WaSxgCfAw4GuoCTIuLeinWWAsuBccAFEXHpYNtJehFwMTAeWAucHhED6X66gB8CN0TER5sRn5mZVdfsms1bgf6IOAI4B/hY+UJJ3cD5wHHAEcBySROG2O5iYHlEzAF6gAVluzsN2CvDWMzMrEbNTjYLgFXp6+uB+RXL5wC3R8SjEbEDuBU4stp2kvYCDoqI3nT+KmAhgKRJwF8DF2UUh5mZjUBTm9GAaUAfQETskjRG0piIeKpyeeohYGq17UhqMpurrAtwLvBB4CW1FGrixPGMHTumvohaWE9Pd95FGHWOqTjaMa5STDffuZ5rbvw1v920jf2ndPOmY17IUS+fkXPp6teMY5VZspF0GklTVrmXVll1oOz1ExXLutLllfMBdlVbV9JRwGMR0SuppmSzZcuOWlYrlJ6ebvr6tuVdjFHlmIqjHeMqxdS7bhMXXXv37vn3bdzKeVfcwdatjzF31pQcS1ifRo7VSJJUZskmIlYCK8vnSVoJTE5f7wU8GRH9ZatsLC1PTQW+Vz6/tB1JTWffinU3ACcA8yT9iKT28yxJ/xMRXxm96MysU61ec98g8+8vZLJplmY3o10HvD79fzFwQ8Xy24BDJO0DPAXMBd4FdFduFxH9ku6SNC8i1gBLgU9ExK2lnUk6BZjpRGNmo2XDw9VbQjY+sr3JJSmWZg8Q+BYwTtJa4D2ko8okfSBNGk8AK4CbgVuAsyJi52DbAe8HPiPpDuA35YnGzCwL0yeNrzp/2n57N7kkxdI1MDAw/Fptrq9vW9t9CO3cZt5O2jEmaM+4BuuzKVm2ZHYhm9Ea7LPpqnXdZjejmZkVWimhrF5zPxsf2c60/fZm0bwDCplomsnJxsxshObOmuLkMkJONjnpXbeJ1WvuY8PDO5g+aTyL5s30l9fM2paTTQ4q23zX923fPe2EY2btyHd9zsFQ4/TNzNqRk00OPE7fzDqNk00OPE7fzDqNk00OFs2bOcj8A5pbEDOzJvEAgRx4nL6ZdRonm5x4nL6ZdRI3o5mZWeacbMzMLHNONmZmljknGzMzy5yTjZmZZc7PszEzs8y5ZmNmZplzsjEzs8w52ZiZWeacbMzMLHNONmZmljknGzMzy5yTjZmZZc53fW5BkuYD1wB3p7N+BvwT8GVgX2A9cFJEPC5pKbAcGAdcEBGXShoDfA44GOhK171X0ouAi4HxwFrg9IgYkPQu4OR0/j9ExHdGOZ6DgVXApyPiQkmTmxmLpL2By4AZwHbgxIjYPMoxXQ78GfBIusp5EbG6SDGlcZ0DLAD2BM4FbqL4x6oypkUU+FhJGg9cDkwB9gbOAn5Eix8n12xa100RMT/9927gPOCyiDgcuA84SVI3cD5wHHAEsFzSBOCtQH9EHAGcA3ws3efFwPKImAP0AAskHQgsA44GFgLnS+oarSDSL+UFwI1ls5sdy3Lgjoh4JUmCeE8GMQF8sOyYrS5STGlcRwEvi4h5wLHApyn+saoWExT7WC0B1kbE0cAb03K3/HFysimO+cC16etVJAd+DnB7RDwaETuAW4EjSc7iVqXrXg/Ml7QXcFBE9Fbs42jg+oh4MiI2ARuBF49iuR8Hjgc25BhL+T5K6452TNUUKSaA/wb+Mn39e2Av4FUU+1hVi6na715hYoqIqyLiE+nkDJKazHwullBSAAADx0lEQVRa/Di5Ga11zZJ0HdBNcubRHRE702UPAVOBaUBf2TbPmB8Ru9Jqcw+wucq6WwfZxy9GI4iI2AXsklQ+u9mxlO+7NG+0YwJ4t6T3Ab8D/rZIMZXF9Yd08jTgO8CSNjhWlTFBwY8VgKTedF/HA7e0+nFyzaY1/Ro4m6Rt+WRgJUnbakkXMAA8UbHdYPMBdtW4bml+lsrfsxmxVM7PIr6vAB9KmzZuJ2lHL2RMkk4A3knSNNIWx6oiprY4VhExF1gKXFVRppY8Tk42LSgiHoyIr0ZEf0TcS3L2NSHtGITkLGIDSZV2ctmmz5ifVo+fJDkL2XeodSvmZ2lbk2Mpnz+NDOKLiBsjYm06uRqYXUt5Wi0mSQuBFcBxEfF72uBYVcZU9GMl6RWS9geIiB+T/I5vb/Xj5GTTgiT9laSPpq8nkYw6WQmckK7yBpI/ktuAQyTtk3b8zQVuAa4DXp+uuxi4ISL6gbskzUvnL033cQOwUNKekqYDz42IX2Uc4vVNjqV8H6X3G1WSvibpkHTyz4GfFy0mSfsAnwKOj4jSSK1CH6tqMbXBsXolaYe8pCkkTe3fpsWPk/tsWtO3gRMl3UpyQnA6cCdwpaQzgQCuTttbVwA3A/3AWRGxU9K3gCWS1gI7gDen+30/cJmkscAPIuJWAEmXkjQn9DMKo5rKSfoz4JPATOBJSX8BnAT8exNjuSj97NaStC+fmEFMK4CVknYA24B3RMQTRYkpdSIwEbi6rD/qbcCXinqsBomp6MfqovS9bwGeRfL7cAfN/X0YcUx+no2ZmWXOzWhmZpY5JxszM8uck42ZmWXOycbMzDLnZGNmZplzsjFrMkkzJa0f5X2+WtIPRnOfZqPJycbMzDLnizrNWoSkdwDvAnaSXCh3akQ8KumdJDeRfCz9d2JE/F7S64F/JrlVSORUbLOauGZj1gLSe12dBRwTEUcBvwXOTBc/GzghvXHk/cBb0vkXAm+MiNcATzW5yGYj4mRj1hoOJXn2yNZ0+kbgFenrx4BrJN1E8tyQSZL2A54dEb8sW9+sZTnZmLWGyqejdgH9kp5P8jTFv0xrNt8cZH3/LVtL8xfUrDXcDrxC0nPS6eNInis/EeiLiI3pHcBfQ3LzxUdIbgJaurvkaDyp0ywzvhGnWZNJmknSob+mbPZtJE9HfRfJY6cfIHnY107gGyTPDPkf4BqSvpoTSR4JvAK4B/glcFja32PWcpxszMwsc25GMzOzzDnZmJlZ5pxszMwsc042ZmaWOScbMzPLnJONmZllzsnGzMwy97+deWV4JH27MAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x6acf0410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import leastsq\n",
    "\n",
    "def func(params, L):\n",
    "    a,b,c=params\n",
    "    return a*L*L+b*L+c\n",
    "\n",
    "def error(params,L,S):\n",
    "    return func(params,L)-S\n",
    "\n",
    "def slovePara():\n",
    "    p0=[1,1,1]\n",
    "    Para=leastsq(error,p0,args=(x,y))\n",
    "    return Para\n",
    "\n",
    "def solution():\n",
    "    Para=slovePara()\n",
    "    a,b,c=Para[0]\n",
    "    print('Least squares estimate values:',a,b,c)\n",
    "    plt.plot(x, y, \"k.\",label='Original Data',markersize=12)\n",
    "    plt.plot(x,func(Para[0],x),'r',label='Fitted Line')\n",
    "    plt.xlabel(\"Load\")\n",
    "    plt.ylabel(\"Diflection\");\n",
    "    plt.legend()\n",
    "    plt.title('Least Squares Estimate for Linear Regression')\n",
    "    plt.show()\n",
    "    \n",
    "    residuals=y-func(Para[0],x)\n",
    "    plt.plot(x, residuals, \"o\")\n",
    "    plt.xlabel(\"Load\")\n",
    "    plt.ylabel(\"Residuals\");\n",
    "    plt.title('Residuals versus Load')\n",
    "\n",
    "solution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The residuals randomly scattered around zero, indicate that the quadratic is a good function to describe these data. There is also no indication of non-constant variability over the range of loads."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "[1] https://en.wikipedia.org/wiki/Probability_density_function<br>\n",
    "[2] https://ned.ipac.caltech.edu/level5/March02/Silverman/Silver1.html<br>\n",
    "[3] https://github.com/fonnesbeck/scipy2015_tutorial/blob/master/notebooks/2.%20Density%20Estimation.ipynb<br>\n",
    "[4] https://en.wikipedia.org/wiki/Bayes_estimator<br>\n",
    "[5] http://gradquant.ucr.edu/wp-content/uploads/2013/11/Maximum-Likelihood-and-Method-of-Moments-Estimation.pdf<br>\n",
    "[6] http://research.cs.tamu.edu/prism/lectures/pr/pr_l7.pdf<br>\n",
    "[7] http://www.cmrr.umn.edu/~kendrick/statsmatlab/StatsLecture3Slides_ModelSpecification.pdf<br>\n",
    "[8] http://www.itl.nist.gov/div898/handbook/pmd/section6/pmd61.htm <br>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}