Lecture 03_Maximum Likelihood Estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center><font color=navy>Big Data Economics</font></center>\n",
    "### <center>Maximum Likelihood Estimation (MLE)</center>\n",
    "#### <center>Ali Habibnia</center>\n",
    "\n",
    "\n",
    "<center> Assistant Professor, Department of Economics, </center>\n",
    "<center> and Division of Computational Modeling & Data Analytics at Virginia Tech</center>\n",
    "<center> habibnia@vt.edu </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "### <center>Linear Regression : The Probabilistic Perspective</center>\n",
    "\n",
    "Probability is a great way to express the relationship between the data and the model, concisely. We are familiar with how Linear Regression works. We know the loss function just tries to minimize the quadratic distance between data points and the model (line). Today, we will revisit Linear Regression from a probabilistic perspective, using a method known as the Maximum Likelihood estimation.\n",
    "\n",
    "<img src=\"images/lr2.png\"  width=\"600\">\n",
    "</div>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Readings:\n",
    "\n",
    "1. ***Chapter 8.2*** [The Elements of Statistical Learning: Data Mining, Inference, and Prediction](https://web.stanford.edu/~hastie/ElemStatLearn/printings/ESLII_print12.pdf). \n",
    "2. For a quick review see: ***Chapter 24.1***, Understanding Machine Learning From Theory to Algorithms.\n",
    "3. Nando de Fretais’s [lecture](https://www.youtube.com/watch?v=voN8omBe2r4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In simple terms, MLE lets us choose a model (parameters) that explains the data (sample set) better than all other models. For any given regression model, the objective function can be derived based on the principle of Maximum Likelihood.\n",
    "\n",
    "- Maximum Likelihood Estimation or MLE is a tool based on probability. \n",
    "- Probability is a framework for meauring and managing uncertainty. \n",
    "- Every inference we make, has some degree of uncertainty associated with it. \n",
    "- It is essential for us to quantify this uncertainty."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normal (Gaussian) Distribution\n",
    "\n",
    "A probability distribution is a function that provides us the probabilities of all possible outcomes of a stochastic process. The most commonly occurring distribution is the Gaussian Distribution or the Normal Distribution.\n",
    "\n",
    "The gaussian distribution is a means to measure the uncertainty of a variable that is continuous between −∞\n",
    " and +∞. The distribution is centered at mean, μ. The width depends on the parameter σ, the standard deviation (variance, $σ^2$). Naturally, area under the curve equals 1.\n",
    "\n",
    "<img src=\"images/normal1.png\"  width=\"600\">\n",
    "\n",
    "<img src=\"images/normalpdfcdf.gif\"  width=\"550\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ##### Do you remember what is the \"Standard Normal Distribution\"?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How is Gaussian Distribution relevant to linear regression? \n",
    "\n",
    "The data points in the sample set, do not accurately represent the original data generating distribution or process. Hence we consider the process stochastic and build our model to accomodate a certain level of uncertainty. Every data point can be considered a random variable sampled from the data generating distribution which we assume to be gaussian. By that logic, estimation is basically recreating the original distribution, that generated the sample data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Sampling\n",
    "The process of sampling from a normal distribution is expressed as, $x∼N(μ,σ^2)$. ‘$x$’ is a random variable sampled or generated or simulated from the guassian distribution. As we sample from this distribution, most samples will fall around the center, near the mean, because of higher probability density in the center.\n",
    "\n",
    "<img src=\"images/sampling1.png\"  width=\"600\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Gaussian\n",
    "\n",
    "The 2D gaussian distribution or bivariate distribution, consists of 2 random variables $x_1$ and $x_2$. It can have many different shapes. The shape depends on the covariance matrix, $Σ$. The Multivariate Gaussian Distribution is a generalization of bivariate distribution, for ‘$n$’ dimensions.\n",
    "\n",
    "<img src=\"images/bivariate3d.png\"  width=\"600\">\n",
    "\n",
    "<img src=\"images/bivariatepdf.png\"  width=\"350\">\n",
    "\n",
    "In the above expression, $y$ and $μ$ are vectors of ‘$n$’ dimensions, and $Σ$ is a matrix of shape ‘$n$x$n$’.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum Likelihood\n",
    "\n",
    "Consider 3 data points, $y_1=1$,$y_2=0.5$,$y_3=1.5$, which are independent and drawn from a guassian with unknown mean θ and variance 1. \n",
    "\n",
    "Let’s say we have two choices for θ: {1, 2.5}. Which would you choose? Which model (θ) would explain the data better?\n",
    "\n",
    "In general, any data point drawn from a gaussian with mean θ and and variance 1, can be written as,\n",
    "\n",
    "### $y_i∼N(θ,1)=θ+N(0,1)$\n",
    "\n",
    "θ, the mean, shifts the center of the standard normal distribution (μ=0 and σ2=1).\n",
    "\n",
    "The likelihood of data (y1,y2,y3) having been drawn from N(θ,1), can be defined as,\n",
    "\n",
    "### $P(y_1,y_2,y_3|θ) = P(y_1|θ)P(y_2|θ)P(y_3|θ)$\n",
    "\n",
    "as $y_1,y_2,y_3$ are independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have two normal distributions defined by θ = 1 and θ = 2.5. Let us draw both and plot the data points. In the figure below, notice the dotted lines that connect the bell curve to the data points. Consider the point $y_2=0.5$ in the first distribution( $N(μ=1,σ^2=1)$ ). The length of the dotted line gives the probability of the $y_2=0.5$ being drawn from $N(μ=1,σ^2=1)$.\n",
    "\n",
    "<img src=\"images/toy-eg.png\"  width=\"400\">\n",
    "\n",
    "The likelihood of data $(y_1,y_2,y_3)$ having been drawn from $N(μ=1,σ62=1)$, is given by,\n",
    "\n",
    "### $P(y_1,y_2,y_3|θ=1) = P(y_1|θ=1)P(y_2|θ=1)P(y_3|θ=1)$\n",
    "\n",
    "The individual probabilities in the equation above, are equal to the heights of corresponding dotted lines in the figure. \n",
    "\n",
    " - We see that the likelihood, given by the product of individual probabilities of data points given model, is basically the product of lengths of dotted lines. \n",
    " \n",
    "It is obvious that the likelihood of model θ=1 is higher. We choose the model (θ=1), that maximizes the likelihood.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression\n",
    "\n",
    "We are ready to apply MLE to linear regression, where the objective is to find the best line that fits the data points. But first, let us make some assumptions. We assume each label, $y_i$, is gaussian distributed with mean, $x_i^Tθ$ and variance, $σ^2$, given by\n",
    "\n",
    "<img src=\"images/reg.png\"  width=\"300\">\n",
    "\n",
    "The mean, $x_i^Tθ$ represents the best fit line. The data points will vary about the line, and the second term, captures this variance, $N(0,σ^2)$\n",
    "\n",
    "<img src=\"images/lr1.png\"  width=\"600\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimation\n",
    "\n",
    "If we assume that each point $y_i$ is gaussian distributed, the process of estimation becomes the process of maximizing the product of the individual probabilities, which is equivalent to maximizing the log likelihood. We switch to log space, as it is more convenient and it removes the exponential in the gaussian distribution.\n",
    "\n",
    "\n",
    "<img src=\"images/lr2.png\"  width=\"600\">\n",
    "\n",
    "As the data points are independent, we can write the joint probability distribution of $y,θ,σ$ as,\n",
    "\n",
    "<img src=\"images/ml1.png\"  width=\"500\">\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "rewriting in matrix form,\n",
    "\n",
    "<img src=\"images/ml2.png\"  width=\"400\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Log likelihood,\n",
    "\n",
    "<img src=\"images/ml3.png\"  width=\"400\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first term is a constant and the second term is a parabola, the peak (maxima) of which can be found by equating the derivative of $l(θ)$ to zero. Equating first derivative to zero, we get,\n",
    "\n",
    "<img src=\"images/ml4.png\"  width=\"350\">\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we get,\n",
    "\n",
    "<img src=\"images/ml5.png\"  width=\"200\">\n",
    "\n",
    "Finally, we reach our goal of finding the best model for linear regression. This equation is commonly known as the normal equation. The same equation can be derived using the least squares method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Showcase:  Maximum likelihood estimation to fit a linear (polynomial) model to some data points\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "\n",
    "x = np.linspace(0.0,10.0, num=100)\n",
    "a = 4.0\n",
    "b = -3.5\n",
    "c = 0.0\n",
    "y = (a*(x**2)) + (b*x) + c\n",
    "\n",
    "#   let's add noise to the data\n",
    "#   np.random.normal(mean, standardDeviation, num)\n",
    "noise = np.random.normal(0, 10., 100)\n",
    "y = y+noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X9w3Hed3/HnO8IRUXQ0eIW3xolW3BGuF26GH8rkQmk7sSNKSJkzNz1aZySR0njUJoHmDtxrcvoDrq2mtBXQ3BxJz+cEjKWiywAdMpAjFyVyGRgItSEEB18mJmgdg8+OZOAinDrBfvcPfVd8tfp+d7/7e/e7r8fMjna/+/3u9/ONnPd+9Pm+P++PuTsiIpJeF7W6ASIi0lgK9CIiKadALyKScgr0IiIpp0AvIpJyCvQiIimnQC8iknIK9CIiKadALyKScq9odQMABgYGfGhoqKpjf/GLX3DppZfWt0FtTtfcHXTN3aGWaz58+PCSu7+m3H5tEeiHhoY4dOhQVccePHiQ6667rr4NanO65u6ga+4OtVyzmeWT7KehGxGRlFOgFxFJOQV6EZGUU6AXEUk5BXoRkZRToBcRaYHZ2VmGhobYsWMHQ0NDzM7ONuxcbZFeKSLSTWZnZ5mYmODs2bMA5PN5JiYmABgdHa37+dSjFxFpssnJybUgX3D27FkmJycbcj4FehGRJjt+/HhF22ulQC8i0mSDg4MVba+VAr2ISJNNTU3R19e3bltfXx9TU1MNOZ8CvYhIk42OjrJ3715yuRxmRi6XY+/evQ25EQvKuhERaYnR0VFGR0ebUshNPXoRkSYp5M5fdNFFDc+dD1OgFxFpgkLufD6fx93Xcufn5+cbfm4FehGROovqucflzu/bt6/h7dEYvYhIHcXNei0O8gWnT59ueJvUoxcRqaO4nntPT0/k/lu2bGl4mxToRUTqKG526/nz5yNz53fv3t3wNpUN9Gb2SjP7tpl9z8yeMrM/CbZ/xsx+ZGZPBI83B9vNzP7UzI6Z2ZNm9tZGX4SISLuIm91ayJUvzp0fGRlpeJuS9OjPATvc/U3Am4EbzOza4L1/7+5vDh5PBNveBVwZPCaAe+vdaBGRdlVq1uvo6CiLi4tcuHCBxcXFhk2QKlY20PuqleDlpuDhJQ7ZCXw2OO5bwGVmtrX2poqItL9mz3pNItEYvZn1mNkTwGngEXd/PHhrKhie+aSZ9QbbtgHPhQ4/EWwTEekKreq5xzH3Up3zop3NLgP+N/BBYBn4W+BiYC/wQ3f/j2b2FeC/uPvXg2MeBf7I3Q8XfdYEq0M7ZLPZ4bm5uaouYGVlhf7+/qqO7VS65u6ga+4OtVzz9u3bD7v71WV3dPeKHsBHgD1F264Dvhw8/3PgptB7TwNbS33m8PCwV2thYaHqYzuVrrk76Jq7Qy3XDBzyBHE7SdbNa4KePGZ2CTAC/E1h3N3MDHgPcCQ45EHgfUH2zbXAz939ZNJvKBERqa8kM2O3AvvNrIfVMf0H3P3LZvaYmb0GMOAJ4N8G+z8E3AgcA84C769/s0VEJKmygd7dnwTeErF9R8z+Dtxee9NERKQeNDNWRCTlFOhFRFJOgV5EJOUU6EVEqhCuOT8wMMDAwEDkylGtWlUqTPXoRUQqVFxzfnl5ee29Qv35gqja9EBTZ8sq0IuIVCiq5nzY2bNnmZycXHse9Z4CvYhIG4urOZ90nyTH15PG6EVEKhRXc754n7j9khxfTwr0IiIlRN1Mjao5H1aoP1+qNn0zKdCLiMQo3HTN5/O4+7qbqeGa85lMhkwms6H+fLvUptcYvYhIkdnZWSYnJ8nn8xveO3v2LGNjY+RyubVVo0opBPxWUqAXEQkpTp2M06pUyWpo6EZEJKRc6mRYOI2ynSnQi4jwq5uuUcM1pTQ7VbIaGroRka6XdLgmSrNTJauhHr2IdL1ywzV9fX3ceuutbZEqWY0kSwm+0sy+bWbfM7OnzOxPgu2vM7PHzewZM/tLM7s42N4bvD4WvD/U2EsQEalNqeGXQkrkPffc0xapktVIMnRzDtjh7itmtgn4upn9FfAh4JPuPmdm/xO4Bbg3+PlTd3+9me0C/ivwLxvUfhGRmg0ODkaOzedyORYXF9det0OqZDXK9uiDxcZXgpebgocDO4DPB9v3s7pAOMDO4DXB+9cHC4iLiLSldpnB2ii2usRrmZ1WFwY/DLwe+BTw34Fvufvrg/evAP7K3X/bzI4AN7j7ieC9HwK/4+5LRZ85AUwAZLPZ4bm5uaouYGVlhf7+/qqO7VS65u6ga26u+fl59u3bx+nTp9myZQu7d+9mZGSk4eet5Zq3b99+2N2vLrujuyd+AJcBC8A/Bo6Ftl8BfD94/hRweei9HwKZUp87PDzs1VpYWKj62E6la+4OuubuUMs1A4c8QeyuKOvG3X8GHASuBS4zs8IY/+XAT4LnJ4LAT/D+3wPOVHIeERGpnyRZN68xs8uC55cAI8BRVnv2vx/sdjPwpeD5g8FrgvcfC755RESkBZJk3WwF9gfj9BcBD7j7l83sB8Ccmf1n4LvAfcH+9wEHzOwYqz35XQ1ot4iIJFQ20Lv7k8BbIrY/C1wTsf3/Ae+tS+tERKRmmhkrIl0laiGRtFOgF5GuEbWQyPj4OGaW6qCvQC8iXSOqpk0hV6RQXz6NwV6BXkS6RrmSwp1SX75SCvQi0jWSlBTuhPrylVKgF5GOVsnN1aiaNsU6ob58pRToRaRjRd1cLTXOPjo6ulZqGKC43mKaCpmFKdCLSMeKurlabpx9dHSUxcVF3J0DBw50ZH35SinQi0jHihtPP378+LohnYGBAQYGBjYM7xSC/oULF1hcXExlkAetGSsiHSxuwZDNmzevWwN2eXl57b3C8A6Q2sBeTD16Eek4hd56Pp+PHGcHSq4Bm9Y0yjjq0YtIRyncgC0EcnfHzHB3MpkMsL4HHyeNaZRxFOhFpKPEzW7NZDK8+OKLJXvyYWlMo4yjQC8iHSWuJ56kF1+Q1jTKOBqjF5GOUk1PPJPJkMlkUp9GGSfJClNXmNmCmR01s6fM7I5g+0fN7Mdm9kTwuDF0zF1mdszMnjazdzbyAkSku0TNbu3r61sbny+Wy+VYWlpiaWkp9WmUcZIM3fwS+LC7f8fMfg04bGaPBO990t2nwzub2VWsrir1RuC1wLyZvcHdz9ez4SLSnQpBenJykuPHjzM4OLg2DBO+SQvdN0QTJ8kKUyeBk8HzF8zsKLCtxCE7gTl3Pwf8KFhS8Brgm3Vor4gIo6Ojsb3y4i+Abuu9R6lojN7MhlhdVvDxYNMHzOxJM7vfzF4dbNsGPBc67ASlvxhERACYn5+vafWnbpnpWikrFN0vu6NZP/B/gCl3/6KZZYElwIH/BGx1939tZp8CvunuM8Fx9wEPufsXij5vApgAyGazw3Nzc1VdwMrKCv39/VUd26l0zd2h2655fn6e6elpzp07t7att7eXPXv2MDIy0sKWNVYtv+ft27cfdvery+7o7mUfwCbgYeBDMe8PAUeC53cBd4Xeexh4W6nPHx4e9motLCxUfWyn0jV3hzRf88zMjOdyOTczz+Vya69Z7Tiue+RyuVY3t6Fq+T0DhzxBDE+SdWPAfcBRd/9EaPvW0G6/BxwJnj8I7DKzXjN7HXAl8O2y3zgi0hXiSgtH1ayBX+XNd+Oi3vWSJOvm7cA48H0zeyLY9sfATWb2Zla/dReBfwPg7k+Z2QPAD1jN2LndlXEjIoG40sI9PT2cP78xVAwODm4oe9CNhclqkSTr5uuARbz1UIljpgDlNInIBnEzW8+fP09vb++6MfpCemSpuvMK9OVpZqyINFXczNZcLseePXsiFwIpVXdeylOgF5GmipvZOjU1xcjISGR6ZNyXQzcVJquFAr2INFV43dZStWfCN19XVla4+OKL172vWa/JqXqliDRdqZmtsLHm/PLyMps2bSKTyXDmzBnNeq2QevQi0hSVpEdG3Xx9+eWX6e/v16zXKqhHLyINV2l6pG6+1pd69CLScKXSI6Po5mt9KdCLSMOV6qGHh3R27drF7OxsycwcqZwCvYjUVdRYfFxPfPPmzevKIZw6dWptSCdJZo4kozF6EambqLH48fFx3B0zKxQ6BFjrsccN6eiGa/2oRy8idRM1Fl8I7oVgD6z10M+cORP5ObrpWl8K9CJSN+UCtLuTy+XWeuu66docCvQiUjdJAnT4y0A3XZtDgV5E6iYqcBcLfxkUl0PIZrO66doACvQiUjfhwA2sjckXRPXWw+u8zs3NKcg3gAK9iNQsnFI5OTnJ1NQU7s6BAweUItkGyqZXmtkVwGeBvw9cAPa6+91mthn4S1bXi10E/oW7/zRYevBu4EbgLPCv3P07jWm+iLRaufIGCuytl6RH/0vgw+7+W8C1wO1mdhVwJ/Cou18JPBq8BngXq+vEXglMAPfWvdUi0jYqLW8gzVc20Lv7yUKP3N1fAI4C24CdwP5gt/3Ae4LnO4HPBouUfwu4rGghcRFJERUga38VjdGb2RDwFuBxIOvuJ2H1ywDYEuy2DXgudNiJYJuIpEhhXD482zVMufDtI3EJBDPrB74A/IG7/13x3fTwrhHbNvxLMLMJVod2yGazHDx4MGlT1llZWan62E6la+4O7XzN8/PzTE9Pr1vIO6y3t5exsbGK29/O19woTblmdy/7ADYBDwMfCm17GtgaPN8KPB08/3Pgpqj94h7Dw8NerYWFhaqP7VS65u7Qztecy+Wc1Q7chkcul/OZmZmqPredr7lRarlm4JAniOFJsm4MuA846u6fCL31IHAz8LHg55dC2z9gZnPA7wA/92CIR0TSIW783cxYXFxsbmOkrCRj9G8HxoEdZvZE8LiR1QD/DjN7BnhH8BrgIeBZ4BjwF8Bt9W+2iNRbJUv9qUZNZynbo3f3rxM97g5wfcT+DtxeY7tEpIkqXepvampq3f6gGjXtTDNjRaRkLnxUT7+4Ro1mvbY3LTwiIrFj7oWeffFCImNjY+RyOaamphTcO4B69CISO7be09MTu5BI4Uug1Fi+tAcFehGJrQt//vz5ksep1EFnUKAXkdgx90K54VJU6qD9aYxeRID4SpPF2TXFlFLZ/tSjF5FY1SwkIu1HgV6kSyWdIFVYAcq1kEjHUqAX6UKFCVL5fB5335BBE/clEF72b3FxUUG+Q2iMXqQLlVsspJJZstL+1KMX6SKFnno+n498//jx41oxKoXUoxfpEsX1bKIMDg5qxagUUo9epEtE9dTDChk0qkyZPgr0Il2iVI88nEETN0tWaZSdS4FeJIXCWTMDAwMMDAzEru2ay+XWZdCoMmX6aIxeJGWKx+KXl5dj943rqcfNkpXOVLZHb2b3m9lpMzsS2vZRM/tx0YpThffuMrNjZva0mb2zUQ0XkWjlxuIL1FPvHkl69J8B/gz4bNH2T7r7dHiDmV0F7ALeCLwWmDezN7h76RJ4IlI3SbJjtLZrdynbo3f3rwFnEn7eTmDO3c+5+49YXTf2mhraJyIVSpIdowya7lLLGP0HzOx9wCHgw+7+U2Ab8K3QPieCbRuY2QQwAZDNZjl48GBVjVhZWan62E6la+4O1V7z2NgY09PTnDt3LvL93t5exsbG2vK/p37PDeLuZR/AEHAk9DoL9LD6F8EUcH+w/VPAWGi/+4B/Xu7zh4eHvVoLCwtVH9updM3doZZrnpmZ8Vwu52bmmUzGM5mMm5nncjmfmZmpXyPrTL/nygCHPEEMryq90t1Puft5d78A/AW/Gp45AVwR2vVy4CfVnEOkWyStIlnJ50xOTjI1NcWFCxdYWlpiaWlJhci6WFVDN2a21d1PBi9/Dyhk5DwI/C8z+wSrN2OvBL5dcytFUqo4FTJcQGzbtshRz4o/R4FdkqRXfg74JvCbZnbCzG4B/puZfd/MngS2A38I4O5PAQ8APwC+CtzuyrgRiZW0gFjUBKjwXwAqRCallO3Ru/tNEZvvK7H/FKvj9iJSRpICYqUmQOXzecbHx2NnvaoQmYBKIIi0VJICYuUmQMUF+VKfL91FgV6khZIUEKu2V65CZFKgQC/SQkkKiFXTK1d5AwlTUTORFitXQGxqaqrsgiFhhWqUIgXq0Yu0QCW588W9/kwmQyaTAVZr1oRpuEaiKNCLNFkhiyafz+Puaznvs7Oz674Adu3atfYFMDo6yuLi4roJUO7OgQMHVDdeytLQjUiTxeW833HHHbz44otr7506darspCfVjZck1KMXabK4LJrl5WVNepKGUKAXabJKs2g06UlqpUAv0mRxufOFG6zFNOlJaqVAL9Jkcbnzd999d9nJUyLV0M1YkRYodRN1cnKS48ePs2XLFj7+8Y/rZqvUTD16kTYSTqOcm5tTkJe6UKAXEUk5BXqRNjU/P1+XladEFOhFmqSSsgezs7NMT09Hzp4VqVSSFabuN7PTZnYktG2zmT1iZs8EP18dbDcz+1MzO2ZmT5rZWxvZeJF2UxzMb7vtNoaGhjAzxsfHEwfuyclJzp07t26bJk9JtZL06D8D3FC07U7gUXe/Eng0eA3wLlbXib0SmADurU8zRdpfVA2be++9l3w+D2xcIKRU4E6y8pRIUmUDvbt/DThTtHknsD94vh94T2j7Z33Vt4DLzGxrvRor0gyVDLGElVsJKkpc4E6y8pRIUtXm0Wfd/SSAu580sy3B9m3Ac6H9TgTbThZ/gJlNsNrrJ5vNcvDgwaoasrKyUvWxnUrX3Djz8/NMT0+vDZvk83luueUWjh49ysjISMljq+ltb9myJfK6xsbG1rUDoLe3l7GxsVT/7vVvu0HcvewDGAKOhF7/rOj9nwY/vwL8o9D2R4Hhcp8/PDzs1VpYWKj62E6la26cXC7nwIZHLpdbt9/MzIzncjk3M8/lcmuvo46Ne/T19fnMzExsWyYnJzecI+30b7sywCFPEMOrzbo5VRiSCX6eDrafAK4I7Xc58JMqzyHSdEnGxuPqyd94440bShgUKywUkqR2/MjIyNrkqcXFRU2ekqpVG+gfBG4Ont8MfCm0/X1B9s21wM89GOIR6QRJxsbj6sk/9NBDG2rY3HrrreteHzhwAHdX4JamKjtGb2afA64DBszsBPAR4GPAA2Z2C3AceG+w+0PAjcAx4Czw/ga0WaRhotZnLS4sVqrXr4VApB2VDfTuflPMW9dH7OvA7bU2SqRVCkG6UFhscHCQqampdcF7cHBwLWUyTBkx0q40M1akSLiw2NTUFJOTk+tSLePqyaucsLQrBXqRGHE3XYHIevIaspF2pUAvEiPupuvY2BiTk5NMTU3F9vpF2okCvQjRs2FLTYAq9O5vu+22yF6/gr20EwV66VqF4B5XcGzz5s0ljz979ix79+6N7PWr+Ji0Ey0lKF2pMP5eCNIeUXDskksuoa+vr2T9mvPnz0duV/ExaSfq0UtLVVtArFZJCpCdOXNm7aZrnJ6ensjtSrWUdqJALy0Tl9VSbbCv5EsjSY97cHBwLdVyZmYmMqVyYmJCqZbS9hTopWXislqqGd+u9EujXI+7OFiPjo5GplTec889SrWUtqdALy1Tz8U1Kv3SiJr0VK7gWHgiVbhWTdx2kXahQC8tU8/FNSr90ojqoavgmKSVAr20TD1LCVTzpaGeuHQLBXppmbhx72oCrurPiMRToJeWStKrnp2dZdeuXSWzaer5pSGSNgr00hTV5ssXsmlOnTpVNptGQzEi0WoK9Ga2aGbfN7MnzOxQsG2zmT1iZs8EP19dn6ZKp6olX75cNk3cF0irJmKJtKN6lEDY7u5Lodd3Ao+6+8fM7M7g9X+ow3mkQ5UK1uV63XFZM/l8noGBAV544QVeeumltW0TExN84xvfYP/+/WvnDJcXVi9fulEjhm52AvuD5/uB9zTgHNJGyvWea8mXL5U1s7y8vBbkC1RoTGSjWgO9A39tZofNbCLYli0sCB783FLjOaSNJRmWqSVfPiqbppy4QmOFvwIGBgY0pCNdxYqr9lV0sNlr3f0nZrYFeAT4IPCgu18W2uen7r5hnD74YpgAyGazw3Nzc1W1YWVlhf7+/qqO7VTtdM27du3i1KlTG7Zns1kKv9P5+Xmmp6c5d+5c5H67d+9mZGQk9hzz8/Ps3buX559/PlGbLrroIi5cuJBo397eXvbs2VPy/K3STr/nZtE1V2b79u2H3f3qsju6e10ewEeBPcDTwNZg21bg6XLHDg8Pe7UWFhaqPrZTtdM1m5mz+pfduoeZrdtvZmbGc7nc2nvhffv6+nxmZqbkeRYWFtaOL/Xo6+vzW2+91fv6+sruW3jkcrkG/heqXjv9nptF11wZ4JAniM9VD92Y2aVm9muF58A/BY4ADwI3B7vdDHyp2nNI+0s6LFNIfczlcpG135OMn0cN42zatIlMJhNbaCwJ1Y6XtKtljD4LfN3Mvgd8G/iKu38V+BjwDjN7BnhH8FpSqtIZqbXcmI2aFPXpT3+apaWl2EJjSYK9asdL2lUd6N39WXd/U/B4o7tPBduX3f16d78y+Hmmfs2VdlHItBkfH+eSSy7Z0KuOS2OstZBZpZOiyt3MVZkE6QaaGSsVK860WV5e5sUXX+TAgQNlg2+pvwAaMcmp+K+ATCaT+EtJJC0U6LtIvQJpLQuGxNWkAeq62lTxOQt/BSwtLUUO9YikmRYH7xLFi2HXMlu01gVDRkdHN5xzaGio6tmzIlKaevQpUa63XkkvvNxnJR1nr8carvl8nqGhIebn52OPFZEykuRgNvqhPPrKFF/zzMzMhrzx4tz0SvLdy31WvfYJK5cj39vbWzbXPm30b7s7NCOPvuVB3hXoK1Z8zXFBMjwRKMk+5fYrTHoyM89kMp7JZNzM1t5zXz8xKsn5CqK+GJIem1b6t90d2nrClLSPJGPmSfPdSw2hjI+Pl8y0CWfjVNrW8E3aSo8VkdIU6FMgyZh50hWYSuWze5kZrVH3AZK2tdDGUpOcNLFJpDoK9CmQtLceN9kofNN0ZWWFiy++OPG5w73scj3upJOToq6nt7dXE5tEqqRAnwK1rJcaNfnJ3clkMonOHe5ll+pxV9KmqOvZs2eP0ixFqpTKQF/LxKBmLEHXqBmg5XrrUUvtjY2NbRhuefnll+nv7y9bJ6a4hx73l8XMzEzFk5OKr6cdywiLdIwkd2wb/ahn1k2laX31OjapWs8xMzPj2Wx2Q7ZLJedKUsbXzCKPL6Rpxp07nJlTrn2VUDZGd9A1V4ZuTa9MmkYYpZZj3ZMFuaTniPqsSgNv3Ll6enoS12hvVOCulAJAd9A1V6ZrA33SiUFRajk2aU89yTniPiuTyZQMzkknSZV71PuvmHpQAOgOuubKJA30qRujr2Z6fmEd0dX/bsk/MyxpiYEk7Yv7rOXl5ZJtKD5f3Ll6enpiP0MVHUXSJ3WBPkkZXDPbMPknLogmTQlMWugrSSpkLRODjh8/vnad+XweM9twromJibrdNBWR9tewQG9mN5jZ02Z2zMzubNR5iiUpgwsbJ/9EqaR3W8mSeuVSIWuZGLR58+YN11kI9lFL7akuu0gXSDK+U+kD6AF+CPw6cDHwPeCquP2bUesmycLS4UeScfmwUmP0ld7QTFL3pdDG4vPFjeN3ep0Yjd12B11zZWjxGP01wDFfXW7wJWAO2Nmgc5UUHsaoRKW96nouqJGk7ksul+PAgQMbznfmTPTKjaoTI9LFknwbVPoAfh/YF3o9DvxZ3P6N6tEn7RkXP6KyTqpNM6wkZTNpSmWpkr21poi2K/X0uoOuuTIk7NGbJxirrpSZvRd4p7vvDl6PA9e4+wdD+0wAEwDZbHZ4bm6uonPMz8+zb98+Tp8+zZYtW9i9e/eG2ZO7du3i1KlTZT/rVa96FQAvvPBC5GfNz88zPT3NuXPnNhybzWYjz12wY8eOyPsBZsZjjz1W8hy9vb3s2bMHYN21jo2N8e53vzvyfKU+p5Nnl66srNDf39/qZjSVrrk71HLN27dvP+zuV5fdMcm3QaUP4G3Aw6HXdwF3xe1faY++1px1YiYYRZ0nydh+NQtqVFIHvli5HkC7THKqJ/X0uoOuuTK0csIUq2vRPgu8jl/djH1j3P6VBvpyC1sUAlstwxiVDvtUsqBG+Iuh3JdJ1E1h/c/QHXTN3aFjb8a6+y+BDwAPA0eBB9z9qXp9fqkbi+GbnUnL90ZJUls9SZtKpVMmWahDNdhFpFYNy6N394fc/Q3u/hvuXtdC4uWCX2GGaC3leyvNUnH32EqUcZUly32ZJP1SEhEppSNnxkb11IsVAnVckC2n1JdJ8WzTgiSpk1FtjKJJTCJSLx0Z6JPkmdc65FGqtnohfz1KVH2bStuYy+VUikBE6qYjAz38qqc+OTlZ9Th8uc+PG/YpnDuuZx+uN1NqcZFa7iGIiCTVsYG+YGRkpGF1W8oN+8T1yMP1ZrzEbNha7iGIiCTV8YEeqh+Hr1VcjxxIVLIYWtd2EekeqQj0rRLXI1e9GRFpJwr0NYrqkSctWSwi0gwK9A2gm6wi0k4U6BtAN1lFpJ28otUNSKtCGqaISKupRy8iknIK9CIiKadALyKScgr0IiIpp0AvIpJyDVkztuJGmD0PxK++UdoAsFTH5nQCXXN30DV3h1quOefurym3U1sE+lqY2SFPsjhuiuiau4OuuTs045o1dCMiknIK9CIiKZeGQL+31Q1oAV1zd9A1d4eGX3PHj9GLiEhpaejRi4hICR0d6M3sBjN72syOmdmdrW5Po5nZFWa2YGZHzewpM7uj1W1qBjPrMbPvmtmXW92WZjGzy8zs82b2N8Hv+22tblMjmdkfBv+mj5jZ58zsla1uUyOY2f1mdtrMjoS2bTazR8zsmeDnq+t93o4N9GbWA3wKeBdwFXCTmV3V2lY13C+BD7v7bwHXArd3wTUD3AEcbXUjmuxu4Kvu/g+AN5Hi6zezbcC/A652998GeoBdrW1Vw3wGuKFo253Ao+5+JfBo8LquOjbQA9cAx9z9WXd/CZgDdra4TQ3l7ifd/TvB8xdY/Z9/W2tb1Vhmdjnwz4B9rW5Ls5jZq4B/AtwH4O4vufvPWtuqhnsFcImZvQLoA37S4vY0hLt/DShea3QnsD94vh94T73P28liAktfAAAB40lEQVSBfhvwXOj1CVIe9MLMbAh4C/B4a1vScP8D+CPgQqsb0kS/DjwPfDoYstpnZpe2ulGN4u4/BqaB48BJ4Ofu/tetbVVTZd39JKx25oAt9T5BJwd6i9jWFSlEZtYPfAH4A3f/u1a3p1HM7N3AaXc/3Oq2NNkrgLcC97r7W4Bf0IA/59tFMCa9E3gd8FrgUjMba22r0qWTA/0J4IrQ68tJ6Z97YWa2idUgP+vuX2x1exrs7cDvmtkiq0NzO8xsprVNaooTwAl3L/y19nlWA39ajQA/cvfn3f1l4IvAP2xxm5rplJltBQh+nq73CTo50P9f4Eoze52ZXczqzZsHW9ymhjIzY3Xc9qi7f6LV7Wk0d7/L3S939yFWf7+PuXvqe3ru/rfAc2b2m8Gm64EftLBJjXYcuNbM+oJ/49eT4pvPER4Ebg6e3wx8qd4n6Ng1Y939l2b2AeBhVu/S3+/uT7W4WY32dmAc+L6ZPRFs+2N3f6iFbZLG+CAwG3RingXe3+L2NIy7P25mnwe+w2pm2XdJ6QxZM/sccB0wYGYngI8AHwMeMLNbWP3Se2/dz6uZsSIi6dbJQzciIpKAAr2ISMop0IuIpJwCvYhIyinQi4iknAK9iEjKKdCLiKScAr2ISMr9f+JO1b9/J5jIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#   plot my data\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure()\n",
    "plt.plot(x,y, 'ok')\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "plt.savefig('myData.png', dpi=300, type='png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   define a function to calculate the log likiehood\n",
    "def calcLogLikelihood(guess, true, n):\n",
    "    error = true-guess\n",
    "    sigma = np.std(error)\n",
    "    f = ((1.0/(2.0*math.pi*sigma*sigma))**(n/2))* \\\n",
    "        np.exp(-1*((np.dot(error.T,error))/(2*sigma*sigma)))\n",
    "    return np.log(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   define my function which will return the objective function to be minimized\n",
    "def myFunction(var):\n",
    "    yGuess = (var[2]*(x**2)) + (var[1]*x) + var[0]\n",
    "    f = calcLogLikelihood(yGuess, y, float(len(yGuess)))\n",
    "    return (-1*f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#    Let's pick some random starting points for the optimizaiton    \n",
    "nvar = 3\n",
    "var = np.zeros(nvar)\n",
    "var[0] = -1\n",
    "var[1] = 1.5\n",
    "var[2] = -1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 380.089609\n",
      "         Iterations: 14\n",
      "         Function evaluations: 130\n",
      "         Gradient evaluations: 26\n"
     ]
    }
   ],
   "source": [
    "#   let's maximize the liklihood (minimize -1*max(likelihood)\n",
    "#   In numerical optimization, the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm\n",
    "#   is an iterative method for solving unconstrained nonlinear optimization problems.\n",
    "\n",
    "from scipy.optimize import minimize\n",
    "res = minimize(myFunction, var, method='BFGS',\n",
    "                options={'disp': True})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 380.0896087488445\n",
       " hess_inv: array([[ 7.0125222 , -2.84080428,  0.23827251],\n",
       "       [-2.84080428,  1.71106895, -0.16597606],\n",
       "       [ 0.23827251, -0.16597606,  0.01730676]])\n",
       "      jac: array([0.00000000e+00, 0.00000000e+00, 7.62939453e-06])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 130\n",
       "      nit: 14\n",
       "     njev: 26\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([-1.41118235, -3.24773431,  4.00964233])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   pefrom least squres fit using scikitlearn\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.pipeline import Pipeline\n",
    "model = Pipeline([('poly', PolynomialFeatures(degree=2)),\n",
    "    ('linear', LinearRegression(fit_intercept=False))])\n",
    "\n",
    "model = model.fit(x[:, np.newaxis], y)\n",
    "coefs = model.named_steps['linear'].coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.09796836, -3.25735478,  3.95873605])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coefs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XlclNX+wPHPYRFFXBKVayqD3jQzMxWX1FJQXLLSbt3KQrO8XtosK1vscsvqSsu93tJ+ZWVqmpC23Gy1MhLCFjO3zCVzAzJNhUxBRAW+vz8YaIAZGJZhmJnv+/WalzPPc57nOQfqy5nznOd7jIiglFLKe/m5uwJKKaVcSwO9Ukp5OQ30Sinl5TTQK6WUl9NAr5RSXk4DvVJKeTkN9Eop5eU00CullJfTQK+UUl4uwN0VAGjdurVERETU6NgTJ07QtGnTuq1QA6dt9g3aZt9QmzZv2LAhS0TaVFWuQQT6iIgI1q9fX6NjU1NTiYqKqtsKNXDaZt+gbfYNtWmzMSbDmXI6dKOUUl5OA71SSnk5DfRKKeXlGsQYvT1nzpxh//795OfnV1quRYsW7Nixo55q1TA05DY3btyYDh06EBgY6O6qKKWsGmyg379/P82aNSMiIgJjjMNyOTk5NGvWrB5r5n4Ntc0iQnZ2Nvv376dTp07uro5SyqrBDt3k5+cTGhpaaZBXDYsxhtDQ0Cq/hSml6leDDfSABnkPpL8zpRqeBh3olVLKWyUlJfFEiw7cEx1NREQESUlJLruWBvpKhISEuOS8v//+O/PmzXO4PyEhgfPPP5+ePXvSq1cvvv32W5fUQynlHklJSbwyKZV5x9cynAgyMjKIi4tzWbDXQO8GlQX6b775hg8//JCNGzeyZcsWkpOT6dixY62vKSIUFRXV+jxKqdp7adprrC38P9ryI8/yMwB5eXnEx8e75Hoa6KvpyJEjXH311fTr149+/frx1VdfAbBu3ToGDRpE7969GTRoEDt37gRg27Zt9O/fn169etGzZ0927drFjBkz2LNnD7169eL+++8vc/6DBw/SunVrgoKCAGjdujVnn302AJ988gndunVj5MiR3HXXXVx++eUAPProo8yePbv0HD169CA9PZ309HTOO+88br/9dvr06cPPP//MqlWrGDhwIH369OGaa64hNzcXgBkzZtC9e3d69uzJfffd59ofolI+LGtnNvuyXyaUw2RwPUUUlu7LzMx0yTWrnF5pjGkMpAFB1vJvi8hMY8xiYChwzFr0JhHZbIrvxs0FxgB51u0ba1XLu++GzZvt7mpSWAj+/tU/Z69eMGdOtQ+bNm0a99xzDxdffDGZmZmMGjWKHTt20K1bN9LS0ggICCA5OZl//OMf/O9//+Oll15i2rRpxMbGcvr0aQoLC3nqqafYunUrm+20aeTIkTz++ON07dqVmJgYrrvuOoYOHUp+fj5///vfWb16NWFhYUyZMsWp+u7cuZNXX32VefPmkZWVxaxZs0hOTqZp06Y8/fTTPPPMM0ydOpUVK1bw448/Yozh999/r/bPRSlVtdeXLuXFm1qQxUg6cTEHyCqzPzw83CXXdWYe/SlgmIjkGmMCgS+NMR9b990vIm+XK38p0MX6GgC8aP3XKyQnJ7N9+/bSz8ePHycnJ4djx44xadIkdu3ahTGGM2fOADBw4EASEhLYv38/V111FV26dKn0/CEhIWzYsIE1a9aQkpLCddddx1NPPUWvXr3o1KkTXbp0IScnhwkTJjB//vwq62uxWLjooosAWLt2Ldu3b2fw4MEAnD59moEDB9K8eXMaN27MlClTuOyyy0q/KSilai8pKYn4+HgyMjJ4CniR87mLpaSwoUy54OBgEhISXFKHKgO9iAiQa/0YaH1JJYeMA16zHrfWGNPSGNNORA7WuJaV9LxP1vPDQ0VFRXzzzTc0adKkzPY777yT6OhoVqxYQXp6emk2uhtuuIEBAwbw0UcfMWrUKBYsWEDnzp0rvYa/vz9RUVFERUVxwQUXsGTJEnr16uVw6mJAQECZ8Xfbeey26U9FhBEjRrBs2bIK51i3bh2ff/45y5cv5/nnn2f16tVV/iyUUvbZBndjDCLCpXTmAfbyEttIYVuZ8haLhYSEBGJjY11SH6fG6I0x/saYzcBh4DMRKZkGkmCM2WKMedYYE2Td1h6sdxeK7bdu8wojR47k+eefL/1cMvxy7Ngx2rcvbubixYtL9+/du5fOnTtz1113MXbsWLZs2UKzZs3Iycmxe/6dO3eya9euMue3WCx069aNffv2sWfPHoAywToiIoKNG4tHxzZu3Mi+ffvsnvuiiy7iq6++Yvfu3UDxzZ+ffvqJ3Nxcjh07xpgxY5gzZ47dISWllHOSkpKIi4sjI6M4g7CIcA49+IItTOU+ppUrb4whPT3dZUEenEyBICKFQC9jTEtghTGmB/AQ8CvQCJgPPAg8Dtjrdlb4BmCMiQPiAMLCwkhNTS2zv0WLFg6Doa3CwkKnytVEXl5eafAGmDp1Kk888QTTp0+nR48eFBQUMHjwYObMmcMdd9zBrbfeyn/+8x+GDBmCiJCTk8Nrr73GG2+8QWBgIG3btuWee+6hUaNG9O/fn+7duzNixAhmzZpVeo1Dhw5x//33c+zYMQICAujcuTPPPfccZ86cYc6cOVx66aW0atWKQYMGcfToUXJychg5ciSLFi2iZ8+e9OnTh3POOaf0JmtRUVHpz6dx48bMmzePa6+9ltOnTwPw8MMPY4xh/PjxnDp1ChHhiSeeqNXPND8/v8Lvs7Zyc3Pr/JwNnbbZM02fPp28vLzSz81pyWlW0IxjvEUiZ8qVb926tevbLCLVegEzgfvKbYsCPrS+fxm43mbfTqBdZeeMjIyU8rZv315hmz3Hjx93qpw3OX78uKSkpMhll13m7qrY5ezvrjpSUlLq/JwNnbbZMxljhOLOrRj8JJKPJZBT0oOBpdtLXsHBwRIfH1/jawHrxYm4XeXQjTGmjbUnjzGmCRAD/GiMaWfdZoArga3WQ94HbjTFLgKOSW3G55VSyoPYzpwZQgIbGM0A7mAr3wB/pAmxWCzMnz+fmJgYl9fJmTH6dkCKMWYL8B3FY/QfAknGmB+AH4DWQMn4w0pgL7AbeAW4vc5rrYiKiuLDDz90dzWUUuUkJCQQHBzMtcC9fM1wZvOVWQgUB/elS5ciIi4fl7flzKybLUBvO9uHOSgvwB21r5pSSnmGklk2mZmZhIeHc//oK3lwxTLWywdkhH/P0ieW1ltQt6fB5qNXSilPUDLLpuQGbE7GCRZnPEa7pmdxy+5/sutPf3JzDTUFglJK1Up8fHxpkPcngI68xSHak9h4NzSAIA8a6JVSqlZs89MMYg7fE0VvpvDVb6vcWKuyNNBXwhjDxIkTSz8XFBTQpk2bOkkRkJqaavc8U6ZMKU2xUN00ySXlDxw4wF//+leg+OGtqVOn1rK2FbnqvEp5mpJZNhfzd9ZwB0P5N9/wusvy1tSEBvpKNG3alK1bt3Ly5EkAPvvsszIPULnCggUL6N69e63OcfbZZ/P22+VTECmlXCEhIYGRQUGMoSX9+Ig0HnJp3pqa0EBfhUsvvZSPPvoIKE47cP3115fuc5Sa+JlnnmHy5MkA/PDDD/To0aPMk3KViYqKYv369WW2ZWVlMXDgwNJ6zJ07l379+tGzZ09mzpxZ4Rzp6en06NGj9POBAwcYPXo0Xbp04YEHHijdvmzZMi644AJ69OjBgw8+WOX2V199la5duzJ06NDS9MxK+bobBg7ivcBArg14lp+4gnBLR+bPn+/WWTblecysG2uOsDKuvRYmToS8PBgzpuL+m24qfmVlgXUko5SzTxyPHz+exx9/nMsvv5wtW7YwefJk1qxZA+AwNfHdd99NVFQUK1asICEhgZdffpng4OBqtPYPhw4dYuzYscyaNYsRI0awatUq9uzZw7p16xARxo4dS1paGkOGDHF4js2bN7Np0yaCgoI499xzufPOO/H39+fBBx9kw4YNnHXWWYwcOZJ3332X/v37290+YMAAZs6cyYYNG2jRogXR0dH07l1h1q1SPiX311yu6HmER4giett/+b1rV3dXyS6PCfTu0rNnT9LT01m2bBljyv01cZSa2M/Pj8WLF9OzZ09uueWW0rTA1XXmzBmGDx/OCy+8wNChQwFYtWoVq1evLg2yubm57Nq1q9JAP3z4cFq0aAFA9+7dycjIIDs7m6ioKNq0aQNAbGwsaWlpGGPsbgfKbL/uuuv46aefatQupbxBUUERE3pvI+1EX8488Sg00CAPHhToHfXAc3IgOLjyHnrr1s734O0ZO3Ys9913H6mpqWRnZ5duf/jhh+2mJgbYtWsXISEhHDhwoMbXDQgIIDIykk8//bQ00IsI9957L9Omlc+B51jJalVQnAK5oKCgJA9RBY62Aw7TJCvli+IvSeO9X6OYe9UXjHxoqLurUykdo3fC5MmTeeSRR7jgggvKbHeUmvjYsWNMmzaNtLQ0srOza3xj1BjDokWL+PHHH3nqqacAGDVqFEuXLi3NTvnLL79w+PDhap97wIABfPHFF2RlZVFYWMiyZcsYOnRopdtL/tCdOXOGt956q0ZtUsobLL31S55aG0VctzTufMvxt+mGQgO9Ezp06GC3B/3AAw/w0EMPMXjwYAoL/1j38Z577uH222+na9euLFy4kBkzZtgNxp9//jkdOnQofX3zzTcVyvj7+7N8+XJSUlKYN28eI0eO5JprrmHgwIFccMEF/PWvf61RSuF27drx5JNPEh0dzYUXXkifPn0YN25cpdsfffRRBg4cSExMDH369Kn2NZXydElJSVzVrh2rXt5Nb5PKu4euwT/An4iICJKSktxdPcecSXHp6pemKa6eht5mTVNcN7TNDUtiYqKc27ix/AryE8hZNK2QcjgxMbHa561Nm6mrNMVKKeVrkpKSiIiIwM/Pr7S3PuvBBCz5SzhIJ64AjnKizDF5eXnEx8e7p8JV0ECvlFI2bJcCFBEyMjK4ZcothPzybz7nKibTmZ0OjrVNh9CQNOhAL5XMAFENk/7OlKezTVJWIjL/X6zncgYxlU187vDYhpT2wFaDDfSNGzcmOztbA4cHERGys7Np3Lixu6uiVI2V75Vfwq2kcQ9DmMOXZr7D4xpa2gNbDXYefYcOHdi/fz9HjhyptFx+fr7PBZaG3ObGjRvToUMHd1dDqRoLDw8nIyMDgBH48Ts30ZcPWcN0RARjDCJCaGgoAL/99hvh4eEkJCQ0qLQHtqoM9MaYxkAaEGQt/7aIzDTGdAKWA62AjcBEETltjAkCXgMigWzgOhFJr27FAgMD6dSpU5XlUlNTfe5RfF9ss1L1JSEhgbi4OCLy8niLInYwjBj8EIqA4m+uFouF9PR091a0GpwZujkFDBORC4FewGjrot9PA8+KSBfgKPA3a/m/AUdF5BzgWWs5pZTyCLGxscz95/P0MvM5TAjXkMcJcsuUaag3XR2pMtBbp2uWtDLQ+hJgGFDyyOcS4Err+3HWz1j3Dzf67LxSykPkZeUxP6E/78oNTA4dwH47ZRrqTVdHnLoZa4zxN8ZsBg4DnwF7gN9FpMBaZD9Qkqi9PfAzgHX/MSC0LiutlFKuUFRQxISeW1h/4jxef2grt869uULm2YZ809URU51ZLcaYlsAK4BHgVevwDMaYjsBKEbnAGLMNGCUi+6379gD9RSS73LnigDiAsLCwyOXLl9eoAbm5udVeicnTaZt9g7a5/i27JYf5P11B/MVvE/Ov1gAkJyezYMECDh8+TNu2bZkyZQoxMTF1ds3atDk6OnqDiPStsqAzj8/avoCZwP1AFhBg3TYQ+NT6/lNgoPV9gLWcqeyc9lIgOKshPzLtKtpm36Btrl+fXjdVwjgoQ3hOwjuG1yidQU00iBQIxpg21p48xpgmQAywA0gBSpbzmAS8Z33/vvUz1v2rrRVSSqkGKXX6dIa/8TxP0oc13E3mz5nExcU17ERl1eDMGH07IMUYswX4DvhMRD4EHgTuNcbspngMfqG1/EIg1Lr9XmBG3VdbKaXqxvrXtvPxM61Yj2EqB0unUTbk3DXVVeU8ehHZAlSYtC0ie4H+drbnA9fUSe2UUsqF9qZmctlNbWjEDSziv+RxtMx+T5tG6UiDTYGglFKulL3rNy4dWcgZAmjfZhJZ5YI8eN40Skc00CulfM7J304yNnI/GWfa0bPlJL49sqbCUpmeOI3SEQ30Simf8vrSpTx39jg25nShNxP54ugHAKV5bAAsFgvz589vsLlrqqvBJjVTSqm6lpSYyLHJk3mwoIBd/JmFHCyzXzwwj40ztEevlPIZSX//kZYF1zAbKgT5Et5yA9aW9uiVUj4h6fav+Dh/Fsd4na9Z5rCct9yAtaU9eqWU13v6uoXc9GJ/erGa77jZYTlvugFrSwO9Usqrzb3lFf715nV0Zht7+AtnOF1mv7fegLWlQzdKKe+1axfZr+zlLI5wlEvJ4XiZ3RaLpUGvDFVXNNArpbzTwYMwahR3yj6SeJH9HCuz2xjjdbNrHNGhG6WUx0lKSiIiIgI/Pz8iIiIqJB/7PeMY0Z138eW+VlwG7C0X5ME7b7o6oj16pZRHSUpKIi4ujry8PAAyMjKIi4sr3f/og48S/MtCdnARd9GKTXbO4a03XR3RQK+U8ijx8fGlQb5EXl4e06ZN49SJU3TPX8p3XMwArmctn1U43lfG5W1poFdKeRRHDzRlZ//GYBbyFVdyCVNZw5sVyvjSuLwtHaNXSnkUR2PrTxFIC9oylJms4YVqHevtNNArpTxKQkJChQW74wOCeZDTjGAcX/C43eN8bVzelgZ6pVSDkZycXDqbpnXr1rRu3brCzJrY2Fjmz5+PxWLBGMPo4Ol8XJDGxn6j+WeTRmXO5wsPQznDmTVjOxpjUowxO4wx24wx06zbHzXG/GKM2Wx9jbE55iFjzG5jzE5jzChXNkAp5R2SkpKYPXs2GRkZiAjZ2dlkZ2cjIqUza2yDfXp6Oom3reHTvH/TyO8QA75LpnFwMKGhoRhjsFgsLF26FBEhPT3dZ4M8ONejLwCmi8h5wEXAHcaY7tZ9z4pIL+trJYB133jgfGA0MM8Y4++CuiulvEh8fDynTp1yuL/8Gq4fPrKOG+cN4EK+YFPR1RRQQHZ2NidPnmTp0qU+H9xtVRnoReSgiGy0vs8BdgDtKzlkHLBcRE6JyD5gN3bWllVKKVvOpAcuKZM2dxPX/OsCuplN7GIsp8gvLeNNi3rXlWpNrzTGRFC8UPi3wGBgqjHmRmA9xb3+oxT/EVhrc9h+7PxhMMbEAXEAYWFhpKamVr/2QG5ubo2P9VTaZt/ga21u27Ythw4dqrLMhnnz6PDgHGKC/4+v8mI5QW6FcpmZmR7zs6uX37OIOPUCQoANwFXWz2GAP8XfChKARdbtLwATbI5bCFxd2bkjIyOlplJSUmp8rKfSNvsGX2tzYmKiBAUFCWD3FRgYKP2b9ZAj+Mm+gAB5+//+TywWi92yFovF3c1xWm1+z8B6cSJ+OzXrxhgTCPwPSBKRd6x/IA6JSKGIFAGv8MfwzH6go83hHYAD1fvzo5TyNbGxsdx3332ls2lCQ0NLb6yGhoYSIV3Zm5PMQ/yHoQUF3Pjgg4wZM6bCVEtfnkbpiDOzbgzFvfIdIvKMzfZ2NsX+Amy1vn8fGG+MCTLGdAK6AOvqrspKKW8VExNDeno6RUVFZGVlkZWVRVFREZ38u5BT8AkAybxMJsVj8StXriwz1dLXp1E64swY/WBgIvCDMWazdds/gOuNMb0o/qqUDtwCICLbjDFvAtspnrFzh4gU1nXFlVK+4edvD/Dr4dc5RRNCiWI3P5Xuy8zMJDY2VgN7FaoM9CLyJWDs7FpZyTEJFI/bK6VUjRUd+JUrhhzjOGdzNsP4sXTgoJivpjSoLn0yVinlVra55cePH/9HbvkjR/AbMZw5fvfy6M1vkRn8Y5njdCzeeZq9UinlNuVzyx86dIi4uDhOHTxNyP9t59oj+4j6+AWioqJoO7wJ8fHxZGZmEh4e7nOphmtDA71Sym3s5ZYPyAtk9gMXsk/GM/j1K2gfNQRAx+JrQYdulFJuU/5p2GY0pwOr2C09+N+jW2l//RA31cy7aKBXSrmN7c3UEJoRzif8RG+GNv8bY2b2c2PNvIsGeqWU25Tklg8B7uUKdtKXAQE3cNO80e6umlfRQK+Uqhe2s2tKcs1PnDiR1o0bsyoggId5nb+E9Oe2xVfqWHwd05uxSimXKz+7Jjs7G4BgmtL+tyUUMotvpvbj9quvJioqyo019U7ao1dKuZy92TXBNOUcPmYdo7mfjkz84AM31c77aY9eKeVy5WfXNCWEP7OSbQykH7Gs5W1Mpr0H8FVd0B69Usrlys6uCaETn1qD/PWs5c0KZVTd0kCvlKpTtjddSxb1Lpld0xxYST7nkUkk17KWtwFNZ+BqGuiVUnWm5KZryQLfJYt6Azw/6wXeCezAAAoIahbHntA0TS1cT3SMXilVa0lJScTHx5ORkVFhX15eHk/d/zQBv7+Jf8G7rHv3AEvHXeGGWvouDfRKqVopP3WyvNa0Jf9gEvux8O6sbfhpkK93GuiVUrVib+pkiTDOpinJHCCcj2b/yLDpfeu5dgqcW0qwozEmxRizwxizzRgzzbq9lTHmM2PMLuu/Z1m3G2PMc8aY3caYLcaYPq5uhFLKfcpPnSxhAbrzEkdoz+M3LmPY9N71WzFVypmbsQXAdBE5D7gIuMMY0x2YAXwuIl2Az62fAS6leJ3YLkAc8GKd11op1WDYmxZ5DpAGvGhuZVbcG0xfMqXe66X+UGWgF5GDIrLR+j4H2AG0B8YBS6zFlgBXWt+PA16TYmuBluUWEldKeSh7+WoyMjIw5o+Hnc7hfAbzAq1CWnLuxo+46+W/u7HGCqo5Rm+MiQB6A98CYSJyEIr/GBhj2lqLtQd+tjlsv3XbwdpWVinlPo7y1QCICMYYukkkv/IJq8wpjr43gpBeXdxVXWXDiIhzBY0JAb4AEkTkHWPM7yLS0mb/URE5yxjzEfCkdVFxjDGfAw+IyIZy54ujeGiHsLCwyOXLl9eoAbm5uYSEhNToWE+lbfYNDa3N48eP59ChQw739+QS9vEhof5HmT17M6G9WpCcnMyCBQs4fPgwbdu2ZcqUKcTExDg8R0Nrc32oTZujo6M3iEjVd7hFpMoXEAh8Ctxrs20n0M76vh2w0/r+ZeB6e+UcvSIjI6WmUlJSanysp9I2+4aG1mZjjAB2X30ZLY3Jk85sl/3fHRARkcTERAkODi5TLjg4WBITEx1eo6G1uT7Ups3AenEihjsz68YAC4EdIvKMza73gUnW95OA92y232idfXMRcEysQzxKKc/lKBfNdcAz/E531hHcbjzt+xbfkrM37TIvL4/4+HhXV1WV48ysm8HARGCYMWaz9TUGeAoYYYzZBYywfgZYCewFdgOvALfXfbWVUvWtJF+NravoxetAEWvZ3/hSZvzngdJ9jqZdOtquXKfKm7FSPNbuKH/ocDvlBbijlvVSSjUwJbloSlIdxAQ8wjsFj/Eol/NGxy088+STZfLVhIeH202JoFkq658mNVNKOS02Npa9u/dyd+9UkgseY0KnL3n4xAp2ZmZWSEpm7xuAZql0Dw30Simnnck7w03nfs2cTUOZ1usLlvw0iMDgQLtlY2NjmT9/PhaLRbNUupkGeqVUpUoekgoxhmdaXs3SvRczKyaVZzcMwS+g8hASGxtLeno6RUVFpKena5B3Ew30SqkKSoK7MYaJEydyPCOTZOC+Mx9wbUAkETf9gvHTpf88hQZ6pVQZtouHALSTjvyJ9eRyMVcDbxZs1CmSHkYDvVI+xN4yf+XZzn/vQk8K+YYDdOJe/nhYRqdIehbNR6+Ujyifq8Z2mT/bsfOSIN6LaPaygqYcJ5SL+YHtpWV0iqRn0R69Uj7C2SdVw8PDGUkk2/iEtmQiDGSvTZDXKZKeRwO9Uj7CqSdVRVgxYAAfs5GJPMZhLuFXfilNQ6xTJD2TBnqlfISj4ZaS7QX5BcT3/ZTQN78hY9BA1oQnkmOOY7FYWLp0KSKiUyQ9lI7RK+UjEhISKiziXTIMk/trLuMv3M5Hh0fTdvgLTFt1GT/5aT/QW+hvUikf4ehJ1ejzhjOk8898fDiSF69PY1ryFaBB3qvob1MpH1L+SdX+LfpyUf9Cdp3swIePbeTW14e4u4rKBTTQK+WrPv6YP10fTY+AH1mz/ACXPtLP3TVSLqKBXikftOzmVZy47FqadfkTK/d0o9d157q7SsqFNNAr5SWceeq1IL+Auy78ghsWj+SJdrPonpWFX8eODssr76CzbpTyAs489Xos8xjX9dnFp9lDmRTxDnN/jedE/gmH5ZX3cGbN2EXGmMPGmK022x41xvxSbmnBkn0PGWN2G2N2GmNGuariSqk/VPXU697UTAZ1PcLn2Rfyyo1rSJV7S4O8bflJkyZV+o1AeSZnhm4WA6PtbH9WRHpZXysBjDHdgfHA+dZj5hlj/Ouqskop+yp76jU5Pp5jw4Zx6lQhl5w1niYjMx2WLywsRERKe/ga7L1DlYFeRNKA35w83zhguYicEpF9FC8Q3r8W9VNKOcHRU69XNbqCIU88QZDsQehOytF3iIuLo1WrVlWe014eHOWZTPFa3lUUMiYC+FBEelg/PwrcBBwH1gPTReSoMeZ5YK2IJFrLLQQ+FpG37ZwzDogDCAsLi1y+fHmNGpCbm0tISEiNjvVU2mbfUJ02JycnM3v2bE6dOgWAPwEMZi5p3M50YnmF1zluU7558+acOnWqtLwjxhhWr15d0yZUm/6eqyc6OnqDiPStsqCIVPkCIoCtNp/DAH+KvxEkAIus218AJtiUWwhcXdX5IyMjpaZSUlJqfKyn0jb7Bts2JyYmisViEWOMWCwWSUxMrFC+pEworaW3SRMQGcrTYvAToMzLGFPmnP7+/hXKAGKxWOqvwaK/5+oC1osTMbxG0ytF5JCIFIpIEfAKfwzP7Ac62hTtAByoyTWUUsVsV3ySSsbPY2NjeffpT2nqv4kd0pfE274i3TIPoajCOcPDw8s8JbuMaZ/2AAAZD0lEQVRkyRKCg4PLlNF0xN6jRoHeGNPO5uNfgJIZOe8D440xQcaYTkAXYF3tqqiUb3M2jzyvv87vN96FH0WsWbKP2HmDSUhIcCqAO8qDo1MtvYMz0yuXAd8A5xpj9htj/gb82xjzgzFmCxAN3AMgItuAN4HtwCfAHSJS6LLaK+UDKptRk5SUROfwztxiRkJsLN07ZbBzXxA7/TcRERHBxIkTadKkCaGhoVUG8PJ5cDTIe48qH5gSkevtbF5YSfkEisftlVJ1IDw8vHShblutWrVixpR/0Dr/VeYzjLO4gBcz9xD75GMsWbKk9FtAdnY2wcHBLF26VIO3j9IUCEo1cI6GXzrnX0Bh/pf8yEAGM4mn2crxkyeZP3++c0M9ymdooFeqgaswfh4ezrQ+L/D9iU/wp4BwBvEVr5WWLyy0P1rqaAhIeT8N9Eo1UMnJyaVJyuLj40lISKAoN5f06Gg6fvktF/qlkUtffmJzmeP8/e0/jO7ooSrl/TTQK9UAJSUlMXv27DJTKh+b/ARvRtwGr73GrQ+35a5Fv3I6OL/MccHBwcTFxelUSVWGBnqlGqD4+PgyT60O4CoOnf6a+448xun3PsY8/hgTJk20OyVy3rx5OlVSlaFpipVqgErG0wMIZBD/Jo276c63HONaGl3xxwyc2NhYh1MlNbCrEtqjV8rN7C0YEh4eTkeacC5ppHE3Q5jDT1xCgMW4u7rKA2mgV8qNHKU3eKhbN77nJNF8yQD+Shr30Cg4UMfZVY1ooFfKjcqnNwikEf3z/sWAT39ln58fX7dMYJ15R8fZVa3oGL1SbmQ7tz2CrgSxnFR68zC/sqroewJOn9YnWlWtaY9eqXriaCweYDA3coQNHKEj/biCD/kPp9EnWlXd0B69UvXA0eLdt4wfT9GrecyVJVxIKr8Sy3flMnvrE62qtrRHr1Q9sJdquEeecM+SJTzNO1zd5C62MJxDdpZvKP9Eq71vBkpVRgO9Uk6qTYC17ZX7E8BQHuVXdnC4sAVBX3/B23nPsTTxtSqfaHV2ERKlbGmgV8oJNQ2wJX8cxLo2czh/phtr+IKZdCSVaM4QMX48SUlJFZKXhYWFVZhp4/QiJErZ0DF6pZxQWYB1NCOm/Lj8xfydTTzDcU5zEdfyFW8BkJORQ1xcHFD2idbU1FSioqLKnLOyRUiUcsSZFaYWGWMOG2O22mxrZYz5zBizy/rvWdbtxhjznDFmtzFmizGmjysrr1R9cTbA2g7vTJo0iby8PMIoXmOzGVfSlbU0pidrrUG+hLO9ckcZKDUzpaqMM0M3i4HR5bbNAD4XkS7A59bPAJdSvE5sFyAOeLFuqqmUe1UWYEuCuzGGiRMnlg7vFBYWMoCrWYmFGKAj17GJkfzKL3bP5Uyv3Nk1YJWyVWWgF5E04Ldym8cBS6zvlwBX2mx/TYqtBVqWW0hcKY/kKMCOGTOmdOweKB2Lb0krBvI63/I2T3MvfYD55AJSq3zxuoi3qoma3owNE5GDANZ/21q3twd+tim337pNKY/mKMCuXLmywth9Py6nEVtZz9UMJZ63uZcfrfvqIl+8LuKtqquub8baS60ndgsaE0fx8A5hYWGkpqbW6IK5ubk1PtZTaZvdo3379ixevLjMNtvhlrOAvzKJV1hMV76nJZfyBd/j5+eHEUPbtm2ZMmUKMTExtGrVigULFnD48OHS7e3bty/TxobQ5vqmbXYREanyBUQAW20+7wTaWd+3A3Za378MXG+vXGWvyMhIqamUlJQaH+uptM0Nh8ViEUDG0FwOgGTRTGK4TwIIFECCg4MlMTGxRuduqG12JW1z9QDrxYkYXtOhm/eBSdb3k4D3bLbfaJ19cxFwTKxDPEp5o4dvnclA8ybpfM3PNGI4OXxu/ksBZ3T8XDUYVQ7dGGOWAVFAa2PMfmAm8BTwpjHmb0AmcI21+EpgDLAbyANudkGdlXI7KRKW3voVDywYR640JbrxvxmcX0R7i4WlCQka3FWDUmWgF5HrHewabqesAHfUtlJKNRRJSUnEx8eTmZlJeHg4CQkJjO42mKuGHiDtxMX04GuatnmIic/G8UnsGXdXVym79MlYpcopCe4ZGRkYY0qnTP6SkcH2m2/mL0V+FBV+wiXczhpegiPCD3HrAbQnrxokzXWjlA3bnDbwx7z47gygL59y35mmJBeeYi/RrOFFSiaVab4Z1ZBpoFfKRvmcNs1pySW8wA6+JoPz+AudGQd2kglrvhnVcGmgV4o/ctSU9OQBBjGRQHbyFbdwCXM5Tne+YGOtnmxVyh000CufV364pgeQCoTyV/7EHs4hkjTu5QS5dfJkq1L1TQO98ho1XRikZLimGc0ZyjO8yTmcD5zFRLYxmF1mC0DpvPh58+ZpvhnlUXTWjfIKjtZkhapnwmRmZDKYm9nJk6TRhsf5iVXs5jeOY7FYSLAzL942b7xSDZ326JVXcHblpfK9/jlTFnCeWcdXLOJP7OFc+rGcl/iN4h68Jg1T3kB79MorOLMwiG2vvwOQkJHB+oU5/GbaM8R/EmmFr5WW1TF35U20R6+8gjMrL8XHx0OeYSiP8gpDuQpoykxanB1N3JKROuauvJYGeuXRbKdFGlM2S7Ztr7zwdCEdM2IIYRdfMJP5DKcbkEAOO3/ZSXx8PAkJCZrjXXklDfTKo9iOsbdu3ZrJkyeXeYq1JNiX9spvuIHkhG/p3WIPX7KAtuzjfAaygkewHewpuXnr7EwdpTyJBnrlMWznu4sI2dnZnD59ukwZEfnjJuq558Lw4Wz/ZxJ5BUE8NnoRexrHsI21ds+vaQyUt9JArzyGvZk19vhlBHFNh29Y2O9F+OEHbpvTje1H2/HIx5N5ZcErWCwWh8dqGgPljTTQK49RVRBuR0cuZj6ZbOOTX3qQN/pq2LOHN1u3oGuPrvj5+ZWOxTsK9prGQHkjDfTKYzgKwmHAldxDNrtYx43cEP4Be7bmc+fHY0j64IMywz0lY/FjxozRNAbKZ2igVx4jISGhTHAOpQ3/Mi3ZC0xiF4MDXuf/ZibxWsZf+GzzKiIiIpgwYYLdB6lWrlypaQyUz6jVA1PGmHQgBygECkSkrzGmFfAGxQuKpwPXisjR2lVT+ZryKztNmDCBWbNmAfD0/f+h1cEb+I47OCXPEDxxF1c+8ghXnnNO6bG26RDsyczM1DQGymfURY8+WkR6iUhf6+cZwOci0gX43PpZ+bjqJBwrP7smIyOD2bNnszBhERv+257dB79mDdP5S8RmJqyMhddeA2uQB+du2upYvPIlrkiBMI7ixcQBllCc8fVBF1xHeYjqJhwrH6gjgPtPneKDf7bkQy4m9s9r+cfz7Tl39GC716vqpq2OxStfU9sevQCrjDEbjDFx1m1hInIQwPpv21peQ3k4ZxOOlSgJ1J3pzmAW8zF/ZgrQlxn0ajuEkY9lcO7oTg6vV1lvXcfilS8yJWti1uhgY84WkQPGmLbAZ8CdwPsi0tKmzFEROcvOsXFAHEBYWFjk8uXLa1SH3NxcQkJCanSsp/K0Ng8bNgx7/50ZY1i9enWF7fFjn+ZMzu18xxUEc4KruZnVvMUv1v1BQUHcd999xMTE2L1ecnIys2fP5tSpU6XbqjqmIfK033Nd0DZXT3R09AabYXPHRKROXsCjwH3ATqCddVs7YGdVx0ZGRkpNpaSk1PhYT+VpbbZYLELxt78yL4vF8kehggIpfOMtiWq+XkAklMMylIelJa2qPtaOxMREsVgsYowRi8UiiYmJLm2jK3ja77kuaJurB1gvTsTnGg/dGGOaGmOalbwHRgJbgfeBSdZik4D3anoN5R3KT4uEP8bJj+8/zus3fgJduuB33TVE+X/J89d8wX/nrSTd8hq/85vdc1Y1Dh8bG0t6eromKVOK2o3RhwFfGmO+B9YBH4nIJ8BTwAhjzC5ghPWz8jLlk4u1bt3a4Yya2NjYCnPW/zvtOb6b3YEOHSF26Wh+aD4Y3nmHmUemcsebQ5l02yTS09P1CVal6oIz3X5Xv3Topnrc3ebExEQJDg62O6QCSHBwsP2hksJCyVi8Wsa0WSeGQgngtIxr/YlcEjrK4RCLvWs5PL+Xcffv2R20zdWDk0M3usKUqraq5qnbzqiJj48nJyOX8SE9eCp4B6GHc9nnv4lHhqbRYcQ+pj0xtdJplyX/Tp8+ncOHDxMeHm53DVellGOaAkFVmzMZHjMyMphz82LaZ8ziBPtJyZ3HpiNZbLr9Jv654BsWp9/E3/852e60ywkTJpQZAoqNjWX58uU63q5UDWmPXlVbeHh46WIf5bUFBnAVW0lg/ZluNOd3+rOAX3iJoVJE6BtvcHLx4iqfXK3qoSqllPO0R6+qrfwsmgAC6cdYFhHGL8ANBBJCFoO4kdOczRruZC/bAMjOznYqpzzoQiBK1RUN9KqCqmbUxMbG8vJLL3Nx6CiGMJcW/MJ3vEc6sSxo3pzgpyM5bpnA1ywln5O1qosuBKJU7enQjSqjfF6a7Ozs0n0ZGRk8OWUK5/7vA5786F9sPz2BIPIZ13EjN05OZ9SMpwlo/F8ActqfXSGDZHBwME2aNClzzqroNEqlak979B6sOhkhnWVvRk07OjCUe7mG59ian0/fFW9wRdh3zJ+YxsG9+byROYjLHu1HQOM/+g325s7Pnz+fuXPn2n146rbbbtOFQJRyFWfmYLr65Wnz6N39eH1KSorL5pcbYwSQ9nSQIdwr5/O1gAiInM9auZtA6Qi1uoajn19lP1edX+0btM3Vg5Pz6N0e5MXDAn1lAdY2UIWGhkpoaGid/DEoHwDj4+Odyx9j51hHQXXpa0vlh3d+klnNO8gGkGe4W0CkG+tlKA9JOH92Ks+Mq/4IagDwDdrm6tFAXwfsBS1HATY0NLRmT4s6UYfy5w0KCnJ4HWNMpccGBwfLbbfdJsHBwRJAoPRmuAxhjnRkj4DIEibKN35+MpVW0pFOds9fEuzr8ylWDQC+QdtcPRroa8lR0HIUYJ15VZVx0R5Hf1j8/f2rvIa9Y9sTKFNAXqW1hHBcQKQxedKP9+Vi/i4X/qmn3HbbbRW+mdgGeUdB3NlvGTWhAcA3aJurx9lA75U3Y+viJqWjxTL8/f1rXK+qpgraq7ejYwoLCx3evCw5T0ZGBiE0oy+XMYRn6cx2BpHIK0A0WVzCHPoxFgjlO8byJa/w/a9bWLJkCQkJCRQVFZGVlUVWVhYWi6W4Z1Du52E7z91RXXWKpFJu5sxfA1e/6rJHX1fDB+V7r+XPV/5zSa+3sldVY9vVOa/tUJLt0FLiS4vkiqAgmQUyhkXiz5nSXnskH0sMU6R7Db59OPp52A4VaY++bmmbfYMO3Tih/A+proJNZeexG2BrmtGxiuvZG/sPCgoqPdehrYdlxYy1Mj0yRQY0/UFa8pucJkDOgEzndhnKLOlFtDTC8bh+VQHc2Z+rjtHXLW2zb9BA74TyPyRnep7OqEnQqs2sm8q+QZScK4gg6R8cLXcFhcpSkEf87hEo/i0GcVIuab5ZhpIg0TSTptUI6s78YXT256GzbuqOttk3aKB3gqt69CL1O1++fL1DQPrQVi7hNrmYV6QbGySA0wIinzJCfgH5DxfKML8H5ZmbXpH8Y/l2z1PZNwPbGTjO/EFz5/MDGgB8g7a5etwe6IHRFK8fuxuYUVnZ2gb68j3pRo0a1bgnXt0AVtvgl5edJ9+/tVNmjlgo0SZB+vOOzONKEZBNXCgg0ooj0ptVEsWT0p+/yNmcVaPhk5o8qNRQaADwDdrm6nFroAf8gT1AZ6AR8D3Q3VH52gT6+Pj4CoEtMDDQ6WGT2owrO3ts7qFc2fruLvng4W/l2StT5cPLXxQZNUqyO14ohsLS4Rd/zkhntsvt3CjxIJfjL+3oWO3xdE8I3NWlAcA3aJurx9lA76qkZv2B3SKyF8AYsxwYB2yv6wstWLCgwjTIM2fOEBISQlZWVpXHO5pGOWHCBOLj4ytdzeiRf/yDpnlNOZs/04wwgjmbgLx2pE5JI/att+Dnn/nz5v+xtygCOKf0uJsDf+aynlm0uuR8njicQqfuTeg+pDVdR1gIan4eERFfUJzuvRD4ufQ4f39/CgsLK9SjfOKv2NhYzeGulCrlqkDfHtsIBfuBAXV+lY0bGXtoMD8xAMhFyEU4iXCC9hlLYcUKvt/fil+PN8U/0I9Nmzfy8ScrOXb0V8KbF+dHb3f8HNpwAX4E4E8QfjShCSe5kDdompHBihu/5YMHmnPiVGMOHRVOFLWgs0nn7eC/sevECXryA9voUaZajfJXwu4HIDycCb220qRNOhFdGxHRsznnXNKO0C6x4DcBgBl2mpWQkGA38+OkSZNYsmRJhe2a+EspVRlT3Puv45Macw0wSkSmWD9PBPqLyJ02ZeKAOICwsLDI5cuXV/s6oV9+yQsPN2UJt5JPk9LtwZzgBCEATGApSUwoc1xbDnGIPwEwjnd5n3Fl9ndmD3usPfDLeI/v6E9TjhHMMZpwlD+xhcv8ZnAsIIAvTl/KUQI4yWF+5wBHOEjzsGZU1p7k5GQWLFjA4cOHadasGQA5OTm0bduWKVOmEBMTU6ZM+e3z588nKyurzHZvl5ubS0hIiLurUa+0zb6hNm2Ojo7eICJ9qyzozPhOdV/AQOBTm88PAQ85Kl8XY/R++EszmktrwuScoK7y/tNPy4cJCRLZ6DzpwUC5gMHSk0vkQqKkF4NkCMglIL3oIucRKV25UDrRTdpjkTa0kRAQvyrGxh3NZKnuPYHqHC9if0zPG8flbenYrW/QNlcPbr4ZGwDsBTrxx83Y8x2Vr8tZN7ZBztFUw7p6GWOqHWCdqVNVU0Fd9TRwQ6YBwDdom6vH2UDvkjF6ESkwxkwFPqV4Bs4iEdnmimuB45uP1c2xEhoaysmTJ51eFSk8PLzaNz6dqVN16+3ohnJ8fLzelFVKuS6pmYisFJGuIvJnEXHL3cLqLEMXHBzM3Llzq7UqUk1ugjpTp+oun6fJxJRSlfHY7JUlGRqHDRvmMENlQkJChQBtjAGKe++hoaFlAnpJ7zw9PZ2ioiLS09NLt9n7A1CT3rK9OtmqyR8QR38YdL1VpRTgmSkQqjMm3RBvUtZ2JSodo/cN2mbf4LE3Y6v7qm6gd2U6XE+gs258g7bZN3jszVhX0zHpivRpWKWUIx45Ru9pY9J1seKVUkrVlEcGens3NBtqKoCkpCTi4uLIyMhARMjIyCAuLk6DvVKq3nhkoK/LWTCuVtkcd6WUqg8eOUYPf4xJp6amEhUV5e7qOKT3E5RS7uaRPXpP4mn3E5RS3kcDvYt50v0EpZR30kDvYp50P0Ep5Z08dozek+gcd6WUO2mPXimlvJwGeqWU8nIa6JVSystpoFdKKS+ngV4ppbycKc506eZKGHMEyKjh4a2BrDqsjifQNvsGbbNvqE2bLSLSpqpCDSLQ14YxZr2I9HV3PeqTttk3aJt9Q320WYdulFLKy2mgV0opL+cNgX6+uyvgBtpm36Bt9g0ub7PHj9ErpZSqnDf06JVSSlXCowO9MWa0MWanMWa3MWaGu+vjasaYjsaYFGPMDmPMNmPMNHfXqT4YY/yNMZuMMR+6uy71xRjT0hjztjHmR+vve6C76+RKxph7rP9NbzXGLDPGNHZ3nVzBGLPIGHPYGLPVZlsrY8xnxphd1n/PquvremygN8b4Ay8AlwLdgeuNMd3dWyuXKwCmi8h5wEXAHT7QZoBpwA53V6KezQU+EZFuwIV4cfuNMe2Bu4C+ItID8AfGu7dWLrMYGF1u2wzgcxHpAnxu/VynPDbQA/2B3SKyV0ROA8uBcW6uk0uJyEER2Wh9n0Px//zt3Vsr1zLGdAAuAxa4uy71xRjTHBgCLAQQkdMi8rt7a+VyAUATY0wAEAwccHN9XEJE0oDfym0eByyxvl8CXFnX1/XkQN8e+Nnm8368POjZMsZEAL2Bb91bE5ebAzwAFLm7IvWoM3AEeNU6ZLXAGNPU3ZVyFRH5BZgNZAIHgWMissq9tapXYSJyEIo7c0Dbur6AJwd6Y2ebT0whMsaEAP8D7haR4+6uj6sYYy4HDovIBnfXpZ4FAH2AF0WkN3ACF3ydbyisY9LjgE7A2UBTY8wE99bKu3hyoN8PdLT53AEv/bpnyxgTSHGQTxKRd9xdHxcbDIw1xqRTPDQ3zBiT6N4q1Yv9wH4RKfm29jbFgd9bxQD7ROSIiJwB3gEGublO9emQMaYdgPXfw3V9AU8O9N8BXYwxnYwxjSi+efO+m+vkUsYYQ/G47Q4Recbd9XE1EXlIRDqISATFv9/VIuL1PT0R+RX42RhzrnXTcGC7G6vkapnARcaYYOt/48Px4pvPdrwPTLK+nwS8V9cX8Ng1Y0WkwBgzFfiU4rv0i0Rkm5ur5WqDgYnAD8aYzdZt/xCRlW6sk3KNO4EkaydmL3Czm+vjMiLyrTHmbWAjxTPLNuGlT8gaY5YBUUBrY8x+YCbwFPCmMeZvFP/Ru6bOr6tPxiqllHfz5KEbpZRSTtBAr5RSXk4DvVJKeTkN9Eop5eU00CullJfTQK+UUl5OA71SSnk5DfRKKeXl/h94MrICMpAaUwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#   plot the data and model fits\n",
    "plt.figure()\n",
    "plt.plot(x,y, 'ok')\n",
    "plt.plot(x,model.predict(x[:,np.newaxis]), '-r', label='Least Squres')\n",
    "plt.plot(x,(res.x[2]*(x**2)) + (res.x[1]*x) + res.x[0], '--b', label='Max Likelihood')\n",
    "plt.grid(True)\n",
    "plt.legend(loc=2)\n",
    "plt.show()\n",
    "plt.savefig('maxLikelihoodComp.png', dpi=300, type='png')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}