path: root/jax_optimization_tutorial/linear-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3387e989",
   "metadata": {},
   "source": [
    "# Linear Regression with JAX\n",
    "Source: https://danielrothenberg.com/blog/2020/Sep/jax-first-steps-pt1/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "36893279",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import jax.numpy as jnp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "565bff63",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax import random\n",
    "\n",
    "def make_key():\n",
    "    \"\"\" Helper function to generate a key for jax's parallel PRNG \n",
    "    using standard numpy random functions. \n",
    "\n",
    "    \"\"\"\n",
    "    seed = np.random.randint(2**16 - 1)\n",
    "    return random.PRNGKey(seed)\n",
    "\n",
    "n = 100\n",
    "rands = random.uniform(make_key(), shape=(n, ), minval=-1, maxval=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "11fa0b53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]\n",
      " [1. 1.]]\n",
      "[[-0.49573565  1.        ]\n",
      " [-0.75743556  1.        ]\n",
      " [-0.10826278  1.        ]\n",
      " [-0.06045723  1.        ]\n",
      " [-0.6449168   1.        ]\n",
      " [-0.73400044  1.        ]\n",
      " [-0.37836146  1.        ]\n",
      " [-0.01198149  1.        ]\n",
      " [ 0.61107445  1.        ]\n",
      " [ 0.29680896  1.        ]\n",
      " [ 0.9514713   1.        ]\n",
      " [-0.94317484  1.        ]\n",
      " [ 0.5943241   1.        ]\n",
      " [-0.00902605  1.        ]\n",
      " [ 0.28311515  1.        ]\n",
      " [-0.1576097   1.        ]\n",
      " [ 0.50905204  1.        ]\n",
      " [ 0.40418315  1.        ]\n",
      " [-0.13876987  1.        ]\n",
      " [-0.10519576  1.        ]\n",
      " [ 0.4820497   1.        ]\n",
      " [ 0.16001129  1.        ]\n",
      " [-0.78287053  1.        ]\n",
      " [-0.83056164  1.        ]\n",
      " [-0.45683146  1.        ]\n",
      " [ 0.65588903  1.        ]\n",
      " [-0.30958128  1.        ]\n",
      " [-0.57034206  1.        ]\n",
      " [ 0.11413717  1.        ]\n",
      " [ 0.30006003  1.        ]\n",
      " [ 0.3948493   1.        ]\n",
      " [ 0.93246984  1.        ]\n",
      " [ 0.08236265  1.        ]\n",
      " [ 0.5821161   1.        ]\n",
      " [ 0.01212478  1.        ]\n",
      " [ 0.9328327   1.        ]\n",
      " [ 0.09593153  1.        ]\n",
      " [ 0.71609616  1.        ]\n",
      " [ 0.5419252   1.        ]\n",
      " [-0.78604007  1.        ]\n",
      " [ 0.9995785   1.        ]\n",
      " [-0.48301148  1.        ]\n",
      " [ 0.93194747  1.        ]\n",
      " [ 0.6257863   1.        ]\n",
      " [-0.42878604  1.        ]\n",
      " [-0.11879826  1.        ]\n",
      " [-0.60145855  1.        ]\n",
      " [ 0.25006437  1.        ]\n",
      " [ 0.4595716   1.        ]\n",
      " [ 0.48950648  1.        ]\n",
      " [ 0.77740407  1.        ]\n",
      " [ 0.17917514  1.        ]\n",
      " [ 0.52525926  1.        ]\n",
      " [ 0.7628691   1.        ]\n",
      " [-0.6807344   1.        ]\n",
      " [ 0.14769888  1.        ]\n",
      " [ 0.8373101   1.        ]\n",
      " [-0.90767694  1.        ]\n",
      " [ 0.90058756  1.        ]\n",
      " [-0.4039154   1.        ]\n",
      " [-0.9033947   1.        ]\n",
      " [ 0.8437958   1.        ]\n",
      " [-0.75383854  1.        ]\n",
      " [-0.7118697   1.        ]\n",
      " [ 0.2737422   1.        ]\n",
      " [-0.9684801   1.        ]\n",
      " [ 0.20652676  1.        ]\n",
      " [ 0.01346517  1.        ]\n",
      " [-0.9524138   1.        ]\n",
      " [-0.08028603  1.        ]\n",
      " [-0.71781445  1.        ]\n",
      " [ 0.08346677  1.        ]\n",
      " [ 0.6743789   1.        ]\n",
      " [ 0.01168156  1.        ]\n",
      " [ 0.0992105   1.        ]\n",
      " [-0.10572052  1.        ]\n",
      " [-0.07284904  1.        ]\n",
      " [-0.66415024  1.        ]\n",
      " [ 0.37942672  1.        ]\n",
      " [ 0.16580987  1.        ]\n",
      " [ 0.15668988  1.        ]\n",
      " [-0.4482379   1.        ]\n",
      " [-0.9518373   1.        ]\n",
      " [-0.22508144  1.        ]\n",
      " [-0.87991786  1.        ]\n",
      " [ 0.02797484  1.        ]\n",
      " [ 0.5076406   1.        ]\n",
      " [-0.40049028  1.        ]\n",
      " [ 0.55776215  1.        ]\n",
      " [-0.44021344  1.        ]\n",
      " [ 0.17129445  1.        ]\n",
      " [-0.4132347   1.        ]\n",
      " [ 0.4849968   1.        ]\n",
      " [ 0.7109468   1.        ]\n",
      " [ 0.66096044  1.        ]\n",
      " [-0.56528807  1.        ]\n",
      " [-0.55501866  1.        ]\n",
      " [-0.38958454  1.        ]\n",
      " [-0.17273688  1.        ]\n",
      " [-0.00741458  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "# We have to use functions to update an array unlike numpy as JAX does not mutate arrays in-place\n",
    "\n",
    "x = jnp.ones((n, 2))\n",
    "print(x)\n",
    "x = x.at[:, 0].set(rands)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9a1e2fb",
   "metadata": {},
   "source": [
    "# Simple first model\n",
    "\n",
    "$$\n",
    "    y = X\\beta + \\epsilon\n",
    "$$\n",
    "\n",
    "where $\\beta = [\\beta_0, \\beta_1]$ where $\\beta_0$ a slope relating x and y abd $\\beta_1$ is an offset bias, and $\\epsilon$ is an uncorrelated noise modeled as a normal distribution, $\\mathcal{N}(0, 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "7bf67529",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fix our true slope and bias\n",
    "slope, bias = 3.0, 2.0\n",
    "beta_true = jnp.array((slope, bias))\n",
    "eps = random.normal(make_key(), shape=(n,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5a6d3ce8",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = (x @ beta_true) + eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "462c84ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f8b306678b0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAc/klEQVR4nO3dfZAlVXkG8OdhGGBQ46yyURhYd4mKsSTZ1Sk/somG9WOJWrAiBkiR+JXaaBJLygQdgn+YlBabbFU0qVjRLUWTSJAEZCRBXcFdigrlIjNZZPlaXVGEK8r4MUZkxGF588ftCz13uvt23z7dfU7f51e1tTN97+0+0z3z9ulz3nMOzQwiIhKuI5ougIiIlKNALiISOAVyEZHAKZCLiAROgVxEJHBHNnHQ4447ztavX9/EoUVEgjU/P/9DM1vbv72RQL5+/XrMzc01cWgRkWCRvDdpu5pWREQCp0AuIhI4BXIRkcApkIuIBE6BXEQkcI1krYiItM3s/g527j6I7y0u4YTJCVy49RRs2zRVy7EVyEVESprd38FFnzuApeXDAIDO4hIu+twBAKglmDtpWiE5SfJKkneTvIvky1zsV0QkBDt3H3w8iPcsLR/Gzt0Hazm+qxr5PwD4kpmdTfIoAMc62q+IiPe+t7hUaLtrpWvkJJ8K4OUAPgkAZvZLM1ssu18RkVCcMDlRaLtrLppWNgBYAPApkvtJfoLkkxzsV0QkCBduPQUT42Mrtk2Mj+HCrafUcnwXgfxIAC8E8M9mtgnAzwHM9L+J5HaScyTnFhYWHBxWRMQP2zZN4ZKzTsXU5AQIYGpyApecdWptWSssu2YnyWcC2Gdm66PvfwfAjJm9Lu0z09PTpkmzRESKITlvZtP920vXyM3s+wDuI9l7hnglgDvL7ldERPJxlbXyLgCXRRkr9wB4q6P9iojIAE4CuZndCmBVdV9ERKqnuVZERAKnQC4iEjjNtSIirdXkRFZ1UiAXkVYqMpFV6AFfTSsi0kp5J7LqBfzO4hIMTwT82f2dGktbjgK5iLRS3omsmp650AUFchFppbwTWTU9c6ELCuQi0kp5J7JqeuZCFxTIRaSV8k5k1fTMhS4oa0VEnPEt+2PbpqmBx++97lO5i1IgFxEnml63sow8Ad9naloRESfakP0RKgVyEXGiDdkfoVIgFxEn2pD9ESoFchFxog3ZH6FSZ6eIONGG7I9QOQnkJL8D4GcADgN4NGlNORFpv9CzP0LlskZ+mpn90OH+REQkB7WRi4gEzlUgNwBfJjlPcnvSG0huJzlHcm5hYcHRYUVExFXTym+bWYfkrwK4juTdZnZj/A1mtgvALgCYnp42R8cVkSH4NpQ+VL6cRyeB3Mw60f8PkrwawIsB3Jj9KRFpQshD6X3i03ks3bRC8kkkn9L7GsBrANxedr8iUg0NpXfDp/Pookb+DABXk+zt79/N7EsO9isiFdBQejd8Oo+lA7mZ3QPgNx2URURqcMLkBDoJwUZD6Yvx6Twq/VBkxGgovRs+nUcN0RcZMRpK74ZP55Fm9WcCTk9P29zcXO3HFREJGcn5pClQ1LQiIhI4BXIRkcCpjVxEaufLiMi2UCAXkVr5NCKyLdS0IiK1ShsRecEVt2Lzjj2Y3d9pqGThUiAXkVpljXzs1c4VzItRIBeRWg0a+ah5X4pTIBeRWiWNiOwXyrwvs/s72LxjDzbMXNtos5A6O0WkVvERkUlzlQD55itpOvOlaKdtleVVIBeRoZQJTL1FmvuDIZBvvhIfMl+yprHtL0PV5VXTiogU1gtMncUlGIbvpNy2aQqXnHUqpiYnQABTkxO45KxTBwY3H+YCLzKNbdXlVY1cRAorUhsdpFc7L8KHucCLTGNbdXmdBXKSYwDmAHTM7PWu9iviq6bbaJvUdCBtYi7w/ut92vPW4qr5Tq5moarL67Jp5d0A7nK4PxFvuWpaCEV/dsZTJ8YT31fXogp1zwWedL2vmu/gjS+aytUsVHV5ndTISZ4I4HUAPgTgPS72KeIzl00LRTTxFJDUUTc+RowfQSw/9sQ02HUuqlD3XOBp13vv3Qu4aWZL4+V11bTyEQDvBfCUtDeQ3A5gOwCsW7fO0WFFmtFE00JTmRpJQWz5sGHNseM49qgj0Vlcwhi5ovOujiamYdrWh+XieldZ3tJNKyRfD+BBM5vPep+Z7TKzaTObXrt2bdnDijQqrQmhyqaFpjI10oLV4sPLjzcZHI4WqPGlicn1QJ0mrncRLtrINwM4g+R3AHwWwBaSn3GwXxFvNbFeY1MdjFlBzIc0wH5V9F/4tD5nktKB3MwuMrMTzWw9gHMB7DGz80uXTMRjw+Y/l9FUrTAriDWdvZKkiptLE9e7COWRiwypzjZaoBtQhxkFWVZWR13aMPu0m0sdnbVV3Vzqvt5FOA3kZnYDgBtc7lNEuppctT0tiBW5udTVWdtEjnnTVCMXCYhvtcIiN5e6UjabenJpkgK5iJSSdHNJakKpqz296JNLG0boKpCLiFNpTShPnRjH4tLyqvdX0eSR98nFh1kUXVAgl5HUhlpYE/Kct7QmlGPGj8DE+JhXTR5NjdB1TYFcRk5bamF1y3vesgYQfficjYVuoFXfcH1MnxyGArmMnLbUwuqW97xlZY0U6azNunH0ylM2wJfNcPHlyU4LS8jIaUstrG55z5urUZBpN44PXHOHs5GbZcrq0wyYCuQycnyfN8NXec+bq1GQqU00S8vORm6WKatP0xOoaUVGzijmGbtQ5Ly5yHdPa/ZIM+wT1bBl9enJTjVyGTm+z5vhq7rPW1qzx5pjm13UYtDxmniyU41cRlLdIyR96RQrW5Zhz9swx0wb2APAiycqn57sFMhFKuZTumMTZRnmmP2B/8PnbEzMV2/yxtjk3Df9aGaD3+XY9PS0zc3N1X5cEReK1i4379iT2NY7NTmRa5kwl5ooS9Fj9gd+oFvTVfMXQHLezKb7t6tGXiGfHqfFjWFqlz51ijVRlqLHVJ5/cersrIhPOabizjApZz51ijVRlqLH9OnGFwoF8or4lGMq7gwTZHxaJqyJshQ9pk83vlC4WHz5GJJfI/l1kneQ/GsXBQudahXtNEyQ8SndsYmyFD2mTze+ULhoI38EwBYze4jkOID/IflFM9vnYN/BGsVVSkIzTB/GsClnPi0I0URZihzTp2yQUJQO5NZNe3ko+nY8+ld/KoxnfMoxldWGTcNTkKmHTze+EDhJPyQ5BmAewLMBfNTM3pfwnu0AtgPAunXrXnTvvfeWPq7vlLXiL59SAkXyqjT90MwOA9hIchLA1SRfYGa3971nF4BdQDeP3MVxfadahb/UhyFt4jSP3MwWSe4FcDqA2we9X6QpdfVhtOmprE0/S9u4yFpZG9XEQXICwKsB3F12vyJVqiMzok1jCdr0s7SRizzy4wHsJXkbgFsAXGdm/+1gvyKVKZOGN7u/g8079mDDzLXYvGNPajBr01iCtJ/lgituzTwHUg8XWSu3AdjkoCwitRqmD6NItktae3uRObZ9kdV3oDVPm6eRnRK0vLVjV4rUstPa2wkEV4Md1HcQ6pNGWyiQS7CaaLctku1y4dZTwIT3GlA46NV9w+qX1KfQTxk/zVEgl2A10QZdZIj+tk1TqSPjigQ9Hzoa430KaTRquTkK5BKsJnLBi2S7zO7vYIxJdfJiQc+XTtNtm6Zw08wWfOScjZoLxTOaj1yC1cR8NnmH6Pdq0YcTRk4XDXq+DV7SNAX+USCXYDU1n02ebJekWjQAjJGFZxv0cQI2jVr2i5pWJFg+TQ/bL622/JhZ4fJpWlcZRDVyCVpSzdCHoeQua9Gj1pThw/ULjQK5tIovK9a7bvYZlaYMX65faNS0Iq1SNsPDVb62z80+PvMlQyc0qpFLq5TJ8HBdGwy9Ft1EE4dvGTqhUI1cWqXMwr2qDT6hqUFIWnh5OArk4p0yzRtlMjxUG3xCUzc1ZegMR00r4pWyzRtlMjx8zNduSlM3tVHL0HFFgVy8klUTLLIK+zB/+Fow+wlN3tRC71togosVgk4iuZfknSTvIPluFwWT0ZQ1h3fVs/4p0+QJauIIi4sa+aMA/sLM/pfkUwDMk7zOzO50sG8ZMWk1QaB4M8swWRdtrg0WOR9q4giLixWCHgDwQPT1z0jeBWAKgAK5FJbUvBEX73DLCjIaWLLSMOejzTe1tnGatUJyPbrLvt2c8Np2knMk5xYWFlweVlokz7zXncUlXHjl1zNT45RKuJLOR7s5C+QknwzgKgAXmNn/9b9uZrvMbNrMpteuXevqsNJCvXmvs4L58uGV08P2ByWlEq6k89FuTgI5yXF0g/hlZvY5F/sUybO8WFw8KGlgyUo6H+3mImuFAD4J4C4z+/vyRZKqNL3uY1G9Zpa84kFJWRcr6Xy0m4uslc0A/hDAAZK3Rtv+ysy+4GDf4kionX/bNk1h5+6DqZksPf1Bqe1ZF++fPYDLb74Ph80wRuK8l5yED25Lv+m1/XyMOlrCUlRVm56etrm5udqPO8o279iTGAynJidw08yWx7/3cS7o/psQAIwfQTz5mCOx+PCyN+V0Le1avH/2AD6z77ur3n/+S9dlBnMJH8l5M5vu366RnSMiT2eXr7X2Ntcm04J11rW4/Ob7Evd1+c33KZCPKAXyEZFnyLWL4fFVaWNOc1awzroWSQs6A0jdLu2n2Q9HRJ7OLqWo1SsrWGddizEy8bW07dJ+CuQjIs88IiGmqIWWiROXFayzrsV5Lzkp8bW07dJ+aloZIYOaJ0Kb/e/9swdw2b7voteg4Eubfl5ZzV1Z16L3sxXJWpF2UyCXx4XUqTi7v7MiiPf40qafR55gnXYtPrjtVAVueZwC+YgZlF4YSqfizt0HVwXxnlDa9AcF61CuhTRPgXyE+JpeOIysYO1zm34/BWtxQZ2dKULuREvTphnwsoL1zx95tFXXTWQQBfIETa0gXrU2pRemTah1BIHFpeVWXTeRQRTIE7Sp5hoXYnphmqR0yjXHjuOxvobzrOvWxqcuGU1qI0/QppprXGjphYP0ty9vmLk28X1J161N/QUiqpEnaFPNNa7tiwsXuW5tfeqS0aQaeYK21Vzj2pwlUeS6tfWpS0aTauQJ2l5zbasi1y2t9m6A2sslOJqPXEZS0hzncRPjY7p5i3fS5iN3tWbnpSQfJHm7i/2JVC1ee0+i9nIJiaumlU8DON3RvqQGSr3rBvObZrYgbfJXtZdLKJwEcjO7EcCPXexLqtfWAU/DamuWkoyO2jo7SW4nOUdybmFhoa7DSgKl3q104dZTMH7Eynr5+BFsRZaSjIba0g/NbBeAXUC3s7Ou48pqSr1L0N++UnCxHR8XrZbRofTDEaSmhJV27j6I5cMr6xbLhy33E0qRpir1TUgVFMhHUJ71O0dJ2SeUvE1V6puQqrhKP7wcwFcBnELyfpJvd7FfqYYGPK1U9gkl741AfRNSFSdt5GZ2nov9SH3aPFS/qLJTMmStvRmnvgmpippWZOSVfULJ21SlvgmpiibNqoEyGvxX5gkl76LVbZ6MTZqluVYqljSnB9GdnGlKQX3k6KYuZaTNtaIaecWSOrh6t84iixkoALSD+iakCmojr9igjqw8WQtKWxORLArkGVwM3sjTkTUo2CttTUSyKJCncFULTlvtPW5QsFfamohkURt5iqxacJE2znhGQ2dx6fGOzp48WQuTx47jJw8vr9ruIm2tbNv7MJ9vS3t/kZ+jLT+z+EmBPIXLWnC8g6voH/Ts/g4e+sWjq7aPj5Wfna/sSvLDfL4tq9cX+Tna8jOLv9S0kqKqwRu9xQy+veN1uGlmy8A/5J27D2L5sdUpok866sjSQaBs2/swn29Le3+Rn6MtP7P4SzXyFE0M3kiqrac9Afx0aXVTS5H9bts0lbrvzuISZvd3Bt4ohnlqaUt7f5Gfoy0/s/hLNfIUdU8slda5OnnseOL78z4ZZHXaZu0jT8fuME8tbRmmXuTnaMvPLP5SIM9QtBmkjLTHbzOUmnI267E+K6NmafkwLrji1sy0y2Gmw23LFLpFfo62/MziLwVyT2Q1oZR5Msh6rO89dWTJSrsc5qmlLVPoFvk52vIzi78014onNu/YkzgV6tTkBG6a2VLpftPe47IcIlJe2lwrrhaWOJ3kQZKHSM642OeoqerxO89+8wxaUseciL9KZ62QHAPwUQCvBnA/gFtIXmNmd5bd9yjJOxVqFfvtH7SURB1zIv4q3bRC8mUAPmBmW6PvLwIAM7sk7TNtbFppy8i9pGl3J8bH1KYr4oEqp7GdAnBf7Pv7AbwkoQDbAWwHgHXr1jk4rD9CH7nXfxN644umsPfuheBvSiKjorYBQWa2C8AuoFsjr+u4dXA1L0sTkm5CV813VAMXCYiLzs4OgJNi358YbRsZae3KIXQQavi4SPhc1MhvAfAckhvQDeDnAvgDB/sNwuz+zqoZDXtC6CDU8PGutvRxyGgqHcjN7FGSfw5gN4AxAJea2R2lSzaAL394O3cfTAziBIIYuXfC5ETiE0UINyFXQu/jEHGSR25mXzCz55rZr5nZh1zsM4tPS5+l1VwNYQQBDR9X85KEL8gh+j794aXVXKcCqdFq+LialyR8QU5j69MfXhPT3bo26iu7q3lJQhdkjdynaUF9rNG6WDR6lKh5SUIXZI38tOetxWX7vlt47cuq+FSjVcddcVVNjyBSl+AC+ez+Dq6a76wI4gTwxhf5E0ybFPLgpCb5dDMWKSq4ppWkQGUA9t690EyBPONT/4GI1CO4GrkvgcqXPPZ+6rgTGT3B1ch96Oj0KY+9nzruREZPcIG8yUDVywa54Ipbvcljj+s9JSwtH8YYCcCPLBoRqVZwTStNZRgkzdPdr8l26P7yHTZ7/AanIC7SbsEE8qbbpJM6Wfs12Q6tbBWR0RVEIPchN3pQbbvpdmhfOoFFpH5BtJH7MLdKVm3bh3ZoHzqBRaQZQQRyH2qbaZ2sHzlnI26a2dJ484WyVURGVxBNKz7kRvs+jNv38olIdWhW//KZ09PTNjc3l/v9WtldRAQgOW9m0/3bSzWtkHwTyTtIPkZy1c5d8XGGQRERX5RtWrkdwFkAPu6gLJk0qZGISLJSgdzM7gIARqMIRUSkfrVlrZDcTnKO5NzCgmYqFBFxZWCNnOT1AJ6Z8NLFZvb5vAcys10AdgHdzs7cJRQRkUwDA7mZvaqOgoiIyHCCGBAkIiLpyqYfvoHk/QBeBuBakrvdFEtERPIqm7VyNYCrHZVFRESGEMQQ/SKanu5WRKRuQQbytGDtw3S3IiJ1Cy6QZwVrLa4gIqMouECeFayLTnerZhgRaYPg0g+zgnWRxRV6NfvO4hIMT9TsZ/d3XBZXRKRywQXyrGBdZHEFH1YdEhFxIbhAnhWsi0x368OqQ3nM7u9g84492DBzLTbv2KMnBhFZJbg28kEr4fRPd9sLhP3v9WHVoUGUhSMieQSxQtCwslYWAuD9qkObd+xJvNlMTU7gppktDZRIRJpUyQpBvhuUjuj7qkOhNP+ISLOCa1opYlAgLLrqUN3piiE0/4hI81pdIy+SjjhIE+mKRbJwRGR0tTqQuwyETaQrhtD8IyLNa3XTyqAMlyKaaq/WotMiMkirAzngLhCqvVpEfFV2YYmdJO8meRvJq0lOOirX0KoaQKP2ahHxVdk28usAvMDMfgPANwBcVL5Iw6uyQ1Lt1SLiq7IrBH059u0+AGeXK045VU9jq/ZqEfGRy6yVtwH4YtqLJLeTnCM5t7Cw4PCwT9AAGhEZRQMDOcnrSd6e8O/M2HsuBvAogMvS9mNmu8xs2sym165d66b0fVzmjYuIhGJg04qZvSrrdZJvAfB6AK+0JiZuiblw6ymJ86eoQ1JE2qxUGznJ0wG8F8ArzOxhN0Uansu8cRGRUJSa/ZDkIQBHA/hRtGmfmb1j0Ofqmv1QRKRN0mY/LJu18uwynxcRkfJaPdeKiMgoUCAXEQmcArmISOAUyEVEAtfImp0kFwDcm+OtxwH4YcXFGZbKVpyv5QJUtmH4Wi6gvWV7lpmtGlHZSCDPi+RcUqqND1S24nwtF6CyDcPXcgGjVzY1rYiIBE6BXEQkcL4H8l1NFyCDylacr+UCVLZh+FouYMTK5nUbuYiIDOZ7jVxERAZQIBcRCVzjgZzkm0jeQfIxkqkpOSRPJ3mQ5CGSM7HtG0jeHG2/guRRDsv2NJLXkfxm9P+ahPecRvLW2L9fkNwWvfZpkt+OvbaxrnJF7zscO/Y1se1Nn7ONJL8aXffbSJ4Te835OUv73Ym9fnR0Hg5F52V97LWLou0HSW4tW5aC5XoPyTujc/QVks+KvZZ4bWss21tILsTK8Mex194cXf9vknxzzeX6cKxM3yC5GHut6nN2KckHSd6e8jpJ/mNU9ttIvjD2WrlzZmaN/gPw6wBOAXADgOmU94wB+BaAkwEcBeDrAJ4fvfYfAM6Nvv4YgHc6LNvfAZiJvp4B8LcD3v80AD8GcGz0/acBnF3BOctVLgAPpWxv9JwBeC6A50RfnwDgAQCTVZyzrN+d2Hv+FMDHoq/PBXBF9PXzo/cfDWBDtJ+xGst1Wux36Z29cmVd2xrL9hYA/5Tw2acBuCf6f0309Zq6ytX3/ncBuLSOcxbt/+UAXgjg9pTXX4vucpgE8FIAN7s6Z43XyM3sLjM7OOBtLwZwyMzuMbNfAvgsgDNJEsAWAFdG7/sXANscFu/MaJ959302gC9a9YtsFC3X43w4Z2b2DTP7ZvT19wA8CKCa9f9SfncyynwlgFdG5+lMAJ81s0fM7NsADkX7q6VcZrY39ru0D8CJjo5dumwZtgK4zsx+bGY/AXAdgNMbKtd5AC53dOyBzOxGdCtyac4E8K/WtQ/AJMnj4eCcNR7Ic5oCcF/s+/ujbU8HsGhmj/Ztd+UZZvZA9PX3ATxjwPvPxepfnA9Fj1EfJnl0zeU6ht0Fr/f1mnvg2Tkj+WJ0a1ffim12ec7SfncS3xOdl5+ie57yfLbKcsW9HSsXN0+6tq7kLdsbo+t0JcmTCn62ynIhaobaAGBPbHOV5yyPtPKXPmelFpbIi+T1AJ6Z8NLFZvb5OsqQJqts8W/MzEim5mpGd9ZTAeyObb4I3WB2FLq5o+8D8Dc1lutZZtYheTKAPSQPoBukSnF8zv4NwJvN7LFo89DnrK1Ing9gGsArYptXXVsz+1byHirxXwAuN7NHSP4Juk80W2o8/iDnArjSzA7HtjV9zipTSyC3AQs459ABcFLs+xOjbT9C9/HkyKgm1dvupGwkf0DyeDN7IAo6D2bs6vcBXG1my7F992qmj5D8FIC/rLNcZtaJ/r+H5A0ANgG4Ch6cM5K/AuBadG/m+2L7HvqcpUj73Ul6z/0kjwTwVHR/t/J8tspygeSr0L1BvsLMHultT7m2roLSwLKZ2Y9i334C3b6R3md/t++zN9RVrphzAfxZfEPF5yyPtPKXPmehNK3cAuA57GZbHIXuRbrGuj0Fe9FtmwaANwNwWcO/Jtpnnn2vao+LAlmvXXobgMTe7CrKRXJNr1mC5HEANgO404dzFl3Dq9FtL7yy7zXX5yzxdyejzGcD2BOdp2sAnMtuVssGAM8B8LWS5cldLpKbAHwcwBlm9mBse+K1dVSuvGU7PvbtGQDuir7eDeA1URnXAHgNVj6lVlquqGzPQ7fT8KuxbVWfszyuAfBHUfbKSwH8NKq4lD9nVfbi5vkH4A3otgk9AuAHAHZH208A8IXY+14L4Bvo3kEvjm0/Gd0/rkMA/hPA0Q7L9nQAXwHwTQDXA3hatH0awCdi71uP7l31iL7P7wFwAN1g9BkAT66rXAB+Kzr216P/3+7LOQNwPoBlALfG/m2s6pwl/e6g21xzRvT1MdF5OBSdl5Njn704+txBAL/n+Hd/ULmuj/4meufomkHXtsayXQLgjqgMewE8L/bZt0Xn8hCAt9ZZruj7DwDY0fe5Os7Z5ehmYC2jG9PeDuAdAN4RvU4AH43KfgCxLL2y50xD9EVEAhdK04qIiKRQIBcRCZwCuYhI4BTIRUQCp0AuIhI4BXIRkcApkIuIBO7/AbD57RrvT+FXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x[:,0], y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "547a51fd",
   "metadata": {},
   "source": [
    "# Analytical Model Fitting\n",
    "\n",
    "\n",
    "If the problem is simple enough, we can find a closed form analytical solution to the problem. We seek to find the $\\beta$ that predicts the relationship between x and y. If we use $\\mathcal{L}(\\beta)=||\\mathbf{X}\\beta-\\mathbf{Y}||$, the analytical solution is:\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial \\mathcal{L}(\\beta)}{\\partial \\beta} = -2 \\mathbf{Y}^T\\mathbf{X}+2\\beta^T\\mathbf{X}^T\\mathbf{X}\n",
    "$$\n",
    "\n",
    "And solving for 0 gives:\n",
    "\n",
    "$$\n",
    "\\beta = (\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{Y}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "61ad3658",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.9267178 1.855265 ]\n"
     ]
    }
   ],
   "source": [
    "beta_ols = jnp.linalg.inv(x.T @ x) @ (x.T @ y)\n",
    "print(beta_ols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "2430ba06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7jElEQVR4nO3dd3gUVdv48e9JgYQSOghBSJAaQEBpioIKgiACIj7CY0NRFJVXEaQl78/yEhIEQVAsYEGsCAIqFpD2PIogAqFICYSe0EsCgSSknN8fuwmbsJvsZie7s8n9uS4uktnZmTOzcO+ZM/e5R2mtEUII4bv8vN0AIYQQ7pFALoQQPk4CuRBC+DgJ5EII4eMkkAshhI8L8MZOa9asqcPCwryxayGE8FmbN28+o7WuVXC5VwJ5WFgYmzZt8sauhRDCZymlDttbLkMrQgjh4ySQCyGEj5NALoQQPs4rY+T2ZGZmkpiYSHp6urebYgpBQUHUr1+fwMBAbzdFCGFypgnkiYmJVK5cmbCwMJRS3m6OV2mtOXv2LImJiYSHh3u7OUIIkzNNIE9PT5cgbqWUokaNGpw+fdrbTRFCOGlpXBJTl8dzLDmNelWDeaVXMwa0C/XIvk0TyAEJ4jbkXAjhO5bGJTFh8Q7SMrMBSEpOY8LiHQAeCeaG3OxUSlVVSi1SSu1RSu1WSt1ixHaFEMIXTF0enxfEc6VlZjN1ebxH9m9Uj3wm8KvWepBSqhxQwaDteszZs2fp3r07ACdOnMDf359atSwTqDZu3Ei5cuW82TwhhIkdS05zabnR3A7kSqkqQFdgKIDW+gpwxd3telqNGjXYunUrAK+99hqVKlVizJgxea9nZWUREGCqkSghhEnUqxpMkp2gXa9qsEf2b0RkCgdOA58qpdoAm4EXtdaXDNi2Vw0dOpSgoCDi4uLo0qULISEh+QJ8q1atWLZsGWFhYXzxxRfMmjWLK1eu0KlTJ9577z38/f29fARCCE94pVezfGPkAMGB/rzSq5lH9m/EGHkAcBPwvta6HXAJGF9wJaXUcKXUJqXUpqKyMZRSJfKnOBITE/nzzz+ZPn26w3V2797NggULWLduHVu3bsXf358vv/yyWPsTQvieAe1CiRnYmtCqwSggtGowMQNb+1TWSiKQqLX+y/r7IuwEcq31HGAOQPv27X3mQaEPPvhgkT3rVatWsXnzZjp06ABAWloatWvX9kTzhBAmMaBdqMcCd0FuB3Kt9Qml1FGlVDOtdTzQHdjl5jbdbZZhKlasmPdzQEAAOTk5eb/nzkLVWvP4448TExPj8fYJIYRRtVZGAl8qpbYDbYHJBm3XVMLCwtiyZQsAW7Zs4eDBgwB0796dRYsWcerUKQDOnTvH4cN2q00KIYThDEnD0FpvBdobsS0ze+CBB5g/fz4tW7akU6dONG3aFICIiAgmTZpEz549ycnJITAwkNmzZ9OwYUMvt1gIURYobwxjtG/fXhd8sMTu3btp0aKFx9tiZnJOhBC2lFKbtdbXdJqljK0QQvg4meEihCjVvFnMylMkkAshSi1Xiln5csCXoRUhRKnlbDGr3ICflJyG5mrAXxqX5MHWFp8EciFEqeVsMStvVy90lwRyIUSp5ahoVcHl3q5e6C4J5DYSExPp378/TZo04YYbbuDFF1/kyhVLIce1a9fSt2/fa96zbNky2rVrR5s2bYiIiODDDz+8Zp2MjAx69OhB27ZtWbBgAU899RS7dlkmv06eXCrnTglhCq/0akZwYP4SG/aKWTkb8M1KArmV1pqBAwcyYMAA9u3bx969e0lNTSUyMtLhezIzMxk+fDg//vgj27ZtIy4ujjvuuOOa9eLi4gDYunUrDz30EB999BERERGABHIhSpKzxaycDfhmJVkrVqtXryYoKIgnnngCAH9/f2bMmEF4eDivv/663fdcvHiRrKwsatSoAUD58uVp1iz/B3/q1CkeeeQRTp8+Tdu2bfnuu+8YNmwY06ZNY9GiRaSlpdG2bVtatmwpFRNFqWC27A9nilnlvm6mdrtCArnVzp07ufnmm/MtCwkJoUGDBiQkJNh9T/Xq1enXrx8NGzake/fu9O3blyFDhuDnd/VCp3bt2nz00UdMmzaNZcuW5Xt/bGws7777bt4DLYTwdd5+dqU7vFm90F3mDeRflcDDh/9tfDmCjz76iB07drBy5UqmTZvGb7/9xrx58wzfjxC+oLDsD18Nkr7AvIG8BIJuYSIiIli0aFG+ZRcuXODIkSM0btyYjRs3Onxv69atad26NY8++ijh4eESyEWZ5evZH75KbnZade/encuXLzN//nwAsrOzGT16NEOHDqVCBfvPkk5NTWXt2rV5v2/dutXlioeBgYFkZmYWu91CmImvZ3/4KgnkVkoplixZwsKFC2nSpAlNmzYlKCgoX1bJqlWrqF+/ft6fuLg43nzzTZo1a0bbtm159dVXXe6NDx8+nBtvvJGHH37Y4CMSwvN8PfvDV0kZWxOTcyJ8kdmyVkoTR2VsDRkjV0odAi4C2UCWvR0JIcoGX87+8FVG3uy8U2t9xsDtCSFEqZKYmEjt2rUpV66cods11Ri5mR667G1yLoQoPfbv38/w4cNp1KhRXkKFkYwK5BpYoZTarJQabm8FpdRwpdQmpdSm06dPX/N6UFAQZ8+elQCGJYifPXuWoKAgbzdFCOGGXbt28eijj9K0aVPmzp1LVlYW8fHGV1Q05GanUipUa52klKoN/AaM1Fr/19H69m52ZmZmkpiYSHp6utvtKQ2CgoKoX78+gYGB3m6KKMXkxqR7HJ2/uLg4oqOjWbx4MVpr/P39efTRRxk/fvw1ZTxc4ehmp+FZK0qp14BUrfU0R+vYC+RCCM8qOJ0eLKmC9opKiWvZO3/6RDy19v/E5j9WA1CuXDmGDRvG2LFjCQsLc3ufJZa1opSqCPhprS9af+4JvOHudoUQJUum07sn9/xprUk/sp2UPxeQcWQ7R4AKFSrw7LPPMnr0aOrVq1fibTEia6UOsEQplbu9r7TWvxqwXSFECZLp9O5JOn+Zywc2kfLnN1w5Zhn3VuUqEHLzfez7/h1q1ap17ZtSD0KlcMPb4nYg11ofANoY0BYhhAfVqxpMkp2gLdPpC5eTk8PixYs588UELh2zVEb1Cw4hpH1/Kt90L9dfVyt/ENcaTv0XdsVAyi64dycEVja0TeYtmiWEKFGv9Gpmd4xcptPbl5WVxddff01MTAy7d+8GwL9SNUI63E+ltr3xKxec//zpHEhaBjtj4MpZiBgHYY+Af3nD2yaBXIgyytcfpuApGRkZfPbZZ0yZMoUDBw4A0KBBA8aPH0+Ndj2ZufZw/vPXpjYc/Bx2TQG/8tByAtS/H/z8i9hT8Zmm1ooQQpjJ5cuXmTt3LlOnTiUpKQmAJk2aMHHiRB5++OFrU4OzLsP+T2DPNKjUCCImwHU9QBn3bIUSrbUihBClxYULF3jvvfeYPn06uZMXW7duTWRkJIMGDcLfv0DP+sp52Dsb9r4DNW+FLgugZiePtlkCuRBCAOfOnWPmzJnMmjWL5ORkADp06EBUVBR9+/bN9whHAC4fg/gZll546H3QfS1U8U61UgnkQgivMMus0pMnTzJ9+nTee+89UlNTAejatSuRkZHcfffdqIJDIxcTYNebcHQRhD0KveOgYgOPt9uWBHIhhMeZ4SHNR48eZerUqcydOzevNEjPnj2JjIyka9eu177hXBzsioWTq6HJCOgbD0F2csW9QAK5EMLjHM0qfWnBVqYujy/R3nlCQgKxsbHMnz8/7zGL/fv3JzIykg4dOuRfOS8HPBaSt0Pzl6HTR4bngbtLArkQwuMKmz1aUr3znTt3EhMTw9dff01OTg5+fn4MHjyYiRMn0rp16/wr5+aA74qF9NMQMRa6Li2RHHAjSCAXQnico1mluYys+bJly5a8SoQAAQEBPP7444wfP56mTZvmXzknEw4vsARwv3LWHPCBJZoDbgQJ5EIIj7M3q7Qgd2u+rFu3jkmTJvHrr5bST+XLl8+rRNiwYcP8K2elwYFPYPdUqBgON02H6+4uMgfcLDdsJZALITzOdlapo555UTVf7AXR/m3rsWrVKqKjo1m7di0AFStWzKtEWLdu3fwbuZIM+96D+FlQszN0+cbytxNcuWFb0gFfZnYKIYrNiABVnLroBd+jtSb70CYq7vqBvf/EAVClShVGjhzJiy++SM2aNfNvIO047JkB+z+G0L6WOihVIlxqd5fY1Xa/hEKrBrNu/F1uHZ8jMrNTCGEoo1IIi1PzJa8WeE42l/euJ2X9AjJPHQSgZs2ajBo1iueff54qVarkf+PFBMvwyZGFlgJWvbdAxYZ29lA0Z8sAe6LuuwRyIUSxGBmgBrQLdek9SWcvkrr7P6SsX0jWuUQA/CtVJ6TjQA798DYVK1bM/4bzW2FnLJxcCY2fNSQH3NkywJ6o+25YIFdK+QObgCStdV+jtiuE2ZnlhpeneePBFBkZGcybN48TH79OxvnjAPiH1KZK50FUat2D+jWrXA3iWsPp3y0BPHkbNB8FneZAYEix9l3wc76zeS2+25xUZBlgT9R9N7JH/iKwGyjeWRLCB5lhhqIn2QYzP6XItnOPrSQeTHHp0qW8SoTHjh0DoFz1UCp3/hcVI7qh/AOuBlGdA8d+ttQBTz9pzQFfDP5Bxd6/vc/5u81JPHBzKGv2nC70S9wTdd8NCeRKqfrAvUA08LIR2xTCF3jjuZfeugIoGMzsBXGjA9SFCxeYPXs2M2bMyKtEeOONNxIZGYl/o85MX5mQdx7G9ryB/lXXws+x4BcAEePh+kGG5IA7+pzX7Dmd78amPZ6o+25Uj/xtYCzgcN6qUmo4MBwsRdmFKA08PbzgzSsAe8EMwN/aM/dXKu9LzN32nD17lpkzZ/LOO+/kVSLs2LFjXiXC3EJWD7RvYM0B/xR2D4NzDaHdNKjb09A64O5+zq7eA3CV24FcKdUXOKW13qyUusPRelrrOcAcsKQfurtfIczA08+99MYVQC5HQStba4ID/Q35cjlx4kReJcJLly4B0K1bN6KioujevXv+SoRXkmHf+xA/E2p0glu/glq3AMZftZj9+aZ+Ra9SpC5AP6XUIeAb4C6l1BcGbFcI03ulVzOCA/Nfupfkcy+9cYMxl6OgldsTt2XbM3fGkSNHGDlyJOHh4UydOpVLly5xzz338Pvvv7N27Vp69OhxNYinnYCt4+GHGyBlN9y1Erp9ny+IT1i8g6TkNDRXv1iWxiUV67jB85+zq9wO5FrrCVrr+lrrMGAwsFpr/YjbLRPCBwxoF0rMwNaEVg1GYZkMUpyJHs5yFEw90TN0FMzsjZWDc18u+/btY9iwYdxwww28++67pKenM2DAAP7++29++eUXbrvttqsrX9wPG0fATxGQdQnu2Qy3zoeqrfJts7CrluLy9OfsKskjF8JNJT3+acsTGRCOOLpp52iavaMvl6VxSbz22a8krPicS3t+B22pRDhkyBAmTpxIq1b5AzPnt1mKWJ34zZoDvgeCajtsZ0ldtXjyc3aVoYFca70WWGvkNoUQV3kiA6Ko/dvbl7NfLtO+/JnX/28SqfHrLQv8/KnSpheTX4vkuf6351/51O+WFMLkrdBsFHT80KkccLOPZ5cE6ZEL4WPM1jN05svljz/+YNKkSSxfvtyywD+Qym16EdJpIAEhtflydybP9ccyiefYT5YeeNqJYuWAe/OqxVskkAsh3Gbvy2XJlkSiZn/F/t8+J+PoPwCowCAqt+tDSIf78a9ULW/dk8mpcPBL2D0FlL9bOeDFuWrx9dm5EsiFEIbSWhM5cx4zpk4h/ZjlBqNf+YpU69if6265n1R1dYijvMrgweorea7OEtjfDNq+CXV7uZ0D7spVS2mYnSuBXJR5vt4b8zRH5ys7O5tFixYxefJktm/fDoBfcAghHe+ncrt78StfgcAKgQRn5hCYncIjNX7miZo/sD29Oftu+JB6Xfp55Xi8mZtvFAnkokwrDb0xT7J3vsYvjGP1jwtZ/uX77N27F7BWIuz0AJXa9MIv8Or4dsCVU3x/+3rqnPqcVRfaM+rMNAbddY9Xhz28mZtvFAnkokwrDb0xT7I9XzrrCqk7VpL413fEp5wEICwsjPHjxzPvdEOOp149r9eXO8Eztb6jX7XfCan+OHTZzsBKYQwsYn+FfdHmtsfdAO9OlotZruYkkIsyrTT0xjzpWHIaOVfSSd32Kxc2LiY79RwAAdXr8/GMaIYMGUJgYCB1rAG4oX8CI2ot4vbKcXyb3Js/mv1Bnw5tnN6foy/a137YSUZWjiFXUsXNcjHT1ZwEclGmlcWc4+JKSUkhJ24xSb8vIiftAgCBtcKocutgmnS8i8ceuztv3QGhB+ncYRoByVv56PR9vJs6mud73kQfFwOcoy/U5LTMa5a581ALcL13b6arOQnkokwriznHrjpz5kxeJcKUlBQAytVtRpVbHyL4hg5UKBfA2N4R1hzwn2FXDKQd57oWY6HRMsb7BzG+mPt29EXrSHGvpIqTm2+mqzkJ5KJM8/ZMSTM7fvw4b731Fu+//z6XL18G4I477uDOh57l1/M1OZ6SfrUOeLX/wC+xgJ8lB7zBIEtNcDc5+qINCvTj/OVre+WevJIy09WcBHJR5nlypqRZbo4V1o7Dhw/z5ptv8vHHH5ORkQFA7969iYyMpEuXLgD8P4DsdEsd8F1Pwdn60HYK1L3HYQ54cY7d0RctOF8WoKSY6WpOArkQHmKWm2OO2pF0+ACbf/iUzz//nKysLAAGDhzIxIkTufnmm69u4ErK1Trg1dvDrV9ArVuLtU9wfOwFA/+Mh9pes643vxTNdDWntIMSlCWpffv2etOmTR7frxBGcrWH2SV2td1L8dCqwUU+LsxIBdtx5fQhUtZ/y+U9f+RVIhw8eDATJ06kZcuWV9+YdhLi34aEOVCvN0SMg6qti7XPXI6OvWDgB0tv10ylY71BKbVZa92+4HLpkXuQWS6rhfuK08M0y82x3P1lHN9LyvpvSdu3wfKCXwBPDXuScePG0bhx46tvSD0Iu6fCoa8h7N9wzyaoFF6sfTq73EwZIb5AArmHmOWyWhijOIHGLDfHKp3fR8Jvn5N+cAsAKqAcldr0ommPfzM3dvDVFZN3wM5YOP4rNH7GUgc8uE6x9unqsZvlS89XGPGoN+GEknhqifCe4gQabz4uTGvNihUr6NatG//MGUX6wS2ocsGEdHqA0Gc/JrT3c0Q9ZK0HfnodrO0Lq3tCtRuh3wFoO7nYQRxcP3ZvPgnJFxnx8OUg4L9Aeev2FmmtX3V3u6WN9DBKl+L0rr1xcywnJ4cff/yR6Oho/v77bwCqVq1Kr389wYE6XTl9JdDSjp5NGVB7G/w2GC4nQcQrcNtCCDAmcLp67GbKCPEFRgytZAB3aa1TlVKBwB9KqV+01hsM2HapYZbLanGt4ty7KG6g8VSqY3Z2NgsXLmTy5Mns2GEZwqtVqxajR49mxIgRhIRYn7STkwVHFsKuUXACaw74g4bkgBfkyrGbKSPEF7j9aWlL2kuq9ddA6x/Pp8KYnPQwzKm49y7MGmgyMzP54osviI2NzatEGBoaytixY3nqqaeoUKGCZcXsdDgwz3ITMzgU2sRYMlHcrANuJLM9CcnMDEk/VEr5A5uBxsBsrfU4O+sMB4YDNGjQ4ObDhw+7vV9fI1kr5mOWlEB3paen8+mnnzJlyhRy/2+Fh4czYcIEHnvsMcqXL29ZMfOCJQd8z9tQ/WZoOQFqdfFew4VLSjT9UGudDbRVSlUFliilWmmt/ymwzhxgDljyyI3Yr6+RHob5+Pq9i9TUVD788EPeeustjh8/DkDz5s2JjIxk8ODBBARY/4unnbRM4Nk/xzL78q4VTueAC/MzdCBMa52slFoD3AP8U9T6Qnibp+5dGH01lpyczOzZs5kxYwZnz54FoE2bNkRFRTFw4ED8/KwJaakHYfc0OPw1NBwCvTZCpUamOQ5hDCOyVmoBmdYgHgzcDUxxu2VCeIAn7l0YOYfgzJkzvP3227zzzjtcuGApJdu5c2eioqLo06cPKneMO3kH7JoCx36BxsPh3t1upQ8afRzCWG6PkSulbgQ+A/yx5KV/q7V+o7D3yBR9YSbF7WU6+z4jxuGPHTvGtGnT+PDDD/MqEd55551ERUVx5513Xg3gp/+EnTFwbhM0exGajIByVZzaR1EcHYe/UuRoLT10DyixMXKt9XagnbvbEcJbinPvwpXeqaPxdmfqbB86dIgpU6bwySefcOXKFQD69OlDZGQkt95qLVSltaXnvTMGLidac8C/NSwHPJej48i2dgalh+49MrNTlApL45LoErua8PE/0SV2NUvjkkp0f67M1HU03q7AYTvj4+MZOnQojRs35oMPPiAzM5NBgwaxZcsWfvrpJ0sQz8mCQ9/AL+1g6zjLNPr79lp64QYH8cKOw5bMVvYOCeTC5+X2jpOS09Bc7RmWZDB3JdvllV7NsJedreGaoLd9+3YGDx5MixYt+OyzzwB49NFH2blzJwsXLuQwtbkz9lcmvvkCiZ835GzcdGgTDb23QfjDJTKRx/Y4Ck6zt8dXMn5KEwnkwud5o46NK7VABrQLdThDLjfobdy4kf79+9OmTRsWLFhAQEAATz/9NPHx8cyfP58WLVqwbNMeEtb+P76p+zDdQzby0uEXuW3LGyw91dYjE3kGtAslZmBrQqsGo7CMjdsjs5U9T6ofCp/njVxwV7JdlsYl4a9U3liyrYrn99GzZ09+++03AIKCghg+fDivvPIK9evXt6yUfgriZ3Hb7nfILteOxw++zp703DKyni3tans/wVHNcJmt7HkSyIXP80YdG2en6OcGO9sgrrUm/eAWLm74lsNHd7ITqFSpEs8//zyjRo2iTh1rmmDqIWsO+FfQcDD99k3nyJXrrmmLt4YyzFqmoCySQC58nrfq2DiT7WI77KN1Dmn7/iJl/QKunEgALJUIX3rpJUaOHEn16tUtb0r+x5oD/jM0fhru3QXB15G9ajVcMVfhNZmtbA4SyIXPM3PP8FhyGjonm8t7/iBl/QIyzxwBwK9CVWJencCIESOoXLmyZeXT62FXDJzdaMkBb/8OlKuaty0pvCYckUAuSgV7PUNvTye/cuUK/glrSVz9JVnnLXVQ/CvXJKTTAzS9vR9jx/ax5oD/agngl45Ai1egywK76YNm/sIymrc/O18jgVyUSt6cTp6WlsYnn3zClClTOHr0KAABVa8jpPODVGp1FxWCghjfNwIOL4BdsZZ88Ijx0PChItMHy8JQhpQCcJ0EclEqGfHwXld7hampqXzwwQdMmzaNkydPAtCiRQt6PTyCv2jO8YtXaFjVnxkdt9Eu8QU4Uwdu/D+od6+p6oB7mzx42XUSyEWp5G5Koiu9wuTkZN555x3efvttzp07B0C7du2IiopiwIABlkqEmRdg34cQPwP82kHnT6H27cU9PI/xxhCHr5cW9gYJ5KJUcjcl0Zle4enTp5kxYwazZ8/Oq0R4yy23EBUVRe/evS2FrKw54CR8ANf1hDt+gWpt3Dw6z/DWEIc8FtF1MrNTmFpxa6i4+8T6wnqFx44d4+WXXyYsLIyYmBguXLhA9+7dWb16NevWrbOUk710GDaNhGXNIeMM9PwLunzlM0EcvDNjFtz/7Moi6ZEL03KnR+huhoe9XmFWykmy4pYSPmN5XiXCvn37EhkZSefOnS0rJe+05oD/BDc8BffuhOC6Th+zmXhriKMsZecYRQK5MC13b3q5k+Fhm7OdeTaRlA0LubRzDegclFIMGjSIiRMn0q6dtYLz6fWWDJSzf0Gz/4H2s/LlgPsibw5xlIXsHCO5PbSilLpeKbVGKbVLKbVTKfWiEQ0TorA63iVdpnZAu1CGt/In9eepHPtoBJf+WYWfn+Kxxx7Lq0TYrm1bSw74yjvgzyFQtyf0OwAtJ/p8EAcZ4vAlRvTIs4DRWustSqnKwGal1G9a610GbFuUYY56hIBLN91czbzYsGED0dHRLFu2DIBy5crxxBNPMHbsWBo1agQ52XD4W2sO+BWbHPDAYhylZ7lyLmSIw3cY8YSg48Bx688XlVK7gVBAArlwi70p6blsb7oVFmicHWfXWvOf//yHSZMmsWrVKgCCg4N55plnGDNmDKGhoZCdAQlzYdebEFQbbnwD6vUB5Rs5A8W55yBDHL7B7Wd25tuYUmHAf4FWWusLBV4bDgwHaNCgwc2HDx82bL+i9Foal8RLC7Y6fF1BvlrfwYH+xAxsnRd8inpeptaaX3/9lejoaNatWwdA5cqVeeGFF3jppZeoXbs2ZF6EhA9hzwyo2gZajodat/vcJB4jnh0qvMvRMzsN60oopSoB3wEvFQziAFrrOVrr9lrr9rVq1TJqt6KUG9AulNBCbq4V7IYUTI9zOM5+/hJLliyhQ4cO9OnTh3Xr1lG9enVef/11Dh8+zOTJk6kdomBbFPwQDuc2wx0/wZ0/Q+2uPhfEQSbalGaGZK0opQKxBPEvtdaLjdimELkKG2KxxzYwFRxnt1Qi/J1LGxcx8M1DANSuXZsxY8bw7LPPWioRXjoMm16FQ19Ag39Bzw1QubGhx+QNMtGm9HI7kCulFPAxsFtrPd39JglP8KXqcrY33Zx58rxtYMr9Ericnk7qP6u5sGERWcmWSoT169dn3LhxDBs2jODgYEsO+Po3IWkZ3DDMp3PA7ZEyuKWXET3yLsCjwA6l1Fbrsola658N2LYoAb5YXS73ppujcd5cBQNTr+bV+fnKJj6b+w5XUk4BcF39hkx67X959NFHKVeuHJzZAH/Hwpn10PR/oF8ClKtW4sdkhKilO/j6r6Nka42/UgzpdD2TBrS2u65koZReht7sdFb79u31pk2bPL5fYeHsTS8z9trtPScy94ZnqE0bL168yAcffMBbb72VV4kwIiKCyMhI/vWvfxHg7w8nfoOdMXDpIDQfAzc8CQEVvHNgDhT2GUQt3cEXG45c855HOjdwGMyFb3N0s1NmdpZBztz0Mmuvvahe5fnz53n99deZOXMm58+fB+Cmm24iKiqK/v3744eGxMWwMxZy0q054IO9ngNuL2ADhX4GX/911O62vv7rqATyMkYCeRnkzE0vM9eEtpfbfOrUqbxKhBcvXgSgS5cuREVF0atXL1TOFTjwCex+E8rVgNavQmhfU+SAO/rSLB/gV+hnkO3gatrRclF6SSAvg5y56eUrqWqJiYlMmzaNOXPmkJZmaVuPHj2Iioqia9euqKxU2DPdmgPeCjrONV36oKMvTUdZOrmfgb9SdoO2v4mOTXiGBPIyyJmbXmZPVTtw4ABTpkxh3rx5eZUIq7W4haD2g7gc0ZZLFaqhdrwK+96HOt3hjmVQra13G+2Aq1+OuZ/BkE7X2x0jH9LpekPaJXyHBPIyqqip12ZNVduzZw8xMTF8+eWXZGdno5Siy933caRBT/xqhlMv8BRPVXibW3ev4WDN/oT3XG/6HHBHX5rVKgSSnpnj8DPIHQd3NmtFlF6StSIcMlPWytatW5k8eTKLFi1Ca42/vz8PP/wwEyZMYMi3R6iRvZ9na31Hj5C/WHC+Jx+f7k9mueuI+389vdJeV9jLxMktNQCSLiiukqwVcY2iArUZCiatX7+e6OhofvrpJ8BSifDJJ59k7NixhIeHw5mNxNZ6lZsq7OGzs33pGv8RF7IrWd6clenFljuvqKEub38GwvwkkJdRZk0vBEslwrVr1xIdHZ1XibBChQo888wzjB49mtB69eDESlg1DC7u58/U3rx4ZAzpOsir7XaHGb40he+SQO4CMw01uMuM6YVaa3755Reio6P5888/AUslwpEjR/LSSy9Rq0Z1SFwCy2Mh67IlBzxsCN//3xrS9bW9bwWEj//J5z8rIYoigdxJZu7BFoeZ0gtzcnJYsmQJ0dHRxMXFAVC9enVGjRrFCy+8QNXKwZYCVhvetEydb/W/EHpfXg74a/1a8srCbWTm5L/fk/ubr39WQhRFArmTzNiDdYcZ0guzsrL45ptvmDx5Mrt37wagTp06eZUIK5UHEubAmunWHPAPoXa3a3LAC44x+9nJry7ssypNV1qibJJA7iQz9WCN4M30woyMDObPn09sbCwHDhwA4Prrr2fcuHE8+eSTBKtLsHcq7HsP6twF3X6E6u0K3abtGHP4+J/srmPvsyptV1qibPL+/GQf4ainapYJMq4a0C6UmIGtCa0ajMJScMr2yTol4fLly8yaNYvGjRszfPhwDhw4QOPGjfn4449JSEjg+aH3EbxrAixrCmnH4e4/4bYFRQbxglz5rAq70hLCV0iP3ElmnSDjDk9lSly8eJH333+ft956i1OnLKVkW7VqxcSJE3nwwQcJuLQPtjwDiT9YKhD2+Qcq1Cv2/lz5rErblZYomySQO0lqObvu3LlzzJo1i1mzZuVVIrz55puJioqiX79++J3bBH/+C86sgyYvGFYH3JXPytG9Ao2l3K98xsIXyMxOYbiTJ08yffp03nvvPVJTUwG47bbbiIqKoufdd6NOroJdsXAxAVrk1gGv6JW22ptVaavgw5yF8KYSndmplPoE6Auc0lq3MmKbwvckJiYydepU5syZQ3p6OgA9e/YkMjKSrrd1gcSlsKIjZKdBi3EQNsTrdcCLeoycL2cmibLDqKGVecC7wHyDtic8yN30u/379+dVIszMtEzM6devH1FRUXS4qQ0c+hx+etpuDrgZ5N4rCB//E/auT2W8XJidIYFca/1fpVSYEdsSnuVO+t2uXbuIiYnhq6++IicnB6UUDz30EBMnTuTGFo0sOeA/3A9VWjrMATcTM+TWC1EcHusWKaWGK6U2KaU2nT592lO7FUUoTvpdXFwcgwYNolWrVnzxxRcopRg6dCi7d+/mm3nvciPfwQ/hcHYDdPsB7loOde4wdRAHS7ZLoF/+Ngb6KZ/OTBJlg8eyVrTWc4A5YLnZ6an9isK5kn63fv16Jk2axM8//wxYKhEOGzaMsWPHElbLH/a8Bcvmw/UD4e51ENK0RNteIgp+17j43SOzRIU3SPphGVfUcILWmjVr1jBp0iTWrFkDWCoRPvvss4wePZp6FS/A7jdg01Jo9CT02QEVfDNwTV0eT2Z2/j5GZrZ2+manK8NUEvCFkcxzx0l4xSu9mhEc6J9vWXCgP2N6NmXZsmXceuutdO/enTVr1hASEkJkZCSHDx/mrYmDqbd/JKzsChXD4b4EuGmazwZxcH9ykLPDVLkBPyk5Dc3VgL80LqlY7RbCqPTDr4E7gJpKqUTgVa31x0ZsW5SsgpNn6oaU57bA/bwxbAJbt24FoEaNGowaNYrnn3uOqhmbYdtguLgXmo+GW+Z7LQfcaO7e7HT2i6C0FWAT3mdU1soQI7YjvGNAu1D6tq7D119/zeTJk5m6Zw8AdevWZcyYMTwz/Gkqnv8N/uoJWZcgYhw0HAL+5bzccmO5W4bB2S8CKQsgjCZj5GVcRkYGn332GbGxsRw8eBCAhg0bMm7cOJ54/GGCjn8HaztAYBVoGQn1+5kqB9xI7pZhcPaLQNIchdEkkHuBGW50Xb58mblz5zJ16lSSkixjs02bNmXChAk8/FB/Ag9/BitaQkgLaP8e1LnT9OmDRnCnkJizXwSlsQCb8C6pteJh9mp7KCxFmkI9ENQvXLiQV4kwN5+/devWREZGMui+O/Hf/z7snW2ZvNNyPFS/ucTaUpaZ4ctc+J4SrbUinGfvRldxHknmaiDIrUQ4c+ZMkpOTAejQoQNRUVH0vastfvFvw88jrDngv0OI9A5LkjxsWRipdA52mlhRN7SceaiBK+lrJ0+eZNy4cTRs2JDXX3+d5ORkbr/9dlasWMFfv82nX+3v8fu1rWXYpM926PSRBHEhfIz0yF1gxOWwoxtdtooK9s6krx09epSpU6cyd+7cvEqEvXr1IjIyktsjgi1lZFf+F5q+APftg/I1XDoOIYR5SCB3klHPdrR3o6ugorIXCktfS0hIYMqUKXz22Wd5lQgHDBjAxAkT6NDgIux6A36Pt+aAf1ZqcsCFKMskkDvJqEkcBetf597ozFVU9sLSuCS7T4m/cvowmZsX02zqGnJycvDz82Pw4MFMnDCe1tX2w84X4NRFaw74v+3mgLt7xVGc95eGm36uHENpOF5hPhLInWTkJA7bG12uBoEJi3fkC+IZJxJIWb+AtL3rAQgICGDo0KGMH/syTQI2wu6H4FgItJwA9fs7zAF394qjOO8vDU+wd7W+iq8frzAnCeROKqlJHK5kL9heFaQn7iLlzwWkH9wMQGC58gx/+inGvvw8Da6sgN29LTctncwBd/eKozjvLw1T1V05htJwvMKcJJA7yRuTOAr21hPPXyb98DZS1i8g44ilJ6cCy1O5bR/iF/0f16UsgrhuULsr3L4YalyTbmp3u6/0aubwyiIpOY2lcUlFBpriXLGUhqnqrhxDaTheYU4SyJ3k7vRtV9lehmut2bdpLSl/fsuV45bURFWuAiE330eTW7swuvF/uG7j7VD//iJzwB1d3letEMj5y5l23+PM5X9xrlhKw1R1V46hNByvMCcJ5C7w5CSOqcvjuZxxhcvxf5Ky4VsyT1nqoPgFhxDSYQBtbm3HiNCf6VVlImdqDYGu26FCfae2a+/yvnyAH8GB/nazadIysxn97TZGLdhq6LTz0jBV3ZVjKA3HK8xJArkJZWZmsvf3ZSRvWEjWuUQA/CtVJ6TjQG65pQljwpbRutz3LE0dwLrm67m3Qyunt+3oMj4lLZMZD7XlpQVb7b6ee4PV0Q264lyxePoqpyS4cgyl4XiFOUmtFRPJyMjg008/ZcqUKRw6dAgA/yp1qNJpID1vuY7n6i6haXASdTpOhMZPFysHvEvsaruX96FVg1k3/i6HrztaXwjhOY5qrRgyRV8pdY9SKl4plaCUGm/ENsuSS5cu8fbbb9OoUSNGjBjBoUOHCA27gXr3jeLJMU+w6sFVTGowl18u3smGVn9B85eKPZHH0ROBci/v7b1uj9ygE8I83B5aUUr5A7OBu4FE4G+l1A9a613ubru0u3DhArNnz2b69OmcOXMGgBtvvJH/nTiOgTenc2lrDIkX/Zl5fBD/+N3JmF4t6O/mZXhRl/cFX7c3+QjkBp0QZuL20IpS6hbgNa11L+vvEwC01jGO3lMWhlYKm+hz9uxZZs6cyTvvvJNXibBjx468GjmG3k2TUHumWzJPIsZDnbu8WgfcXtnd4EB/Yga2lrFdITysJMvYhgJHbX5PBDrZacBwYDhAgwYNDNiteTlK8Tt/5hS7VnzF+++/z6VLlwDo1q0br038H7rV3YHa+zycvh1u/w5qdPBq+22/hB64OZQ1e07LDTohTMpjWSta6znAHLD0yD21X28omOKXdeEUSX99x7DYFegsS652r169eGP8M3QM+QMOPAUhA6DHf6FKcy+12sLel9B3m5OkBy6EiRkRyJOA621+r29dVmblZn1knksiZcMiLu1cDTmWwHj//ffzxiuP0MrvZzg6DCo/Dr23QcXrC9ukx8g0ciF8jxGB/G+giVIqHEsAHwz824Dt+qSlcUlknj5E8vqFXN7zO+gcUH5UiOjGff278k3v3XDsGWjyHPTdC0E1vd3kfMr6NHKpTih8kduBXGudpZR6AVgO+AOfaK13ut2yYvD2f8JNmzYx7PGXOLdrnWWBnz8VW/Wgd+8bGdloDR2qzYVaY6HzpxBYyWPtckVZnkYu1QmFrzIkj1xr/bPWuqnW+gatdbQR23SVK48/M9off/zBPffcQ4cOHSxB3D+QkJvu5YkJz/H7yKPENl/Aj8m3E/zAIWg+yrRBHIrOMy/NChtWEsLMSs0UfU+P7WqtWblyJdHR0fznP/8BoGLFitTo0JuBd1Tnuet/IS0nnvdOPciKC52pW7US+Jc3vB1GK8vTyMv6sJLwXaUmkHvqP6HWmmXLljFp0iQ2btwIQJUqVRj94rO83L8KKmE221Jq8dqxZ1iX2gZQPtejLatPeC/Lw0rCt5WaQF7S/wmzs7P57rvviI6OZvv27QDUrFmTiWOeZUT3HIKOfASXboO7FnPyaCiHlsejMEeP1tv3DnyFVCcUvqpUBPKlcUlcvpJ1zXIj/hNmZmby1VdfERMTQ3y8Zay0Xr16vDbuaYZ2Ok/g0dmQ1R+6r4UqLQAYUNM8N8fkBp7zyvKwkvBtPh/I7U0hB6gaHMhr/VoW+z9heno68+bNy1eJMCwsjJgJT/Bgq0P4H5sFAY+ZKgfcHskLd01ZHVYSvs3nA7m9QAVQsXxAsf5DXrp0iQ8//JBp06Zx/PhxAJo3b87UCUPoE7Ydv9PvQCVz5oDbIzfwhCj9fD6QGxWoUlJSmD17NjNmzMirRNimzY28PeF+utX6E5UyB2q9DLfMc5g+aMaxaLmBJ0Tp5/OB3N1AdebMmbxKhCkpKQB07tSRWePupn2FVagrX0PDsRD2SKHpg2Ydi5YbeEKUfj4fyIsbqI4fP85bb73FBx98kFeJsMdd3Xh7VCciWIYK+BVajLc80NjP8YMWcnvh9r5MvD0Wndu2tMxs/K11xUNNcqUghDCOzwdyVzMNDh8+zJtvvsnHH39MRkaGZRt9e/LWcxE0ylgMwf4Q8TZc16PIOuCObrTa8tZYdMG2ZWud9wUnQVyI0sWnA3nBMekZD7V1GKT27t1LbGwsn3/+OVlZllTFR/51L7HDric0dTFUrAgdvoWa15RSd8jRjVZb3hqLlmwVIcoOnw3kzo5J79ixg8mTJ/Ptt9+Sk5ODn58fzz95P68NrkLNlO+hcj/otDYvB9wVRfW2vTkWLdkqQpQdPhvIi+px/v3330RHR/P9998DEBgYyLjnH2DcfX5USVkB1R6FW7dCxeI/rcjRjVbA62PRkq0iRNnhs4HcUc/ywI6/6dVrCitWrAAgKCiI1/6nPyO7X6LCxTVQZwTc/q4hOeCObrSa4Wk6kq0iRNnhs4HctseptSb9UBwpfy4gI3EnJ4BKlSoy9ZX7eKLjScqn/Q7XvwyNv4LAyoa1wcxTus3cNiGEsZTWnn98Zvv27fWmTZvc2sbSuCTGf7eNc7vXk/LnAq6c2AdApZAQ5kT1YVDL/QRmJ0OLVyD8MZ8oISuEEIVRSm3WWrcvuNytHrlS6kHgNaAF0FFr7V50dlJ2djbp8b+TtuA1TidYClkFhVRhxpjbebrtfvwD4qHlBKg/sNAccCGEKA3cHVr5BxgIfGhAW4qUmZnJF198QUxMDPv2WXrgN4TV5aMJXeha+2/8Kl+y5oDfXWQOuBBClBZuBXKt9W4A5YGgefnyZVq1asXBgwcBuLFFAz4a2472lTegamVBxDdQs3OJt0MIIczGYzc7lVLDgeEADRq4nvJXoUIFOnXqRP0a/rz/UlMiym9AhVaFiNVQJcLg1gohhO8oMpArpVYC19l5KVJr/b2zO9JazwHmgOVmp9MttDH3jQepuG05KqwxtHgPKjYszmaEEKJUKTKQa617eKIhzqgU1hOuj4egWt5uihBCmIZv5ZEHVnJYC1wIIcoqP3ferJS6XymVCNwC/KSUWm5Ms4QQQjjL3ayVJcASg9oihBCiGHxraKUYzPj4NSGEMFKpCOSOgrVZH78mhBBG8vlAXliwlocrCCHKAp8P5IUF6+I8XEGGYoQQvsatrBUzKCxYO3qIgqPlub37pOQ0NFd790vjkoxqrhBCGM7nA3lhwfqVXs0IDsxf/bCwhysU1rsXQgiz8vlAXliwHtAulJiBrQmtGozC8vi1wp7e4wvPuVwal0SX2NWEj/+JLrGr5WpBCOH7Y+RFPQlnQLvQfIE7NxDaW9fsz7mULBwhhD0+H8jh2mDtSFGB0OzPuZQsHCGEPT4/tOKKosbAXR2K8TRfGPoRQnheqeiRO8uZQOhs7x48n6po9qEfIYR3lKkeuavpiIXxRqqiq1k4QoiyoUwFciMDoTdSFc0+9COE8I4yNbRSVIaLK7w1Xu3K0I8QomwoU4EcjAuEMl4thDALdx8sMVUptUcptV0ptUQpVdWgdhmqJCbRyHi1EMIs3B0j/w1opbW+EdgLTHC/ScYqqZuSMl4thDALd58QtMLm1w3AIPeaY7ySnEQj49VCCDMwMmvlSeAXRy8qpYYrpTYppTadPn3awN0WTibRCCFKuyIDuVJqpVLqHzt/+tusEwlkAV862o7Weo7Wur3Wun2tWrWMab0TjMwdF0IIMypyaEVr3aOw15VSQ4G+QHettTaoXYYxe/0UIYRwl1tj5Eqpe4CxQDet9WVjmmQsI3PHhRDCjJQ7nWilVAJQHjhrXbRBa/1sUe9r37693rRpU7H3K4QQZZFSarPWun3B5e5mrTR25/1CCCHcV6ZqrQghRGkkgVwIIXycBHIhhPBxEsiFEMLHuZW1UuydKnUaOFzMt9cEzhjYHKNIu1wj7XKNtMs1Zm0XuNe2hlrra2ZUeiWQu0Mptcle+o23SbtcI+1yjbTLNWZtF5RM22RoRQghfJwEciGE8HG+GMjneLsBDki7XCPtco20yzVmbReUQNt8boxcCCFEfr7YIxdCCGFDArkQQvg4UwZypdSDSqmdSqkcpZTDNB2l1D1KqXilVIJSarzN8nCl1F/W5QuUUuUMald1pdRvSql91r+r2VnnTqXUVps/6UqpAdbX5imlDtq81tZT7bKul22z7x9slnvzfLVVSq23ft7blVIP2bxm6Ply9O/F5vXy1uNPsJ6PMJvXJliXxyulernTjmK062Wl1C7r+VmllGpo85rdz9RD7RqqlDpts/+nbF573Pq571NKPe7hds2wadNepVSyzWsleb4+UUqdUkr94+B1pZSaZW33dqXUTTavuXe+tNam+wO0AJoBa4H2DtbxB/YDjYBywDYgwvrat8Bg688fACMMatebwHjrz+OBKUWsXx04B1Sw/j4PGFQC58updgGpDpZ77XwBTYEm1p/rAceBqkafr8L+vdis8xzwgfXnwcAC688R1vXLA+HW7fh7sF132vwbGpHbrsI+Uw+1ayjwrp33VgcOWP+uZv25mqfaVWD9kcAnJX2+rNvuCtwE/OPg9T5YHoepgM7AX0adL1P2yLXWu7XW8UWs1hFI0Fof0FpfAb4B+iulFHAXsMi63mfAAIOa1t+6PWe3Owj4RZf8QzdcbVceb58vrfVerfU+68/HgFNASTwL0O6/l0Lauwjobj0//YFvtNYZWuuDQIJ1ex5pl9Z6jc2/oQ1AfYP27Va7CtEL+E1rfU5rfR74DbjHS+0aAnxt0L4LpbX+L5aOmyP9gfnaYgNQVSlVFwPOlykDuZNCgaM2vydal9UAkrXWWQWWG6GO1vq49ecTQJ0i1h/Mtf+Ioq2XVTOUUuU93K4gZXkA9obc4R5MdL6UUh2x9LL22yw26nw5+vdidx3r+UjBcn6ceW9JtsvWMPI/5NzeZ+rJdj1g/XwWKaWud/G9JdkurENQ4cBqm8Uldb6c4ajtbp8vtx4s4Q6l1ErgOjsvRWqtv/d0e3IV1i7bX7TWWinlMHfT+k3bGlhus3gCloBWDksu6TjgDQ+2q6HWOkkp1QhYrZTagSVYFZvB5+tz4HGtdY51cbHPV2mklHoEaA90s1l8zWeqtd5vfwuG+xH4WmudoZR6BsvVzF0e2rczBgOLtNbZNsu8eb5KjNcCuS7ioc5OSAKut/m9vnXZWSyXLAHWXlXucrfbpZQ6qZSqq7U+bg08pwrZ1L+AJVrrTJtt5/ZOM5RSnwJjPNkurXWS9e8DSqm1QDvgO7x8vpRSIcBPWL7EN9hsu9jnyw5H/17srZOolAoAqmD59+TMe0uyXSilemD5cuymtc7IXe7gMzUiMBXZLq31WZtfP8JyTyT3vXcUeO9aA9rkVLtsDAaet11QgufLGY7a7vb58uWhlb+BJsqScVEOy4f2g7bcPViDZXwa4HHAqB7+D9btObPda8bmrMEsd1x6AGD37nZJtEspVS13aEIpVRPoAuzy9vmyfnZLsIwdLirwmpHny+6/l0LaOwhYbT0/PwCDlSWrJRxoAmx0oy0utUsp1Q74EOintT5ls9zuZ+rBdtW1+bUfsNv683Kgp7V91YCe5L8yLdF2WdvWHMuNw/U2y0ryfDnjB+Axa/ZKZyDF2llx/3yV1B1cd/4A92MZJ8oATgLLrcvrAT/brNcH2IvlGzXSZnkjLP/REoCFQHmD2lUDWAXsA1YC1a3L2wMf2awXhuVb1q/A+1cDO7AEpC+ASp5qF3Crdd/brH8PM8P5Ah4BMoGtNn/alsT5svfvBctQTT/rz0HW40+wno9GNu+NtL4vHuht8L/3otq10vr/IPf8/FDUZ+qhdsUAO637XwM0t3nvk9bzmAA84cl2WX9/DYgt8L6SPl9fY8m6ysQSv4YBzwLPWl9XwGxru3dgk5Hn7vmSKfpCCOHjfHloRQghBBLIhRDC50kgF0IIHyeBXAghfJwEciGE8HESyIUQwsdJIBdCCB/3/wHk2UcKoifSOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x[:, 0], y)\n",
    "\n",
    "xs = np.linspace(-1, 1)\n",
    "xs = np.stack([xs, np.ones_like(xs)], axis=1)\n",
    "ys_true = xs@beta_true\n",
    "ys_fit = xs@beta_ols\n",
    "plt.plot(xs[:, 0], ys_true, color='k', lw=2, label=\"True\")\n",
    "plt.plot(xs[:, 0], ys_fit, color='orange', lw=1, label='OLS fit')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "380ba33a",
   "metadata": {},
   "source": [
    "# TBC..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbd1bd0d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}