Gauss_Newton_for_Parameter_Estimation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Gauss-Newton for Parameter Estimation.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyO5hKtL1xNowQluSSNEyGQH",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/pkhamchuai/Optimization/blob/main/Gauss_Newton_for_Parameter_Estimation.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Gauss-Newton method for parameter estimation\n"
      ],
      "metadata": {
        "id": "fRNXXitx2tv_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "The task is to estimate the parameters $p_1$ and $p_2$ of the function\n",
        "$f(x) = \\frac{p_1x}{p_2+x}$ using Gauss-Newton method. \n",
        "\n",
        "You can use this code to perform parameter estimation of your interested function. For that you'll need to know the Jacobian matrix, aka, the derivatives of $f(x)$ wrt to $p_1$ and $p_2$."
      ],
      "metadata": {
        "id": "QklQX6tiUNVh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline"
      ],
      "metadata": {
        "id": "KXUIcuLFr60r"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "VtD4t0_nqzTD"
      },
      "outputs": [],
      "source": [
        "# sample data\n",
        "X = np.array([0.200, 0.200, 0.222, 0.222, 0.286, 0.286, 0.400, 0.400, 0.667, 0.667, 2.000, 2.000])\n",
        "Y = np.array([0.0831, 0.0781, 0.0948, 0.0891, 0.1163, 0.1213, 0.1467, 0.1663, 0.2586, 0.2260, 0.5557, 0.4207])"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(10,6))\n",
        "plt.scatter(X, Y, c='red', marker='o', label=\"Data\")\n",
        "plt.legend()\n",
        "plt.xlabel(\"x\")\n",
        "plt.ylabel(\"y\")\n",
        "plt.legend(fontsize=14)\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 388
        },
        "id": "TqY6Fv5TqUMZ",
        "outputId": "ef54ef07-c4a9-4851-a4fe-eed608dc63b4"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 720x432 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFzCAYAAAB2A95GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAdMElEQVR4nO3dfZCdZ3kf4N8tWzIVssGOLMHYlgSMZ8B1U8AC0gSIRNxgUmynxQyGDYmbMgoZnEAZ0iEoQAwRbcmMO+ngpmwSAk0WBIGEKuAYQsiGaVPAH2CwMR/GIykWJtgmJMgL+OvpH2cFK+UcS3J09jnava6ZnbPv8746e59b7579zfvxnGqtBQCAxbWidwEAAMuREAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQwYm9Czhaa9eubZs2bepdxkj33HNPHvnIR/YuY+Loy3D6MpreDKcvo+nNcPoy2mL05vrrr7+rtXb6sHXHXQjbtGlTrrvuut5ljDQ7O5stW7b0LmPi6Mtw+jKa3gynL6PpzXD6Mtpi9Kaq9oxa53QkAEAHQhgAQAdCGABAB0IYAEAHQhgAQAfH3d2RD+XBBx/MXXfdlW9961t54IEHutTwqEc9KrfcckuXn32snHDCCXn0ox+dtWvXZsUKOR0AxmFJhbDbb789VZVNmzZl5cqVqapFr+Hb3/52Tj755EX/ucdKay333Xdf/vZv/za33357NmzY0LskAFiSltRhjnvuuSdnnHFGVq1a1SWALQVVlVWrVuWMM87IPffc07scAFiyllQIS+L02TGijwAwXv7SAgDLy8xMsmlTcv31g8eZmS5lLKlrwgAAHtLMTLJtWzI3N1jes2ewnCRTU4taiiNhAMDysX37DwLYAXNzg/FFJoRNgMsuuyxVlarKypUrs27dumzdujVXXXVV7rvvviN+ntnZ2VRV7rrrrjFWCwDHsb17j258jISwCXH++efnjjvuyO7du/PRj340F154Yd74xjfmWc96lrsUAeBYGTX1UocpmYSwCXHSSSflMY95TM4444w8+clPzqtf/erMzs7mhhtuyFvf+tYkyR/+4R/maU97Wk4++eSsW7cuL3zhC7Nv374kye7du7N169Ykyemnn56qymWXXZYkueaaa/KsZz0rp556ak477bQ897nPPe4nlAWAh2XHjmT16oPHVq8ejC8yIWyYA3dNrFjR9a6Jc889NxdccEE+8IEPJEnuvffeXHHFFbnxxhvzoQ99KHfddVde/OIXJ0nOOuus7293880354477shv/dZvJRnMn/aqV70qn/70pzM7O5tHPepRufDCC3Pvvfd2eV0A0M3UVDI9nWzcOFjeuHGwvMgX5SfujvzHJuiuiSQ555xz8rGPfSxJ8vM///PfH3/84x+f3/7t386TnvSk3H777TnzzDNz2mmnJUnWrVuXtWvXfn/bF7zgBQc95+///u/nlFNOyac//ek885nPXIRXAQATZGpq8DU7m+ze3a0MR8IONUF3TSSDjxE6MPv/DTfckIsvvjgbN27MySefnM2bNydJ9h7mYsKvfvWreclLXpInPOEJOeWUU7J+/fo8+OCDh/13AMD4CGGHmqC7JpLkC1/4Qh7/+MfnnnvuyXOf+9ysXr06f/AHf5Brr70211xzTZIc9rTi85///Nx55515+9vfnk996lP5zGc+kxNPPNHpSADoSAg71ATdNXHTTTflmmuuySWXXJIvfvGLueuuu/KWt7wlz372s/PEJz4x3/jGNw7aftWqVUmSBx544Ptjd999d774xS/mda97Xc4///w86UlPyre//e3cf//9i/paAICDCWGH6nTXxPe+9718/etfz9e+9rXceOONufLKK7Nly5acd955ec1rXpMNGzbkpJNOytve9rbcdttt+fCHP5zXv/71Bz3Hxo0bU1X58Ic/nDvvvDP79+/PqaeemrVr1+Z3fud3cuutt+av/uqv8vKXvzwnnuhyQADoSQg71MK7JqoW7a6Jj33sY3nsYx+bDRs25Cd+4ieya9eu/Pqv/3o+8YlP5JGPfGROP/30vOtd78oHP/jBnHPOObniiity5ZVXHvQcZ5xxRq644ops374969evz+WXX54VK1bkve99bz73uc/l3HPPzSte8Yq8+c1vzkknnTTW1wMAPDSHQ4Y5cNfEInnnO9+Zd77znYfd7kUvelFe9KIXHTTWWjto+fWvf/0/OkL2nOc8JzfddNNBY/v37394xQIAx4QjYQAAHQhhAAAdCGEAAB0IYQAAHSy5EHboheo8PPoIAOO1pELYypUr853vfKd3GUvCd77znaxcubJ3GQCwZC2pELZu3brs27cvc3NzjuQ8TK21zM3NZd++fVm3bl3vcgBgyVpS84SdcsopSZKvfe1rue+++7rU8N3vfjePeMQjuvzsY2XlypVZv3799/sJABx7SyqEJYMg1jM8zM7O5ilPeUq3nw8AHB+W1OlIAIDjhRAGANCBEAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGANCBEAYA0MFYQ1hVXVBVX6qqW6vqtUPWX1ZVd1bVZ+e/XjbOegAAJsWJ43riqjohyVVJ/nWS25NcW1W7WmtfOGTT97bWLh9XHQAAk2icR8KenuTW1tptrbV7k+xMcvEYfx4AwHGjWmvjeeKqS5Jc0Fp72fzyS5M8Y+FRr6q6LMl/TnJnki8n+Y+ttb8Z8lzbkmxLkvXr15+3c+fOsdR8LOzfvz9r1qzpXcbE0Zfh9GU0vRlOX0bTm+H0ZbTF6M3WrVuvb61tHrZubKcjj9CfJnlPa+17VfULSd6V5DmHbtRam04ynSSbN29uW7ZsWdQij8bs7Gwmub5e9GU4fRlNb4bTl9H0Zjh9Ga13b8Z5OnJfkrMWLJ85P/Z9rbW7W2vfm1/83STnjbEeAICJMc4Qdm2Ss6vqcVW1KsmlSXYt3KCqHrtg8aIkt4yxHgCAiTG205Gttfur6vIkH0lyQpJ3tNZurqo3JbmutbYryS9X1UVJ7k/yzSSXjaseAIBJMtZrwlprVye5+pCxNyz4/leT/Oo4awAAmERmzAcA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwDoQAgDAOhACAMA6EAIAwCWl5mZZNOm5PrrB48zM13KOLHLTwUA6GFmJtm2LZmbGyzv2TNYTpKpqUUtxZEwAGD52L79BwHsgLm5wfgiE8IAgOVj796jGx8jIQwAWD42bDi68TESwgCA5WPHjmT16oPHVq8ejC8yIQwAWD6mppLp6WTjxsHyxo2D5UW+KD9xdyQAsNxMTQ2+ZmeT3bu7leFIGABAB0IYAEAHQhgAQAdCGABAB0IYAEAHYw1hVXVBVX2pqm6tqtc+xHYvqKpWVZvHWQ8AwKQYWwirqhOSXJXkeUnOSfLiqjpnyHYnJ3llkk+NqxYAgEkzziNhT09ya2vtttbavUl2Jrl4yHZvTvJfk3x3jLUAAEyUaq2N54mrLklyQWvtZfPLL03yjNba5Qu2eWqS7a21F1TVbJLXtNauG/Jc25JsS5L169eft3PnzrHUfCzs378/a9as6V3GxNGX4fRlNL0ZTl9G05vh9GW0xejN1q1br2+tDb3cqtuM+VW1IsmVSS473Lattekk00myefPmtmXLlrHW9k8xOzubSa6vF30ZTl9G05vh9GU0vRlOX0br3Ztxno7cl+SsBctnzo8dcHKSc5PMVtXuJD+SZJeL8wGA5WCcIezaJGdX1eOqalWSS5PsOrCytfb3rbW1rbVNrbVNST6Z5KJhpyMBAJaasYWw1tr9SS5P8pEktyR5X2vt5qp6U1VdNK6fCwBwPBjrNWGttauTXH3I2BtGbLtlnLUAAEwSM+YDAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdHDYEFZVv1RVpy5GMQAAy8WRHAlbn+TaqnpfVV1QVTXuogAAlrrDhrDW2q8lOTvJ7yW5LMlXquotVfWEMdcGALBkHdE1Ya21luTr81/3Jzk1yfur6q1jrA0AYMk68XAbVNUrk/xskruS/G6SX2mt3VdVK5J8Jcl/Gm+JAABLz2FDWJLTkvy71tqehYOttQer6vnjKQsAYGk7bAhrrb3xIdbdcmzLAQBYHswTBgDQgRAGANDBWEPY/LxiX6qqW6vqtUPWv7yqPl9Vn62q/1NV54yzHgCASTG2EFZVJyS5KsnzkpyT5MVDQta7W2v/orX25CRvTXLluOoBAJgk4zwS9vQkt7bWbmut3ZtkZ5KLF27QWvuHBYuPTNLGWA8AwMSowTysY3jiqkuSXNBae9n88kuTPKO1dvkh270iyauTrErynNbaV4Y817Yk25Jk/fr15+3cuXMsNR8L+/fvz5o1a3qXMXH0ZTh9GU1vhtOX0fRmOH0ZbTF6s3Xr1utba5uHrTuSecLGqrV2VZKrquolSX4tyc8N2WY6yXSSbN68uW3ZsmVRazwas7OzmeT6etGX4fRlNL0ZTl9G05vh9GW03r0Z5+nIfUnOWrB85vzYKDuT/PQY6wEAmBjjDGHXJjm7qh5XVauSXJpk18INqursBYv/JoOPQQIAWPLGdjqytXZ/VV2e5CNJTkjyjtbazVX1piTXtdZ2Jbm8qs5Pcl+Sv8uQU5EAAEvRWK8Ja61dneTqQ8besOD7V47z5wMATCoz5gMAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IITBKDMzyaZNyYoVg8eZmd4VAbCEnNi7AJhIMzPJtm3J3Nxgec+ewXKSTE31qwuAJcORMBhm+/YfBLAD5uYG4wBwDAhhMMzevUc3DgBHSQiDYTZsOLpxADhKQhgMs2NHsnr1wWOrVw/GAeAYEMJgmKmpZHo62bgxqRo8Tk+7KB+AY8bdkTDK1JTQBcDYOBIGANCBEAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGo8zMJJs2JStWDB5nZnpXBMAScmLvAmAizcwk27Ylc3OD5T17BstJMjXVry4AlgxHwmCY7dt/EMAOmJsbjAPAMSCEwTB79x7dOAAcJSEMhtmw4ejGAeAoCWEwzI4dyerVB4+tXj0YB4BjQAiDYaamkunpZOPGpGrwOD3tonwAjhl3R8IoU1NCFwBj40gYAEAHQhgAQAdCGABAB0IYAEAHQhgAQAdCGABAB0IYAEAHQhgAQAdjDWFVdUFVfamqbq2q1w5Z/+qq+kJVfa6q/qKqNo6zHgCASTG2EFZVJyS5KsnzkpyT5MVVdc4hm30myebW2g8neX+St46rHgCASTLOI2FPT3Jra+221tq9SXYmuXjhBq21v2ytzc0vfjLJmWOsBwBgYlRrbTxPXHVJkgtaay+bX35pkme01i4fsf3bkny9tfYbQ9ZtS7ItSdavX3/ezp07x1LzsbB///6sWbOmdxkTR1+G05fR9GY4fRlNb4bTl9EWozdbt269vrW2edi6ifgA76r6mSSbk/z4sPWttekk00myefPmtmXLlsUr7ijNzs5mkuvrRV+G05fR9GY4fRlNb4bTl9F692acIWxfkrMWLJ85P3aQqjo/yfYkP95a+94Y6wEAmBjjvCbs2iRnV9XjqmpVkkuT7Fq4QVU9Jcnbk1zUWvvGGGsBAJgoYwthrbX7k1ye5CNJbknyvtbazVX1pqq6aH6z30yyJskfVdVnq2rXiKcDAFhSxnpNWGvt6iRXHzL2hgXfnz/Onw8AMKnMmA8A0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGANCBEAYA0IEQBgDQgRAGANCBELbczMwkmzYlK1YMHmdmelcEAMvSib0LYBHNzCTbtiVzc4PlPXsGy0kyNdWvLgBYhhwJW062b/9BADtgbm4wDgAsKiFsOdm79+jGAYCxEcKWkw0bjm4cABgbIWw52bEjWb364LHVqwfjAMCiEsKWk6mpZHo62bgxqRo8Tk+7KB8AOnB35HIzNSV0AcAEcCQMAKADIQwAoAMhDACgAyEMAKADIWy58dmRADAR3B25nPjsSACYGI6ELSc+OxIAJoYQtpz47EgAmBhC2HLisyMBYGIIYcuJz44EgIkhhC0nPjsSACaGuyOXG58dCQATwZEwAIAOhDAAgA6EMACADoQwAIAOhDAAgA6EMACADoQwAIAOhDAAgA6EsMU2M5Ns2pSsWDF4nJnpXREA0IEZ8xfTzEyybVsyNzdY3rNnsJyYxR4AlhlHwhbT9u0/CGAHzM0NxgGAZUUIW0x79x7dOACwZAlhi+m0045uHABYsoQwAIAOhLDF9M1vHt04ALBkCWGLacOGoxsHAJYsIWwx7diRrF598Njq1YNxAGBZEcIW09RUMj2dbNyYVA0ep6fNEQYAy5DJWhfb1JTQBQA4EgYA0IMQBgDQwVhDWFVdUFVfqqpbq+q1Q9Y/u6puqKr7q+qScdYCADBJxhbCquqEJFcleV6Sc5K8uKrOOWSzvUkuS/LucdUBADCJxnlh/tOT3Npauy1JqmpnkouTfOHABq213fPrHhxjHQAAE6daa+N54sHpxQtaay+bX35pkme01i4fsu07k3yotfb+Ec+1Lcm2JFm/fv15O3fuHEvNx8L+/fuzZs2a3mVMHH0ZTl9G05vh9GU0vRlOX0ZbjN5s3br1+tba5mHrjospKlpr00mmk2Tz5s1ty5Yti1/EzEyyfXuyd+9ghvsdO4ZONTE7O5su9U04fRlOX0bTm+H0ZTS9GU5fRuvdm3GGsH1JzlqwfOb82PFnZibZti2Zmxss79kzWE7M+QUAPCzjvDvy2iRnV9XjqmpVkkuT7Brjzxuf7dt/EMAOmJsbjAMAPAxjC2GttfuTXJ7kI0luSfK+1trNVfWmqrooSarqaVV1e5IXJnl7Vd08rnr+SfbuPbpxAIDDGOs1Ya21q5NcfcjYGxZ8f20Gpykn22mnJXffPXwcAOBhMGM+AEAHQtiR+OY3j24cAOAwhLCFZmaSTZuSFSsGjzMzg/ENG4ZvP2ocAOAwjot5whbFQ01DsWPHweuSZPXqwTgAwMMghB3wUNNQ7N79g20OM1krAMCREMIOONw0FFNTQhcAcMy4JuyAUdNNmIYCABgDIQwAoAMh7ADTUAAAi0gIO8A0FADAIhLCDtixYzDtxEKmoQAAxkQIO2BqKpmeTjZuTKoGj9PT7ogEAMbCFBULmYYCAFgkjoQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHQghAEAdCCEAQB0IIQBAHRQrbXeNRyVqrozyZ7edTyEtUnu6l3EBNKX4fRlNL0ZTl9G05vh9GW0xejNxtba6cNWHHchbNJV1XWttc2965g0+jKcvoymN8Ppy2h6M5y+jNa7N05HAgB0IIQBAHQghB17070LmFD6Mpy+jKY3w+nLaHoznL6M1rU3rgkDAOjAkTAAgA6EsCNUVRdU1Zeq6taqeu2Q9a+uqi9U1eeq6i+qauOCdQ9U1Wfnv3YtbuXjdwS9uayq7lzQg5ctWPdzVfWV+a+fW9zKx+sI+vLfFvTky1X1rQXrluw+U1XvqKpvVNVNI9ZXVf33+b59rqqeumDdUt5fDteXqfl+fL6q/rqq/uWCdbvnxz9bVdctXtWL4wh6s6Wq/n7B78wbFqx7yN/D49kR9OVXFvTkpvn3ldPm1y3Zfaaqzqqqv5z/m3xzVb1yyDaT8T7TWvN1mK8kJyT5apLHJ1mV5MYk5xyyzdYkq+e//8Uk712wbn/v19C5N5cleduQf3taktvmH0+d//7U3q9psfpyyPa/lOQdy2SfeXaSpya5acT6n0ryZ0kqyY8k+dRS31+OsC8/euD1Jnnegb7ML+9Osrb3a+jYmy1JPjRk/Kh+D4+3r8P15ZBtL0zy8eWwzyR5bJKnzn9/cpIvD/m7NBHvM46EHZmnJ7m1tXZba+3eJDuTXLxwg9baX7bW5uYXP5nkzEWusZfD9uYhPDfJn7fWvtla+7skf57kgjHVudiOti8vTvKeRamss9baJ5J88yE2uTjJ/2oDn0zy6Kp6bJb2/nLYvrTW/nr+dSfL6z3mSPaZUf4p708T7yj7spzeY+5ord0w//23k9yS5IxDNpuI9xkh7MickeRvFizfnn/8H7rQf8ggYR/wiKq6rqo+WVU/PY4COzrS3rxg/pDv+6vqrKP8t8ejI35t86euH5fk4wuGl/I+czijereU95ejdeh7TEvy0aq6vqq2daqpt39VVTdW1Z9V1T+fH7PPJKmq1RkEiQ8sGF4W+0xVbUrylCSfOmTVRLzPnDiuJ16uqupnkmxO8uMLhje21vZV1eOTfLyqPt9a+2qfCrv40yTvaa19r6p+Icm7kjync02T5NIk72+tPbBgbLnvM4xQVVszCGHPXDD8zPn9ZV2SP6+qL84fJVkubsjgd2Z/Vf1Ukg8mObtzTZPkwiT/t7W28KjZkt9nqmpNBsHzVa21f+hdzzCOhB2ZfUnOWrB85vzYQarq/CTbk1zUWvvegfHW2r75x9uSzGaQypeKw/amtXb3gn78bpLzjvTfHseO5rVdmkNOEyzxfeZwRvVuKe8vR6SqfjiD36GLW2t3HxhfsL98I8mfZHAabtlorf1Da23//PdXJ1lZVWtjnzngod5jluQ+U1UrMwhgM621Px6yyUS8zwhhR+baJGdX1eOqalUGO/RBd6xV1VOSvD2DAPaNBeOnVtVJ89+vTfJjSb6waJWP35H05rELFi/K4Px8knwkyU/O9+jUJD85P7YUHLYvSVJVT8zg4s//t2Bsqe8zh7Mryc/O3730I0n+vrV2R5b2/nJYVbUhyR8neWlr7csLxh9ZVScf+D6Dvgy9W26pqqrHVFXNf//0DP623Z0j/D1cyqrqURmcmfnfC8aW9D4zvy/8XpJbWmtXjthsIt5nnI48Aq21+6vq8gz+I07I4C62m6vqTUmua63tSvKbSdYk+aP594K9rbWLkjwpydur6sEM3hj+S2ttyfxBPcLe/HJVXZTk/gwuIr1s/t9+s6renMEbZZK86ZDD5cetI+xLMvijsLPN35Yzb0nvM1X1ngzuZltbVbcneWOSlUnSWvufSa7O4M6lW5PMJfn38+uW7P6SHFFf3pDkh5L8j/n3mPvb4IOH1yf5k/mxE5O8u7V2zaK/gDE6gt5ckuQXq+r+JN9Jcun879TQ38MOL2EsjqAvSfJvk3y0tXbPgn+61PeZH0vy0iSfr6rPzo+9LsmGZLLeZ8yYDwDQgdORAAAdCGEAAB0IYQAAHQhhAAAdCGEAAB0IYQAAHQhhAAAdCGHAslVVT5v/YPlHzM8ifnNVndu7LmB5MFkrsKxV1W8keUSSf5bk9tbaf+5cErBMCGHAsjb/mYLXJvlukh9trT3QuSRgmXA6EljufiiDz309OYMjYgCLwpEwYFmrql1JdiZ5XJLHttYu71wSsEyc2LsAgF6q6meT3Ndae3dVnZDkr6vqOa21j/euDVj6HAkDAOjANWEAAB0IYQAAHQhhAAAdCGEAAB0IYQAAHQhhAAAdCGEAAB0IYQAAHfx/KLjtCcnrRGIAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# initial guesses of parameters\n",
        "P0 = [1, 1]"
      ],
      "metadata": {
        "id": "RXaBDa9bsYcg"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# the function that we want to fit\n",
        "def f(X, P0):\n",
        "  return (P0[0]*X)/(P0[1]+X)"
      ],
      "metadata": {
        "id": "_7oDKqSOymK9"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The **Jacobian matrix**, in case you don't know it, is defined as follows\n",
        "\n",
        "$\\mathbf{J}=\\left[\\begin{array}{ccc}\n",
        "\\frac{\\partial f_1}{\\partial p_1} & \\frac{\\partial f_1}{\\partial p_2} & \\frac{\\partial f_1}{\\partial p_3} \\\\\n",
        "    \\frac{\\partial f_2}{\\partial p_1} & \\frac{\\partial f_2}{\\partial p_2} & \\frac{\\partial f_2}{\\partial p_3} \\\\\n",
        "    \\frac{\\partial f_3}{\\partial p_1} & \\frac{\\partial f_3}{\\partial p_2} & \\frac{\\partial f_3}{\\partial p_3}\n",
        "\\end{array}\\right] \\in \\mathbb{R^{i\\times j}}$\n",
        "\n",
        "The 1st and 2nd columns of our Jacobian matrix are derivatives of $f(x)$ wrt to $p_1$ and $p_2$. For $i > j$, $i$ is number of sample points (pairs of $x$ and $y$), and $j$ is number of parameters.\n",
        "In this case,\n",
        "\n",
        "$\\mathbf{J} = \\left[ -\\frac{x}{p_2+x} \\space\\space \\frac{p_1x}{(p_2+x)^2}  \\right]$\n",
        "\n",
        "For Gauss-Newton, the algorithm is ($k$ is iteration number):\n",
        "\n",
        "1.   compute search direction $d_k = (J^TJ)^{-1}J^Tr  \\in \\mathbb{R}^j$\n",
        "2.   optionally compute step length $t_k \\in \\mathbb{R}$\n",
        "3.   iterate $p_{k+1} = p_k + t_kd_k$ \n",
        "4.   check stopping criterien\n",
        "5.   if criterien is not met, redo from step 1.\n",
        "\n",
        "$p_{k+1}$ are the fitted parameters."
      ],
      "metadata": {
        "id": "P66K5Me4ZhWX"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Jacobian\n",
        "J = np.zeros([len(X), len(P0)])\n",
        "Iter = 0"
      ],
      "metadata": {
        "id": "IqFghB7rxmEr"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "while True:\n",
        "  Iter += 1\n",
        "  j1 = -1*(X/(P0[1]+X))\n",
        "  j2 = (P0[0]*X)/((P0[1]+X)**2)\n",
        "  J[:, 0] = j1\n",
        "  J[:, 1] = j2\n",
        "\n",
        "  r = Y - f(X, P0) # objective function\n",
        "  t1 = np.linalg.inv(np.dot(J.T, J))\n",
        "  t2 = np.dot(t1, J.T)\n",
        "  t3 = np.dot(t2, r)\n",
        "  P1 = P0 - t3 # Gauss-Newton update\n",
        "  t4 = abs(P1-P0)\n",
        "\n",
        "  if max(t4) <= 1e-6:\n",
        "    break\n",
        "  P0 = P1\n",
        "\n",
        "x1 = float('{:.4f}'.format(P0[0]))\n",
        "x2 = float('{:.4f}'.format(P0[1]))\n",
        "\n",
        "print(f'The parameters are:\\nx1 = {x1}\\nx2 = {x2}\\nIterations = {Iter}')\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "2TdRgsMfxxwu",
        "outputId": "28371287-109c-4c48-ff8d-9c58078f34f3"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The parameters are:\n",
            "x1 = 1.0313\n",
            "x2 = 2.2189\n",
            "Iterations = 6\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Plotting"
      ],
      "metadata": {
        "id": "IasSj6Zw0RUR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pred = f(X, P0)\n",
        "plt.figure(figsize=(10,6))\n",
        "plt.scatter(X, Y, c='red', marker='o', label=\"Data\")\n",
        "plt.plot(X, pred, '--m', label=\"Model\")\n",
        "plt.legend()\n",
        "plt.xlabel(\"x\")\n",
        "plt.ylabel(\"y\")\n",
        "plt.legend(fontsize=14)\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 388
        },
        "id": "9DZ1Fz5l0I3A",
        "outputId": "c77db3ed-8264-4de4-c8ef-84bad43ce87d"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 720x432 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFzCAYAAAB2A95GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXzU1b3/8dfJZN9DNpJAEpawiawBtHUBq9UWl1rtdUm94kbd2tpeaxeqVVt6r60/e9tbQaPFpcZSq124ytXWtqlttbKpIKsxApJJQgIEyL7M+f0xAQNOIIHMfGd5Px+PPDJzvt/MfObwzeTN+Z75HmOtRUREREQCK8rpAkREREQikUKYiIiIiAMUwkREREQcoBAmIiIi4gCFMBEREREHKISJiIiIOCDa6QIGKysryxYXFztdRr9aWlpISkpyuoygo37xTf3SP/WNb+qX/qlvfFO/9C8QfbN27dpGa222r20hF8KKi4tZs2aN02X0q7Kykrlz5zpdRtBRv/imfumf+sY39Uv/1De+qV/6F4i+Mcbs6G+bTkeKiIiIOEAhTERERMQBCmEiIiIiDlAIExEREXGAQpiIiIiIA0Lu05HH4vF4aGxspKmpiZ6eHkdqSEtLY/PmzY48dyC5XC7S09PJysoiKkpZXkREZLDCKoTt2rULYwzFxcXExMRgjAl4DQcPHiQlJSXgzxtI1lq6urqor69n165dFBYWOl2SiIhIyAmrIYyWlhYKCgqIjY11JIBFCmMMsbGxFBQU0NLS4nQ5IiIiISmsQhigU2MBpL4WERE5cforKiIiIpGlogKKi2HtWu/3igpHygirOWEiIiIix1RRAQsXQmur9/6OHd77AGVlAS1FI2ERoLGxEWMMlZWVA/6Ze++9l8mTJ/uvKBEREScsWvRRADuktdXbHmAKYUFgwYIFGGO44YYbPrbtm9/8JsYYLrzwQgcqExERCTM7dw6u3Y8UwoLEyJEjee655474tGF3dzdPP/20LgEhIiIyVPr7m+rA31qFsCAxZcoUSkpKeO655w63vfTSS8THxzN37tzDbR6Ph+9///uMHDmSuLg4Tj31VP7whz8c8VirV69m5syZxMfHM336dN58882PPd+mTZuYP38+KSkp5OTkcNVVV1FXV+e31yciIhIUFi+GxMQj2xITve0BFhET89+a+9bH2nL+LYeCWwvoae1h/WfXf2z78AXDyVuQR2djJxsv3/ix7QW3FJBzRQ7tH7az+ZqPrpDf09ND6d9LT6jOG264gWXLlnHdddcBHL5dXV19eJ+f/vSn/PjHP+aRRx6htLSUZ555hs9//vOsXbuWadOm0dzczPz58zn77LN56qmnqKmp4Y477jjieWpraznrrLO44YYbePDBB+nq6mLRokVccsklvPHGG7r0hIiIhK9Dk+8PzQErKvIGsABPygeNhAWVq6++mjVr1vDee+9RV1fHyy+/zIIFC47Y58EHH+TOO+/k6quvZty4cdx///2ceeaZPPjggwA8++yzdHZ28sQTTzB58mTOP/98Fh012XDp0qVMnTqVBx54gIkTJzJlyhSefvppVq1axZo1awL1ckVERJxRVgbbt8PMmd7vDgQwiJCRsOmV0/vd5kp0HXN7bFbsMbfHj4w/YvvBgwdPrEggIyODSy+9lGXLlpGens7cuXOPmA924MAB3G43n/zkJ4/4uTPOOIOVK1cCsHnzZqZMmUJycvLh7aeffvoR+69du5bXXnvtiH0Oef/995k9e/YJvwYREREZmIgIYaHk+uuv59prryU5OZn7779/wD83mGWaPB4P8+fPPzx61ldubu6AH0dEREROnE5HBplPfepTxMbG0tjYyOc+97kjtqWmppKfn88///nPI9r/8Y9/MGnSJAAmTpzIhg0bjviU5b/+9a8j9p8xYwYbN26kqKiIsWPHHvEV7ouPi4iIBAuFsCBjjGH9+vV88MEHxMXFfWz7N77xDR588EF+9atfsW3bNu655x7+/ve/c+eddwLeeWXR0dFcf/31bNy4kT/96U8sPuoTH7fddhv79+/niiuu4M0336S6uppXX32VhQsXntTpVBERERk4nY4MQscajfrKV77CwYMHueuuu6ivr2f8+PG88MILTJ06FYDk5GRefPFFbrnlFmbMmMGECRN44IEHuPjiiw8/xqHRtG9/+9tccMEFtLe3U1hYyKc//WmfwU9ERESGnkJYEHjyyScHvD0qKoq7776bu+++u9/958yZw7p1645os9Yecb+kpITnn3++38e49957uffee49Zl4iIiJw4nY4UERERcYBCmIiIiIgDFMJEREREHKAQJiIiIuIAhTARERERByiEiYiIiDhAIUxERETEAQphIiIiIg5QCBMRERFxgEKYiIiIiAMUwoLAggULMMZgjCEmJoacnBzmzZvHww8/TFdX14Afp7KyEmMMjY2NfqxWREREhoJCWJA499xzqa2tZfv27fzxj3/koosu4nvf+x5nnnkmLS0tTpcnIiIiQ0whLEjExcUxfPhwCgoKmDZtGl//+teprKxk3bp1/OhHPwLgmWeeYdasWaSkpJCTk8MXvvAFampqANi+fTvz5s0DIDs7G2MMCxYsAODll1/mzDPPJCMjg2HDhnH++eezefNmR16niIiIeCmE+VJRAcXFEBXl/V5R4UgZkydP5oILLuCFF14AoLOzk/vuu4933nmHF198kcbGRq666ioARo4ceXi/jRs3Ultby09/+lMAWlpauOOOO1i1ahWVlZWkpaVx0UUX0dnZ6cjrEhEREYh2uoCgU1EBCxdCa6v3/o4d3vsAZWUBL2fSpEm8+uqrAFx//fWH20ePHs3SpUuZOHEiu3btYsSIEQwbNgyAnJwcsrKyDu972WWXHfGYTzzxBKmpqaxatYozzjgjAK9CREREjqaRsKMtWvRRADuktdXb7gBrLcYYANatW8cll1xCUVERKSkplJaWArBz585jPsb777/P1VdfzZgxY0hNTSU3NxePx3PcnxMRERH/UQg7Wn/BxKHAsmnTJkaPHk1LSwvnn38+iYmJ/PKXv2T16tW8/PLLAMc9rXjhhRfS0NDAo48+yptvvslbb71FdHS0TkeKiIg4SCHsaIWFg2v3o3fffZeXX36Zyy+/nC1bttDY2MgPf/hDzjrrLCZMmMDu3buP2D82NhaAnp6ew2179uxhy5YtfOc73+Hcc89l4sSJHDx4kO7u7oC+FhERETmSQtjRFi+GxMQj2xITve1+1NHRQV1dHW63m3feeYeHHnqIuXPnMnPmTO68804KCwuJi4vj5z//OdXV1bz00kvcfffdRzxGUVERxhheeuklGhoaaG5uJiMjg6ysLB577DGqqqr429/+xs0330x0tKYDioiIOEkh7GhlZVBeDkVFYIz3e3m53yflv/rqq+Tl5VFYWMinPvUpVqxYwb333strr71GUlIS2dnZPPXUU/z+979n0qRJ3HfffTz00ENHPEZBQQH33XcfixYtIjc3l9tvv52oqCh+/etfs379eiZPnsxtt93G97//feLi4vz6ekREROTYNBziS1lZQD8J+eSTT/Lkk08ed78rrriCK6644og2a+0R9+++++6PjZCdc845vPvuu0e0NTc3n1ixIiIiMiT8OhJmjLnAGLPVGFNljPmWj+0LjDENxpi3e79u9Gc9IiIiIsHCbyNhxhgX8DBwHrALWG2MWWGt3XTUrr+21t7urzpEREREgpE/R8JmA1XW2mprbSewHLjEj88nIiIiEjLM0XOKhuyBjbkcuMBae2Pv/WuAOX1HvYwxC4D/BBqAbcDXrLUf+nishcBCgNzc3JnLly/3+ZxpaWmMHTt2iF/J4PT09OByuRytIZCqqqrYv3//cfdrbm4mOTk5ABWFFvVL/9Q3vqlf+qe+8U390r9A9M28efPWWmtLfW1zemL+/wK/stZ2GGO+BDwFnHP0TtbacqAcoLS01M6dO9fng23evJnk5OTDV5h3wsGDB0lJSXHs+QPJWkt8fDzTp08/7r6VlZX09+8WydQv/VPf+KZ+6Z/6xjf1S/+c7ht/no6sAUb2uT+it+0wa+0ea21H793HgZkn84QxMTG0tbWdzEPIILS1tRETE+N0GSIiIiHJnyFsNVBijBlljIkFrgRW9N3BGJPX5+7FwOaTecKcnBxqampobW392KUbZOhYa2ltbaWmpoacnBynyxEREQlJfjsdaa3tNsbcDrwCuIBl1tqNxpj7gTXW2hXAV4wxFwPdwF5gwck8Z2pqKgBut5uurq6TeagT1t7eTnx8vCPPHUgxMTHk5uYe7nMREREZHL/OCbPWrgRWHtV2T5/b3wa+PZTPmZqa6mgwqKysHNAcKREREYlsTk/MFxEREQmYnvYeGp5rwNPugXHO1qIQJiIiImGvrboN9yNuapfV0r2nm7Qz0uD7ztakBbxFREQkLFnPRx/Sq/52NR8+9CHpc9OZ+uepTHttmoOVeWkkTERERMJKZ0Mntb+opfbRWk5deSpJE5MY/cPRjH1oLHEFcU6Xd5hCmIiIiIQ8ay0H3jhAzZIaGn7TgO20pJ+TjqfNA0DCmASHK/w4hTAREREJWdZajDH0HOjhnXPfwcQY8m/OJ//mfJImJjld3jEphImIiEjIadncgnupm5ZNLUx7dRrRadFM+b8pJM9MJjo5NOJNaFQpIiIiEc/T5aHxD424l7hp+msTJtaQ/YVsetp6cCW4SD873ekSB0UhTERERELC7md3s2XBFuKK4hj1n6PIuz6P2JxYp8s6YQphIiIiEnSstTT9tYmaJTVknJNBwa0FZH8hm+hh0WR+NhPjMk6XeNIUwkRERCRodDV1Uf90Pe6lblq3tBI9LJr0M72nGV2JLrIuynK4wqGjECYiIiJBY/NVm9n78l5S5qQw4akJZH8hG1eCy+my/EIhTERERBzR095Dw/MN1JbXMum5ScQNj6P4vmJG/WAUKTNTnC7P7xTCREREJKDaPmjD/aibul/U0dXYRUJJAh07OogbHkfq7FSnywsYhTAREREJmM6GTlaNW4X1WLIuySL/1nwyzsnARIX+RPvBUggTERERv+ls7KRuWR3tH7Qzbuk4YrNjGb9sPOlz04kfGe90eY5SCBMREZEhZa3lwJsHcC9xs/u53dgOS/q8dDydHqJioxh+zXCnSwwKCmEiIiIypGofq2Xbl7bhSnGRd2MeBbcUkHRKcK/j6ASFMBERETkpLVu86zimz0sn+3PZZH0uC+ux5JblEp2iqNEf9YyIiIgMmqfbw54/7KFmSQ1Nf2nCxBjvEkKfg9icWApuLnC6xKCnECYiIiKDtv6C9TT9uYm4wjhG/bB3Hcfc0F3H0QkKYSIiInJM1lqa/tZE3bI6xi0dhyvJxYivjmDEV0aQOT881nF0gkKYiIiI+NS9v5u6p+u86zhubiU6I5rmW5tJOy0trNZwdIpCmIiIiHxM+652Vk1YhafFQ8rsFMY/MZ6cK3LCdh1HJyiEiYiICJ4ODw0vNNCxq4PCuwqJHxFP4V2FDPvsMFJLI2cpoUBSCBMREYlg7TvacT/qpvbxWroaukiaksTI/xiJcRmK7yl2urywphAmIiISodyPudn2pW1gIPOiTApuKyDjU5G5jqMTFMJEREQiRNeeLmqfqCXtE2mkfSKN9LPTKfxOIfkL84kvjOx1HJ2gECYiIhLOLBxYdYCaJTXsXu5dx7Hou0WkfSKNxHGJjP7BaKcrjFgKYSIiIuHsO7DuX+twJbvIuyGP/FvySZ6c7HRVgkKYiIhIWGnd1kr9M/UU3VNEVHQUnA4l/15C7he1jmOw0b+GiIhIiPN0e9jz4h7cS9zs+9M+TLQh86JMUmelwsVQMFfrOAYjhTAREZEQ1vZBG2+f9TYduzqIGxHHqB+MYvgNw4kbHud0aXIcCmEiIiIhxFrL/r/vp6Omg9yrcokviif9nHSyLs0i88JM7ylICQkKYSIiIiGg+0A39c/UU7OkhtaNrSSUJJBzZQ4myjDxqYlOlycnQCFMREQkyNU+WUvVl6voae4heWYy45d513E0RhdVDWUKYSIiIkHG0+mh4bcNpMxIIXFcIonjEsm6LIuCWwtIna11HMOFQpiIiEiQaN/ZjrvcTe1jtXTt7qLwO4WMXjz68BXuJbwohImIiASBTWWb2L18NwCZ8zPJvzWfYZ8e5nBV4k8KYSIiIg7o2ttF4+8bGX7dcIwxxBfHU/it3nUci7SOYyTQ51hFREQC6MCaA2y5bgtvFLzB1hu20rKhBYDRi0czevFoBbBAqKiA4mJYu9b7vaLCkTI0EiYiIhIAbe+3semqTRxcfZCopCiGLxjuXcdxitZxDKiKCli4EFpbvfd37PDeBygrC2gpCmEiIiJ+0lrVSmdNJ+lnpxObH4uJNYz9n7EMv2Y40Wn6E+yIRYs+CmCHtLZ62xXCREREQpen28Pel/ZSs6SGfX/cR+LERGZtnIUrwcWMf8xwujzZuXNw7X6kECYiIjJE6pfXU31XNR0fdhBbEEvx/cXk3Zini6oGk8JC7ylIX+0Bpon5IiIiJ8haS9Pfm+io6wAgKi6KxPGJnPLbUzht+2kU311MXJ4W0g4qixdDYuKRbYmJ3vYA00iYiIjIIHUf7Ka+oh73EjctG1oovr+Y4ruLyb40m+xLs50uT47l0LyvRYu834uKvAEswPPBQCFMRERkwKy1VN1RRd0TdfQc7CF5ejLjHx9PzpU5Tpcmg1FW5v2qrITt2x0rQyFMRETkGDydHvb/cz8Z8zIwxtC1u4usS73rOKbMTtF8LzlhCmEiIiI+tH/YTm15Le7H3HTVdzH7vdkkjk1k4rMTFbxkSCiEiYiI9NFW3cb7//E+jSsawX60jmPC6AQABTAZMgphIiIS8br2ddFZ10nSxCRcqS4OrjlI4V2F5H0pj4TiBKfLkzClECYiIhHr4NqD1CytYfezu0melsyM12cQmxXLaTtOw0RpxEv8y6/XCTPGXGCM2WqMqTLGfOsY+11mjLHGmFJ/1iMiIgLQ+GIja+esZW3pWnb/aje51+RSsqTk8HYFMAkEv42EGWNcwMPAecAuYLUxZoW1dtNR+6UAXwXe9FctIiIibe+3EZMbQ3RyNB0fdtBzoIexPxvL8H/XOo7iDH+OhM0Gqqy11dbaTmA5cImP/b4PPAC0+7EWERGJQLbH0vi/jaz/zHreHPsm9c/UA5B3Ux6zNs1ixJdHKICJY4y11j8PbMzlwAXW2ht7718DzLHW3t5nnxnAImvtZcaYSuBOa+0aH4+1EFgIkJubO3P58uV+qXkoNDc3k5yc7HQZQUf94pv6pX/qG9/UL/07om88eP/rvwKoB7KAC4H5vbcjiI6Z/gWib+bNm7fWWutzupVj8d8YEwU8BCw43r7W2nKgHKC0tNTOnTvXr7WdjMrKSoK5PqeoX3xTv/RPfeOb+qV/lX+tZPYI77W8ANbdvY6oU6IoeLiAzIsziYqJzOWSdcz0z+m+8WcIqwFG9rk/orftkBRgMlDZe82V4cAKY8zFvkbDREREfOlu7mZ3xW74EazetZrTa04nNiuWqa9OxZXgcro8kX75M4StBkqMMaPwhq8rgasPbbTW7qfPoPCxTkeKiIgcrX1XOx8+8CF1T3nXcWQslPy8BFeiN3gpgEmw89vYrLW2G7gdeAXYDDxnrd1ojLnfGHOxv55XRETCl6fLQ2d9p/d2u4fax2vJuiSL6a9Ph3LIvyn/cAgTCXZ+nRNmrV0JrDyq7Z5+9p3rz1pERCR0ddR04C53U/tYLalzUpn8u8kkjk3kE/WfIDq1909ZpaMligyaPpcrIiJBa/8/9/PhTz6k8feN4IFhnxlG3k15h7cfDmAiIUhHr4iIBJWupi5cyS6ioqPY9+o+miqbGPkfI8n/0keLaIuEA4UwEREJCgffOoh7iZv6Z+uZ+MuJZH8+mxFfH8HIb47EFa95XhJ+FMJERMQxtsdS/2w97iVuDvzrAFEJUeSW5ZI40Xutr+gU/ZmS8KWjW0REAq77QLd3PpeBHffvgCgY+99jyb02l5j0GKfLEwkIhTAREQkI22PZ+/JeapbUcHDVQU7beRquBBfTKqcRmx9L74W7RSKGQpiIiPhV154uah+vxf2Im/bt7cQOjyX/1nxsl4UEiCuIc7pEEUcohImIyJCz1uJp9+BKcNH6XivV36omfW46o380mqzPZUXsOo4ifSmEiYjIkOlp6Tk80T5lTgrjHxlP6pxUZm+bTWJJotPliQQVhTARETlpLVtacC91U/dkHT0HekiakkTa6WkAGGMUwER8UAgTEZET4unyYKINxhhq/qeG2sdqyf5CNgW3FpD6iVRNtBc5DoUwEREZlA53B7WP1eIudzPp15NIPyOdou8WUfy9YmJzYp0uTyRkKISJiMhxWWtpqmzCvcRNw+8aoAeGXTCMqHjvBPu4PH3CUWSwFMJERKRf1mMxUQbbadn0b5uwHsvIr40k/+Z8EsZoHUeRk6EQJiIiH3Pw7YO4l7o58PoBSt8uJSouiimvTCFxYiKuBK3jKDIUFMJERASAnvYeGp5vOBy+ouKjyLk6h57mHqLTokmZkeJ0iSJhRSFMRCTCWWsxxrDvT/vYcs0WEkoSGPPQGIZfO5yYYVrHUcRfFMJERCKQ9Vj2vrIX9xI3yTOTGXXvKDI/m8nUV6eSPi8dE6XLS4j4m0KYiEgE6WzspO6JOu86jtXtxOTGkH5OOgDGZcj4VIbDFYpEDoUwEZEI8t7t79Hw6wbSzkpj9A9Hk3VpFlGxWsdRxAkKYSIiYaqntYfdv9pNzdIaJj4zkaQJSRTfU0zRd4tInpzsdHkiEU8hTEQkzLRuaz28jmN3UzdJk5PoauwCIGlSksPVicghCmEiImGku7mbNdPWYLst2Zdlk39rPmlnpGkdR5EgpBAmIhLCOmo7qH28lua3mpn828lEJ0cz6deTSJ2dSmyu1nEUCWYKYSIiIcZay/7X9lOzpIbG3zZiuy0Z52fQ3dxNdHI0WRdlOV2iiAyAQpiISIjZvXw3m6/eTHRGNAVfLSD/S/kkliQ6XZaIDJJCmIhIkGte34x7qZuU0hQYA1kXZzF+2XhyrsjBlah1HEVClUKYiEgQ8nR4aPhtA+4lbvb/Yz9R8VHEDo+FMeBKcpF3XZ7TJYrISVIIExEJQhuv2MieP+whYWwCY/7fGIYv8K7juL1yu9OlicgQUQgTEXGY9Vj2/Wkf7kfdjFs6jtjcWEb+x0gKbi0g49wMreMoEqYUwkREHNK1p4u6J+uoWVpD+/vtxGTH0LK5hdjcWNLPTHe6PBHxM4UwEREHdO3r4o3CN/C0ekg7I41R3x9F9ueziYrTOo4ikUIhTEQkAHpae9j96920bm5lzI/GEJMRw5gfjSHtzDSSp2gdR5FIpBAmIuJHre+14n7ETd0TdXTv6ybp1CSK7y/GFe+i4LYCp8sTEQcphImI+El9RT2bv7gZE23I+nwWBbcVkHam1nEUES+FMBGRIdJR10HdL+pInp5M5mczSf9UOsX3FZN3Ux5xeXFOlyciQUYhTETkJFhr2f/33nUcX/Cu4zjyGyPJ/GwmccPjKL6n2OkSRSRIKYSJiJyEjZdvpPG3jbjSXBTcXkD+zfkkjtc6jiJyfAphIiKD0PxuM7WP1zJ68WhcSS5yrswhc34mOVdqHUcRGRyFMBGR4/B09lnH8e/7MXGG7EuzST87nZwv5DhdnoiEKIUwEZFj6KzvZPXU1XTVdxE/Op7RPx5N3nV5xGTGOF2aiIQ4hTARkT6sx7Lv1X20VbVRcGsBsbmx5F6dS8Z5GQw7f5jWcRSRIaMQJiICdO31ruPoXuqmraqNuKI48m7KIyomirEPjXW6PBEJQwphIhLx6p+tZ+sNW/G0e0j9ZCrF9xWTfVk2UTFax1FE/Oe4IcwY82XgGWvtvgDUIyLidz1t3nUck05JInVWKsnTk8m9NpeCWwpInqp1HEUkMAYyEpYLrDbGrAOWAa9Ya61/yxIRGXqtVb3rOC7zruNY8JUCUmelkjQxifGPjHe6PBGJMMcNYdba7xpj7gY+DVwH/NwY8xzwC2vt+/4uUERkKGy+djP1T9d713G8NIv8W/NJPzvd6bJEJIINaE6YtdYaY+qAOqAbyACeN8b8yVp7lz8LFBE5EZ27O6l/tp4RXx6BcRmSpyWTMDrBu45jvtZxFBHnDWRO2FeBfwcagceBb1hru4wxUcB7gEKYiAQFay0HXj9AzZIaGn7TgO2ypJSmkH5GOiO/NtLp8kREjjCQkbBhwOettTv6NlprPcaYC/1TlojI4LTvamfD/A20rG/BleYi/9Z88m/OJ2lCktOliYj4NJA5Yd87xrbNQ1uOiMjAtWxqoe39NrIuyiIuL4744ngKvlxA7lW5uJK0jqOIBDddJ0xEQoqn00Pj7xupWVLD/r/tJ25EHJnzMzEuw6l/ONXp8kREBkwhTERCxu7nd1P15So66zqJL45n9AOjGX7dcC0lJCIhya8hzBhzAfBTwAU8bq39r6O23wzcBvQAzcBCa+0mf9YkIqHDWsu+P+8jYVQCCWMSiM2NJXlmMgW3FnjXcXQpfIlI6PLbmhzGGBfwMPAZYBJwlTFm0lG7PWutPdVaOw34EfCQv+oRkdDR1dQFz8OqCatYf956apbWAJB+ZjpTXpxC5mczFcBEJOT5cyRsNlBlra0GMMYsBy4BDo90WWsP9Nk/CdCV+EUi3Ht3vEdteS20QczpMRT9sojsy7OdLktEZMgZf61AZIy5HLjAWntj7/1rgDnW2tuP2u824OtALHCOtfY9H4+1EFgIkJubO3P58uV+qXkoNDc3k5ysteeOpn7xTf0CdAKrgE8CBu/4eRu0ntdK4tRER0sLRjpm+qe+8U390r9A9M28efPWWmtLfW1zfGK+tfZh4GFjzNXAd4FrfexTDpQDlJaW2rlz5wa0xsGorKwkmOtzivrFt0jul7bqNtyPuKldVkv3nm5mrplJyswUmOvdHsl9cyzql/6pb3xTv/TP6b7xZwirAfpeonpEb1t/lgNL/ViPiASB9p3tbLt5G3tf3gtRkPW5LApuKSB5hv6nLiKRxZam8UUAABrISURBVJ8hbDVQYowZhTd8XQlc3XcHY0xJn9OP8/EugyQiYaazoZP27e2kzkolJiuG9p3tFN1TRP5N+cQVaB1HEYlMfgth1tpuY8ztwCt4L1GxzFq70RhzP7DGWrsCuN0Ycy7QBezDx6lIEQlN1loOvPHROo7xxfHM3jIbV6KLWRtmYYw+3Sgikc2vc8KstSuBlUe13dPn9lf9+fwi4ow9K/fwwaIPaH67GVeqi/wveddxPBS8FMBERIJgYr6IhIeWzS3EZMcQmxVLT2sP1mMZ9+g4cq7OITpZbzUiIkfTO6OInDBPl4fGPzTiXuKm6a9NFN9fTPHdxWR/Ppvsy7I14iUicgwKYSIyaNZadvxgB+6lbjprO4krimPUf44i7/o8AK3lKCIyAAphIjIg1lqa324mZXoKxhia1zWTPC2Z/PJ8Mj+jZYRERAZLIUxEjqmrqYv6p+txL3XTurWVOe/NIWFMApN+M4moaL8tPysiEvYUwkTEp46aDrbft536ino8rR5S5qQw4ckJxBbEAiiAiYicJIUwETnM0+Gho7aDhOIETIyh4TcN5FyVQ8EtBd4lhUREZMgohIkIbR+04X7UTd0v6kgYn8CMf8wgNieW02tPxxXvcro8EZGwpBAmEsGa/t7Ezgd2snflXjCQdUkW+bfkH96uACYi4j8KYSIRprOxE1eiC1eii5YNLTSvbabou0Xk3ZRH/Mh4p8sTEYkYCmEiEcBay4E3D+Be4mb3r3dT8rMS8r+UT94NeeTdmEdUrCbZi4gEmkKYSBiz1lL7i1rcS9w0v9WMK8VF3k15pM9NByAqTuFLRMQpCmEiYahzdyexObEYY6h9vBbbbRn3SO86jin6tRcRCQZ6NxYJE55uD3tW7KFmSQ0H/nmA03edTkxmDFNemkL0sGit4ygiEmQUwkRCXGdjJ+4lbtzlbjprOokbGUfR3UXQ+8HGmMwYZwsUERGfFMJEQpC1lp4DPUSnRdO9t5vt920n47wMxi0ZR+Z8reMoIhIKFMJEQkj3/m7qflmHe4mbxImJTH5hMonjEjn9w9OJy49zujwRERkEfTRKpD8VFVBcDFFR3u8VFY6V0vxuM1tv3srrBa9T9eUqXCkusi7NOrxdAUxEJPRoJEzEl4oKWLgQWlu993fs8N4HKCsLSAmeDg8m2mBchobnG6h/qp6cq3LIvzWf1NLUgNQgIiL+o5EwEV8WLfoogB3S2upt97P2He1Uf6eaN0a+QeP/NgIw4qsjOL3mdCYsm6AAJiISJjQSJuLLzp2Daz9J1mPZ+8e9uJe42fPSHgAyL8o8fJoxJkOfcBQRCTcKYSK+FBZ6T0H6ah9K3R/dfO/29+hp7qHoO73rOBZqHUcRkXCmECbiy+LFR84JA0hM9LafJGstB1cfpGZJDbwEPR/24Ip3MWXlFOKL47WOo4hIhFAIE/Hl0OT7RYu8pyALC70B7CQm5fe09rB7+W5qltTQvLYZV7ILPgWeFg+ueBeJ4xKHqHgREQkFCmEi/SkrG5JPQtoei3EZmt9uZusNW0k8JZGSJSXkfjGXf6z9h65oLyISoRTCRPzA0+1hz4t7cD/sJqEkgXFLxpF6eioz/jWDlNkpWsdRREQUwkSGUkddB7WP11L7aC0duzqIGxHHsM8MA8AYQ+ocXV5CRES8FMJETpK1FvCGrB337cD9iJuM8zIo+XkJw+YPIypaE+1FROTjFMJETlD3gW7qn6mnZkkN48vHk/aJNAq/VciIr43QJHsRETkuhTCRQWre0Ix7qZv6X9bT09xD8oxkPJ0eAOKLdG0vEREZGIUwkUHwdHp4e97b9DT3kHNlDgW3FpAySxPtRURk8BTCRI6hw91BzZIamv7axPS/TycqNorJL0wmaXKSLi0hIiInRSFMxAdrLXXL6qj6ehU9zT1kXphJ975uYjJjSD873enyREQkDCiEiRyla18Xm67axL5X9pE+N51xj40jcawm2ouIyNBSCBM5iivFhafdQ8nPS8i/JR8TpfleIiIy9HQBIxGgfVc7m6/ZTNeeLqKio5j212kU3FagACYiIn6jECYRzVpL7bJaVp+ymobfNnBw3UEAfdpRRET8TiFMIlb7rnY2fHYDW2/YSvL0ZGZtmMWw84Z9tENFBRQXQ1SU93tFhVOliohIGNKcMIlY1XdV0/Rak++5XxUVsHAhtLZ67+/Y4b0PUFYW+GJFRCTsaCRMIkr7rnbad7QDMOb/jWHWhlm+534tWvRRADuktdXbLiIiMgQUwiQi9J37te2WbQDE5cWRMDrB9w/s3Dm4dhERkUFSCJOwd/Tcr5L/KTn+DxUWDq5dRERkkBTCJKztf30/q09ZTdNrTYz9n7FM+8s0Esb0M/rV1+LFkHjUBVoTE73tIiIiQ0AhTMKStRaApClJZF6Uyaz1sxhx+4iBX/errAzKy6GoCIzxfi8v16R8EREZMvp0pIQVay11T9RR+3gtU/8ylejkaCY9M+nEHqysTKFLRET8RiNhEjb6zv0yMYbupm6nSxIREemXRsIk5B0a/ar6WhW22zL2Z2O15JCIiAQ9hTAJebbHUrOkhuRpyUxYNmFgE+9FREQcphAmIclaS/0v68mcn0lMZgxT/m8KMZkxGv0SEZGQoTlhEnIOzf3acu0W3I+4AYjNjlUAExGRkKKRMAkZ1lrqnuyd+9X10dwvERGRUKQQJiFjx+IdbL97O2lnpWnul4iIhDyFMAlq1lp6DvYQnRpN3vV5xAyLIf/mfJ16FBGRkOfXOWHGmAuMMVuNMVXGmG/52P51Y8wmY8x6Y8yfjTFF/qxHQktHTQcbLtzAhos3YD2WuPw4Cm7VpSdERCQ8+C2EGWNcwMPAZ4BJwFXGmKMvXf4WUGqtnQI8D/zIX/VI6LDWUvtELatOWUVTZRPZl2U7XZKIiMiQ8+fpyNlAlbW2GsAYsxy4BNh0aAdr7V/77P8v4It+rEdCQOfuTrZct4W9K/dq7peIiIQ1c2ih4yF/YGMuBy6w1t7Ye/8aYI619vZ+9v85UGet/YGPbQuBhQC5ubkzly9f7peah0JzczPJyclOlxF0BtwvbcCtwIXApYT9RVR0vPRPfeOb+qV/6hvf1C/9C0TfzJs3b621ttTXtqCYmG+M+SJQCpzta7u1thwoBygtLbVz584NXHGDVFlZSTDX55Rj9UtHTQc7friDMf9vDK54F573PERFh3n66qXjpX/qG9/UL/1T3/imfumf033jz790NcDIPvdH9LYdwRhzLrAIuNha2+HHeiTIWGupfdI796vuiToOrjkIEDEBTEREIps/R8JWAyXGmFF4w9eVwNV9dzDGTAcexXvacrcfa5Eg01HTwdaFW71zv85MY/yy8SSOTXS6LBERkYDxWwiz1nYbY24HXgFcwDJr7UZjzP3AGmvtCuDHQDLwG2MMwE5r7cX+qkmCx5brtrD/H/sZ+99jKfiyLjshIiKRx69zwqy1K4GVR7Xd0+f2uf58fgkuHTUdRMVHEZMZQ8nPSyAKjX6JiEjE0uQb8TtrLbwMq05ZRdV/VAGQOC5RAUxERCJaUHw6UsJXR00HW7+0FV6C5DOTKfquFkUQEREBhTDxo31/3se7l72L7bRwG0z72TTN/RIREeml05HiN4mTEkk/K53S9aVwOQpgIiIifSiEyZCx1lL3VB0bPte74HZeHKeuOFVzv0RERHxQCJMh0eHuYMNFG9iyYAvde7rpbup2uiQREZGgpjlhclKstdQ/XU/VHVV4Ojy67peIiMgAKYTJSfG0edh+33aSJid5r3pfolOPIiIiA6EQJoNmraXh+QYyL8zElehiWuU04griMC6NfomIiAyU5oTJoHS4O3j34nfZ9G+bqP1FLQDxhfEKYCIiIoOkkTAZEGst9b+sp+qr3rlfY34yhoJbC5wuS0REJGQphMmAfPCdD9j5XztJ/WQqE56YoLlfIiIiJ0mnIyNNRQUUF0NUlPd7RUW/u1pr6WnvASD3mlzG/GQM0/82XQFMRERkCGgkLJJUVMDChdDa6r2/Y4f3PkBZ2RG7drg72LpwK9Gp0Ux6dhJJk5JImpQU4IJFRETCl0bCIsmiRR8FsENaW73tvay11D1dx+pTVtP05yZSZqVgrQ1woSIiIuFPI2GRZOfOY7Z31HWw7aZt7Hlxj+Z+iYiI+JlGwiJJYeGx2y0cfOsgYx7S3C8RERF/UwiLJIsXQ+KRwaojfgTV0352eMHtOVVzGPm1kbrul4iIiJ8phEWSsjIoL4eiIiyGusyrWB31FLteSad5fTMArniXw0WKiIhEBoWwSFNWRsfrW3n3wnfYsmchidPTKV1fSsq0FKcrExERiSiamB9hrLVsuHADrZtbGfPQGEZ8ZYROPYqIiDhAISxC2B6L9ViiYqIY//h4XMkuEsdp4r2IiIhTdDoyAnQ2dLL+M+up/mY1ACkzUhTAREREHKYQFub2v76fNdPX0PRaE0mn6Ir3IiIiwUIhLExZa/nwvz/k7bPfJiouihlvzCDvhrxBrR0pIiIi/qM5YWGqfXs7H3z7A4bNH8aEJycQkx4zqLUjRURExL80EhZmOtwdACSMSmDGqhlM/t1kbwCDAa0dKSIiIoGhEBZG6p6p482SN6mvqAcg+dRkjOlz+YnjrB0pIiIigaMQFgZ62nvYevNWtlyzhZTSFNLPSfe94/HWjhQREZGAUQgLcW3b23jrjLeofbSWkXeNZOqfpxKXF+d7Zx9rR5KY6G0XERGRgFIIC3HNa5tpf7+dyb+fzJgHxhAVfYx/0j5rR2KM93t5uSbli4iIOECfjgxBtsdyYPUB0k5LI/uybNLnpRMzLGZgP1xWptAlIiISBDQSFmI66zt559Pv8PZZb9NW3QYw8AAmIiIiQUMjYSGk6R9NbLpiE917uxlXPo6E0QlOlyQiIiInSCEsRHz43x/y/p3vkzAqgSkrp5A8NdnpkkREROQkKISFiK7GLrIuyWLCsglEp+mfTUREJNTpr3kQa17fTM/BHtI+mcao+0ZBFEdefFVERERClibmB6m6p+pYd9o63vvKe1hrMS6jACYiIhJGNBIWZHrae6j6ShW1j9WSPi+dSb+apPAlIiIShhTCgkjX3i7eOe8dmtc1U/jtQorvLz72xVdFREQkZOkvfKBVVEBxMURFeb9XVBzeFJ0RTdKpSUxeMZnRPxytACYiIhLGNBIWSBUVsHAhtLZ67+/Ygeemm9nxQgrDHzqPhOIEJj450dkaRUREJCA01BJIixZ9FMCADjJY33Y/O36XSsNzDQ4WJiIiIoGmEBZIO3cevtnEqaylnANMZAL/ReFdhQ4WJiIiIoGmEBZIw4YBsIc5vM1PcNHGDG5leOYahwsTERGRQNOcMAek8Q4jeJ5iniaaViDT6ZJEREQkwDQSFiAd7g627LmBbuKJpp2xPNIbwIC9e50tTkRERAJOISwADqw+wNpZa9lt5tHC6I/vUKj5YCIiIpFGIczP6p+t560z3yIqNooZP9xPWuL2I3dITITFix2pTURERJyjEOZHNQ/XsLlsM6lzUpmxagbJ37oCysuhqAiM8X4vL4eyMqdLFRERkQDTxHw/GjZ/GCM+GOG9+n1sb94tK1PoEhEREY2EDblaeP+b72M9loTiBMY+OPajACYiIiLSS+lgCO2r3Ac3Q+1jtbRVtzldjoiIiAQxv4YwY8wFxpitxpgqY8y3fGw/yxizzhjTbYy53J+1+Jv7UTfrz1sP6TDjzRkkjk10uiQREREJYn4LYcYYF/Aw8BlgEnCVMWbSUbvtBBYAz/qrjkCo/m41227eRsa5GfAwJJYogImIiMix+XMkbDZQZa2tttZ2AsuBS/ruYK3dbq1dD3j8WIffDTtvGCO/MZJTXzwVkp2uRkREREKBsdb654G9pxcvsNbe2Hv/GmCOtfZ2H/s+CbxorX2+n8daCCwEyM3Nnbl8+XK/1DwoO4B1wKVHNjc3N5OcrCR2NPWLb+qX/qlvfFO/9E9945v6pX+B6Jt58+attdaW+toWEpeosNaWA+UApaWldu7cuYEvoqICFi2CnTvZk30hmw7cgSs9nln3zSImPebwbpWVlThSX5BTv/imfumf+sY39Uv/1De+qV/653Tf+DOE1QAj+9wf0dsWeioqYOFCbGsrH3IF1bsXkmyqmfztZGLSP+F0dSIiIhKC/DknbDVQYowZZYyJBa4EVvjx+fxn0SJsaytb+QbV3Ew2rzHdfpn4hz72gU8RERGRAfHbSJi1ttsYczvwCuAClllrNxpj7gfWWGtXGGNmAb8DMoCLjDH3WWtP8VdNJ2znTgyQyhbiqaeIpzG97SIiIiInwq9zwqy1K4GVR7Xd0+f2arynKYNW69ZW2pLPIfPgn8nnf4/cOGyYM0WJiIhIyNMV849hX+U+1p2+jvdabsITGp9hEBERkRChENaP2idrWf/p9cQOj2Wq5+tE0f3xnfbuDXxhIiIiEhYUwvqqqMAWFVNtbmTrdVtJG9fG9Nenk1AU43v/wsLA1iciIiJhQyHskN7LULBzB90kk8eLTKn+N2Jeeg4WL4bEo5YiSkz0touIiIicAE10OmTRImhtxQAl/A/gwbT1tm/f/tE+O3d6R8AWL4ayMufqFRERkZCmEHZIn8tNmL5LWR5qLytT6BIREZEho9ORh/R3uQldhkJERET8QCFMRERExAEKYYf0d7kJXYZCRERE/EAh7JD+Ljehy1CIiIiIHyiEHaLLUIiIiEgAKYQdUlYG5eVQVATGeL+Xl+sTkSIiIuIXukRFX7oMhYiIiASIRsJEREREHKAQJiIiIuIAhTARERERByiEiYiIiDhAIUxERETEAQphIiIiIg5QCBMRERFxgEKYiIiIiAMUwkREREQcoBAmIiIi4gBjrXW6hkExxjQAO5yu4xiygEaniwhC6hff1C/9U9/4pn7pn/rGN/VL/wLRN0XW2mxfG0IuhAU7Y8waa22p03UEG/WLb+qX/qlvfFO/9E9945v6pX9O941OR4qIiIg4QCFMRERExAEKYUOv3OkCgpT6xTf1S//UN76pX/qnvvFN/dI/R/tGc8JEREREHKCRMBEREREHKIQNkDHmAmPMVmNMlTHmWz62f90Ys8kYs94Y82djTFGfbT3GmLd7v1YEtnL/G0DfLDDGNPTpgxv7bLvWGPNe79e1ga3cvwbQLz/p0yfbjDFNfbaF7TFjjFlmjNltjHm3n+3GGPOz3n5bb4yZ0WdbOB8vx+uXst7+2GCMed0YM7XPtu297W8bY9YErurAGEDfzDXG7O/zO3NPn23H/D0MZQPol2/06ZN3e99XhvVuC9tjxhgz0hjz196/yRuNMV/1sU9wvM9Ya/V1nC/ABbwPjAZigXeASUftMw9I7L19C/DrPtuanX4NDvfNAuDnPn52GFDd+z2j93aG068pUP1y1P5fBpZFyDFzFjADeLef7Z8F/g8wwGnAm+F+vAywXz5x6PUCnznUL733twNZTr8GB/tmLvCij/ZB/R6G2tfx+uWofS8C/hIJxwyQB8zovZ0CbPPxdyko3mc0EjYws4Eqa221tbYTWA5c0ncHa+1frbWtvXf/BYwIcI1OOW7fHMP5wJ+stXuttfuAPwEX+KnOQBtsv1wF/CoglTnMWvsasPcYu1wCPG29/gWkG2PyCO/j5bj9Yq19vfd1Q2S9xwzkmOnPybw/Bb1B9kskvcfUWmvX9d4+CGwGCo7aLSjeZxTCBqYA+LDP/V18/B+0rxvwJuxD4o0xa4wx/zLGfM4fBTpooH1zWe+Q7/PGmJGD/NlQNODX1nvqehTwlz7N4XzMHE9/fRfOx8tgHf0eY4E/GmPWGmMWOlST0043xrxjjPk/Y8wpvW06ZgBjTCLeIPFCn+aIOGaMMcXAdODNozYFxftMtL8eOFIZY74IlAJn92kustbWGGNGA38xxmyw1r7vTIWO+F/gV9baDmPMl4CngHMcrimYXAk8b63t6dMW6ceM9MMYMw9vCDujT/MZvcdLDvAnY8yW3lGSSLEO7+9MszHms8DvgRKHawomFwH/tNb2HTUL+2PGGJOMN3jeYa094HQ9vmgkbGBqgJF97o/obTuCMeZcYBFwsbW241C7tbam93s1UIk3lYeL4/aNtXZPn/54HJg50J8NYYN5bVdy1GmCMD9mjqe/vgvn42VAjDFT8P4OXWKt3XOovc/xshv4Hd7TcBHDWnvAWtvce3slEGOMyULHzCHHeo8Jy2PGGBODN4BVWGt/62OXoHifUQgbmNVAiTFmlDEmFu8BfcQn1owx04FH8Qaw3X3aM4wxcb23s4BPApsCVrn/DaRv8vrcvRjv+XmAV4BP9/ZRBvDp3rZwcNx+ATDGTMA7+fONPm3hfswczwrg33s/vXQasN9aW0t4Hy/HZYwpBH4LXGOt3danPckYk3LoNt5+8flpuXBljBlujDG9t2fj/du2hwH+HoYzY0wa3jMzf+jTFtbHTO+x8Atgs7X2oX52C4r3GZ2OHABrbbcx5na8/xAuvJ9i22iMuR9YY61dAfwYSAZ+0/tesNNaezEwEXjUGOPB+8bwX9basPmDOsC++Yox5mKgG+8k0gW9P7vXGPN9vG+UAPcfNVwesgbYL+D9o7Dc9n4sp1dYHzPGmF/h/TRbljFmF/A9IAbAWvsIsBLvJ5eqgFbgut5tYXu8wID65R4gE1jS+x7Tbb0LD+cCv+ttiwaetda+HPAX4EcD6JvLgVuMMd1AG3Bl7++Uz99DB16CXwygXwAuBf5orW3p86Phfsx8ErgG2GCMebu37TtAIQTX+4yumC8iIiLiAJ2OFBEREXGAQpiIiIiIAxTCRERERBygECYiIiLiAIUwEREREQcohImIiIg4QCFMRERExAEKYSISsYwxs3oXlo/vvYr4RmPMZKfrEpHIoIu1ikhEM8b8AIgHEoBd1tr/dLgkEYkQCmEiEtF61xRcDbQDn7DW9jhckohECJ2OFJFIl4l33dcUvCNiIiIBoZEwEYloxpgVwHJgFJBnrb3d4ZJEJEJEO12AiIhTjDH/DnRZa581xriA140x51hr/+J0bSIS/jQSJiIiIuIAzQkTERERcYBCmIiIiIgDFMJEREREHKAQJiIiIuIAhTARERERByiEiYiIiDhAIUxERETEAQphIiIiIg74/y+2cyD4DLt6AAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# coefficients of determinations\n",
        "ssr = np.sum((Y-pred)**2)\n",
        "sse = np.sum((Y-np.mean(Y))**2)\n",
        "r2 = 1-(ssr/sse)\n",
        "print(r2)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "B4PSnA9q0owo",
        "outputId": "0acc1533-16d8-4e4b-b633-20530919036c"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.9598634601190857\n"
          ]
        }
      ]
    }
  ]
}