exponential-window-nutpie-mass-matrix.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1ab59106-67a7-42af-b324-9b48d17f21a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from scipy import stats, linalg\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6c94f3e2-6a62-47b0-be1f-b14f39dfb8f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b983b441-dfb6-4a11-8383-ca4c98a4ba72",
   "metadata": {},
   "source": [
    "## Generate a covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3dca9dbf-e97e-4fd2-ac55-5bb4530202f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "\n",
    "# Lognormal eigenvalues, but one very small and one very big\n",
    "sigma = np.exp(rng.normal(0, 1, size=N))\n",
    "sigma = np.sort(sigma)\n",
    "sigma[-1] *= 5\n",
    "sigma[0] /= 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "92293f1c-4a32-4df2-b85b-79d1b0d3da5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate matrix of eigenvectors\n",
    "scale = 0.05 # Small scale -> more similar to identity\n",
    "A = rng.normal(0, scale, size=(N, N))\n",
    "# Make skew-symmetric\n",
    "A = (A - A.T) / 2\n",
    "Q = linalg.expm(A) # map to ortho group\n",
    "np.testing.assert_allclose(Q @ Q.T, np.eye(N), atol=1e-14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "af2ae0c9-7393-4c17-9c91-3a8ece3ec015",
   "metadata": {},
   "outputs": [],
   "source": [
    "Sigma = Q.T @ np.diag(sigma ** 2) @ Q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "641e0729-8561-4b83-9e0d-f671f489d5cd",
   "metadata": {},
   "source": [
    "## Compute the convergence targets for stan and nutpie adaptation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4c063342-2eb4-4c64-a58e-e194cdbfd1b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The mass matrix that stan adaptation converges to\n",
    "stan_target = np.diag(Sigma)\n",
    "\n",
    "# Covariance in the posterior after (fully converged) stan adaptation\n",
    "Sigma_stan = Sigma.copy()\n",
    "stds = np.sqrt(stan_target)\n",
    "Sigma_stan /= stds[None, :]\n",
    "Sigma_stan /= stds[:, None]\n",
    "\n",
    "# The remaining sqrt eigenvalues \n",
    "sigma_stan = np.sqrt(linalg.eigvalsh(Sigma_stan))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "83258dcd-def4-4f82-84b6-2b43f08084d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The mass matrix that nutpie adaptation converges to\n",
    "nutpie_target = np.sqrt(np.diag(Sigma) / np.diag(linalg.inv(Sigma)))\n",
    "\n",
    "# Covariance in the posterior after (fully converged) nutpie adaptation\n",
    "Sigma_nutpie = Sigma.copy()\n",
    "scale = np.sqrt(nutpie_target)\n",
    "Sigma_nutpie /= scale[None, :]\n",
    "Sigma_nutpie /= scale[:, None]\n",
    "\n",
    "sigma_nutpie = np.sqrt(linalg.eigvalsh(Sigma_nutpie))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "de252391-17a7-400e-bd9e-6baf1349275b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(13462.75393777444)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lower, upper = linalg.eigvalsh(Sigma, np.diag(nutpie_target))[[0, -1]]\n",
    "np.sqrt(upper / lower)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b6aa84a0-bdd8-4231-bc42-38f6f9de1395",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stan condition number:  25223.072679706347\n",
      "nutpie condition number:  13462.753932914298\n"
     ]
    }
   ],
   "source": [
    "print(\"stan condition number: \", sigma_stan[-1] / sigma_stan[0])\n",
    "print(\"nutpie condition number: \", sigma_nutpie[-1] / sigma_nutpie[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8ac5b81e-8f70-4da9-9797-510acc9992ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVtpJREFUeJzt3Xd81PXhP/DXXcZlX8Zl57InI2yQDSpDsIJ+iyLIqF9paUFRflqltEoHRVv1S7WVom3BKgpaFXGWIZuwCStkEkggOyS5yx73/v3xSS45MkhILp9L7vV8PO5xuc997j7vezPulfdUCCEEiIiIiGSglLsAREREZL0YRIiIiEg2DCJEREQkGwYRIiIikg2DCBEREcmGQYSIiIhkwyBCREREsmEQISIiItnYyl2AjhgMBuTk5MDV1RUKhULu4hAREVEnCCGg1+sREBAApbLjNg+LDiI5OTnQarVyF4OIiIjuQnZ2NoKCgjo8x6KDiKurKwDpg7i5uclcGiIiIuoMnU4HrVZr/B7viEUHkabuGDc3NwYRIiKiPqYzwyo4WJWIiIhkwyBCREREsmEQISIiItlY9BiRzhBCoL6+Hg0NDXIXhajfsbGxga2tLafPE5HZ9OkgUltbi9zcXFRWVspdFKJ+y8nJCf7+/rC3t5e7KETUD/XZIGIwGJCZmQkbGxsEBATA3t6ev7UR9SAhBGpra1FYWIjMzExERUXdcWEiIqKu6rNBpLa2FgaDAVqtFk5OTnIXh6hfcnR0hJ2dHa5fv47a2lo4ODjIXSQi6mf6/K83/A2NyLz4b4yIzIn/wxAREZFsGERkMGXKFDz77LNyF6PHhIaGYuPGjXIX464oFArs3LlT7mIQEVktBpF+Ijc3FwsWLEBMTAyUSmW7Qeezzz7DgAEDoFKpMGDAAHzxxRfdvvapU6fw05/+tNvvY422bt0Kd3d3uYtBRCQbBpF+oqamBt7e3li7di2GDBnS5jkJCQl47LHHsGjRIpw/fx6LFi3Co48+ihMnTnTr2t7e3hwwTEREd8WsQWTDhg0YNWoUXF1d4ePjg7lz5yIlJcWcl+yTSkpKsHjxYnh4eMDJyQkPPPAA0tLSTM557733jDOEHn74Ybz55psmv0mHhobiL3/5CxYvXgy1Wt3mdTZu3Ihp06ZhzZo1iI2NxZo1a3DffffdsVvl2LFjmDRpEhwdHaHVavHMM8+goqLC5Not3yM5ORkTJkyAg4MDBgwYgL1797bqArl58yYee+wxeHh4wMvLC3PmzMG1a9eMzy9duhRz587F66+/Dn9/f3h5eWHFihWoq6sDAKxZswb33HNPq7LGx8fjlVdeASC11EybNg0ajQZqtRqTJ0/G2bNn2/2cBw4cgEKhQGlpqfFYYmIiFAqFSdnuVB+3O3/+PKZOnQpXV1e4ublhxIgROH36NA4cOICf/OQnKCsrg0KhgEKhwLp16wAAH374IUaOHAlXV1f4+flhwYIFKCgoaFXWffv2YeTIkXBycsK4ceP474uIOu/mGeDbXwKn/ilrMcwaRA4ePIgVK1bg+PHj2LNnD+rr6zF9+vQO/9PuDiEEKmvrZbkJIe663EuXLsXp06exa9cuJCQkQAiBWbNmGb90jx49iuXLl2PVqlVITEzEtGnTsH79+i5fJyEhAdOnTzc5NmPGDBw7dqzd11y8eBEzZszAI488ggsXLmDHjh04cuQIVq5c2eb5BoMBc+fOhZOTE06cOIF3330Xa9euNTmnsrISU6dOhYuLCw4dOoQjR47AxcUFM2fORG1trfG8/fv3IyMjA/v378f777+PrVu3YuvWrQCAhQsX4sSJE8jIyDCef/nyZVy8eBELFy4EAOj1eixZsgSHDx/G8ePHERUVhVmzZkGv13ep3rpTH01lDQoKwqlTp3DmzBm89NJLsLOzw7hx47Bx40a4ubkhNzcXubm5eP755wFI09N///vf4/z589i5cycyMzOxdOnSVu+9du1avPHGGzh9+jRsbW3x5JNP3vVnIyIrU5gKnNwMpHwrazHMuo7I999/b/J4y5Yt8PHxwZkzZzBp0qQev15VXQMGvPzfHn/fzkj63Qw42Xe9OtPS0rBr1y4cPXoU48aNAwBs27YNWq0WO3fuxLx58/D222/jgQceMH5JRUdH49ixY/j666+7dK28vDz4+vqaHPP19UVeXl67r/nzn/+MBQsWGMecREVF4a233sLkyZOxadOmVutK7N69GxkZGThw4AD8/PwAAOvXr8e0adOM52zfvh1KpRL/+Mc/jIvQbdmyBe7u7jhw4IAxLHl4eOCvf/0rbGxsEBsbi9mzZ2Pfvn1YtmwZBg0ahPj4eHz00Uf4zW9+Y6y3UaNGITo6GgBw7733mpRt8+bN8PDwwMGDB/Hggw92qe7utj4AICsrCy+88AJiY2ONr2miVquhUCiMddWkZaAIDw/HW2+9hdGjR6O8vBwuLi7G59avX4/JkycDAF566SXMnj0b1dXVXO+DiO5MNG6NopB3lEavXr2srAwA4Onp2ebzNTU10Ol0Jrf+7sqVK7C1tcWYMWOMx7y8vBATE4MrV64AAFJSUjB69GiT193+uLNuX31WCNHhirRnzpzB1q1b4eLiYrzNmDHDuLLt7VJSUqDVak2+WG8v65kzZ5Ceng5XV1fje3p6eqK6utqkhWPgwIGwsbExPvb39zfpnli4cCG2bdtm/Bwff/yxsTUEAAoKCrB8+XJER0dDrVZDrVajvLwcWVlZd6qmHqsPAFi9ejWeeuop3H///Xj11VdNPmN7zp07hzlz5iAkJASurq6YMmUKALQqe3x8vPFnf39/ADCpIyKidgmDdK+w6fg8M+u1lVWFEFi9ejUmTJiAQYMGtXnOhg0b8Nvf/vaur+FoZ4Ok382469d3h6Pd3f1Bttel0zIgtBUW7qYryM/Pr1XrR0FBQatWkpYMBgN+9rOf4Zlnnmn1XHBwcIfl7ug9R4wYYQwRLXl7ext/trOzM3lOoVDAYDAYHy9YsAAvvfQSzp49i6qqKmRnZ2P+/PnG55cuXYrCwkJs3LgRISEhUKlUGDt2rEn3T0tNC3e1rNum7rGWZe9KfQDAunXrsGDBAnzzzTf47rvv8Morr2D79u14+OGH2zy/oqIC06dPx/Tp0/Hhhx/C29sbWVlZmDFjRquyt6yjpnpvWUdERO0yWEaLSK8FkZUrV+LChQs4cuRIu+esWbMGq1evNj7W6XTQarWdvoZCobir7hE5DRgwAPX19Thx4oSxa6a4uBipqamIi4sDAMTGxuLkyZMmrzt9+nSXrzV27Fjs2bMHzz33nPHY7t27jddty/Dhw3H58mVERkZ26hqxsbHIyspCfn6+MeCcOnWq1Xvu2LEDPj4+cHNz6/LnaBIUFIRJkyZh27ZtqKqqwv33328Sqg4fPox33nkHs2bNAgBkZ2ejqKio3fdrCkG5ubnw8PAAIA1Wvb3sXamPJtHR0YiOjsZzzz2Hxx9/HFu2bMHDDz8Me3v7VjtHJycno6ioCK+++qrx7//d/HkTEXWoqUVE5tWTe+XqTz/9NHbt2oX9+/cjKCio3fNUKhXc3NxMbv1dVFQU5syZg2XLluHIkSM4f/48nnjiCQQGBmLOnDkApPr79ttv8eabbyItLQ2bN2/Gd99916rlITExEYmJiSgvL0dhYSESExORlJRkfH7VqlXYvXs3XnvtNSQnJ+O1117D3r17O1xc7cUXX0RCQgJWrFiBxMRE45iWp59+us3zp02bhoiICCxZsgQXLlzA0aNHjYNVm8q7cOFCaDQazJkzB4cPH0ZmZiYOHjyIVatW4caNG12qv4ULF2L79u349NNP8cQTT5g8FxkZiQ8++ABXrlzBiRMnsHDhQjg6Orb7XpGRkdBqtVi3bh1SU1PxzTff4I033uhWfVRVVWHlypU4cOAArl+/jqNHj+LUqVPGkBkaGory8nLs27cPRUVFqKysRHBwMOzt7fH222/j6tWr2LVrF37/+993qV6IiO7I2DUj80oewowMBoNYsWKFCAgIEKmpqV1+fVlZmQAgysrKWj1XVVUlkpKSRFVVVU8UtVdNnjxZrFq1yvj41q1bYtGiRUKtVgtHR0cxY8aMVvX17rvvisDAQOHo6Cjmzp0r/vCHPwg/Pz+TcwC0uoWEhJic8+mnn4qYmBhhZ2cnYmNjxWeffXbH8p48eVJMmzZNuLi4CGdnZxEfHy/Wr19vfD4kJET83//9n/HxlStXxPjx44W9vb2IjY0VX331lQAgvv/+e+M5ubm5YvHixUKj0QiVSiXCw8PFsmXLjH/WS5YsEXPmzDEpx6pVq8TkyZNNjpWUlAiVSiWcnJyEXq83ee7s2bNi5MiRQqVSiaioKPHpp5+2KisA8cUXXxgfHzlyRAwePFg4ODiIiRMnik8//VQAEJmZmZ2uj5ZqamrE/PnzhVarFfb29iIgIECsXLnS5O/t8uXLhZeXlwAgXnnlFSGEEB999JEIDQ0VKpVKjB07VuzatUsAEOfOnRNCCLF//34BQJSUlBjf59y5c63K2hP68r81IurA8b8L8YqbEJ8s7fG37uj7+3YKIbox7/QOfvGLX+Cjjz7Cl19+iZiYGONxtVrd4W+mTXQ6HdRqNcrKylq1jlRXVyMzMxNhYWFWOUNg2bJlSE5OxuHDh+Uuyh0dPXoUEyZMQHp6OiIiIuQuDnWRtf9bI+q3Et4B/rsGGPRj4Mc9u5ZIR9/ftzPrgIpNmzYBgHHEf5MtW7a0uSYCte/111/HtGnT4OzsjO+++w7vv/8+3nnnHbmL1aYvvvgCLi4uiIqKQnp6OlatWoXx48czhBARWRLjGJF+PGvGjI0tVufkyZP405/+BL1eb1xX4qmnnpK7WG3S6/X45S9/iezsbGg0Gtx///2txloQEZHMLGQdkb41xcSKffLJJ3IXodMWL16MxYsXy10MIiLqiIWsI8JN74iIiKyRcR2Rjtd+MjcGESIiImvUNHxC5jEiDCJERETWyELGiDCIEBERWSMLWdCMQYSIiMgaGceIsGuGiIiIehtbRKzXlClTOtzfpa8JDQ3Fxo0b5S7GXVEoFNi5c6fcxTCrpUuXYu7cuXIXg4gsTdMYEQ5WpZ6Qm5uLBQsWICYmBkqlst2g89lnn2HAgAFQqVQYMGAAvvjii25f+9SpU/jpT3/a7fehjt1taPrLX/6CrVu39nh5iKiPM7aIcPou9YCamhp4e3tj7dq1GDJkSJvnJCQk4LHHHsOiRYtw/vx5LFq0CI8++ihOnDjRrWt7e3vDycmpW+9B5qNWq+Hu7i53MYjI0hi4oBk1KikpweLFi+Hh4QEnJyc88MADSEtLMznnvffeg1arhZOTEx5++GG8+eabJl8uoaGh+Mtf/oLFixdDrVa3eZ2NGzdi2rRpWLNmDWJjY7FmzRrcd999d+xWOXbsGCZNmgRHR0dotVo888wzqKioMLl2y/dITk7GhAkT4ODggAEDBmDv3r2tfpu/efMmHnvsMXh4eMDLywtz5szBtWvXjM83dSe8/vrr8Pf3h5eXF1asWIG6ujoAwJo1a3DPPfe0Kmt8fDxeeeUVAFJLzbRp06DRaKBWqzF58mScPXu23c954MABKBQKlJaWGo8lJiZCoVCYlO1O9XG7devWYejQofjggw8QGhoKtVqN+fPnQ6/Xt1uHADB06FCsW7fO+DwAPPzww1AoFMbHTe+9efNm49+PefPmmXyG27tmhBD405/+hPDwcDg6OmLIkCH4z3/+0275iaif4hgRMxACqK2Q59aNfXWWLl2K06dPY9euXUhISIAQArNmzTJ+6R49ehTLly/HqlWrkJiYiGnTpmH9+vVdvk5CQgKmT59ucmzGjBk4duxYu6+5ePEiZsyYgUceeQQXLlzAjh07cOTIEaxcubLN8w0GA+bOnQsnJyecOHEC7777LtauXWtyTmVlJaZOnQoXFxccOnQIR44cgYuLC2bOnIna2lrjefv370dGRgb279+P999/H1u3bjV2MSxcuBAnTpxARkaG8fzLly/j4sWLWLhwIQBpz5slS5bg8OHDOH78OKKiojBr1iyTANBVXa2PJhkZGdi5cye+/vprfP311zh48CBeffXVTl/31KlTAKQNI3Nzc42PASA9PR2ffPIJvvrqK3z//fdITEzEihUr2n2vX//619iyZQs2bdqEy5cv47nnnsMTTzyBgwcPdro8RNQPWMgYkf6110xdJfDHAHmu/ascwN65yy9LS0vDrl27cPToUYwbNw4AsG3bNmi1WuzcuRPz5s3D22+/jQceeADPP/88ACA6OhrHjh3D119/3aVr5eXlwdfX1+SYr68v8vLy2n3Nn//8ZyxYsMA45iQqKgpvvfUWJk+ejE2bNrXaFn737t3IyMjAgQMH4OfnBwBYv349pk2bZjxn+/btUCqV+Mc//gFFY9/kli1b4O7ujgMHDhjDkoeHB/7617/CxsYGsbGxmD17Nvbt24dly5Zh0KBBiI+Px0cffYTf/OY3xnobNWoUoqOjAQD33nuvSdk2b94MDw8PHDx4EA8++GCX6u5u66OJwWDA1q1b4erqCgBYtGgR9u3b1+lA6e3tDQBwd3c31muT6upqvP/++wgKCgIAvP3225g9ezbeeOONVudWVFTgzTffxA8//ICxY8cCAMLDw3HkyBFs3rwZkydP7lxFEFHfxxYRAoArV67A1tYWY8aMMR7z8vJCTEwMrly5AgBISUnB6NGjTV53++POUtw2KEkI0epYS2fOnMHWrVvh4uJivM2YMQMGgwGZmZmtzk9JSYFWqzX5Ary9rGfOnEF6ejpcXV2N7+np6Ynq6mqTFo6BAwfCxqY5qfv7+6OgoMD4eOHChdi2bZvxc3z88cfG1hAAKCgowPLlyxEdHQ21Wg21Wo3y8nJkZWXdqZp6rD6ahIaGGkNIW5+lO4KDg40hBADGjh0Lg8GAlJSUVucmJSWhuroa06ZNM/kM//73v03qnoisgIVsete/WkTsnKSWCbmufRdEO106LQNCW2Ghvdd1xM/Pr1XrR0FBQatWkpYMBgN+9rOf4Zlnnmn1XHBwcIfl7ug9R4wYYQwRLTX95g8AdnZ2Js8pFAoYmgZXAViwYAFeeuklnD17FlVVVcjOzsb8+fONzy9duhSFhYXYuHEjQkJCoFKpMHbsWJPun5aUSqXxMzRp6h5rWfau1EdnP4tSqWz1Z3r7tTurqf7b+nNouuY333yDwMBAk+dUKtVdXY+I+iiDZSzx3r+CiEJxV90jchowYADq6+tx4sQJY9dMcXExUlNTERcXBwCIjY3FyZMnTV53+vTpLl9r7Nix2LNnD5577jnjsd27dxuv25bhw4fj8uXLiIyM7NQ1YmNjkZWVhfz8fGPAaTmeoek9d+zYAR8fH7i5uXX5czQJCgrCpEmTsG3bNlRVVeH+++83CVWHDx/GO++8g1mzZgEAsrOzUVRU1O77NYWg3NxceHh4AJAGq95e9q7UR2d5e3sjNzfX+Fin07VqYbGzs0NDQ0Or12ZlZSEnJwcBAVK3ZEJCApRKpbGLqqWmqdtZWVnshiGydk0tIkp2zVi1qKgozJkzB8uWLcORI0dw/vx5PPHEEwgMDMScOXMAAE8//TS+/fZbvPnmm0hLS8PmzZvx3XfftfqNNzExEYmJiSgvL0dhYSESExORlJRkfH7VqlXYvXs3XnvtNSQnJ+O1117D3r17O1xc7cUXX0RCQgJWrFiBxMRE45iWp59+us3zp02bhoiICCxZsgQXLlzA0aNHjYNVm8q7cOFCaDQazJkzB4cPH0ZmZiYOHjyIVatW4caNG12qv4ULF2L79u349NNP8cQTT5g8FxkZiQ8++ABXrlzBiRMnsHDhQjg6Orb7XpGRkdBqtVi3bh1SU1PxzTff4I033uhWfXTWvffeiw8++ACHDx/GpUuXsGTJEpNuKUDq3tm3bx/y8vJQUlJiPO7g4IAlS5bg/PnzOHz4MJ555hk8+uijrcaHAICrqyuef/55PPfcc3j//feRkZGBc+fO4W9/+xvef//9bn0GIupjOEaEmmzZsgUjRozAgw8+iLFjx0IIgW+//dbYnD9+/Hj8/e9/x5tvvokhQ4bg+++/x3PPPddqYOSwYcMwbNgwnDlzBh999BGGDRtmbA0AgHHjxmH79u3YsmUL4uPjsXXrVuzYscNkfMrt4uPjcfDgQaSlpWHixIkYNmwYfvOb38Df37/N821sbLBz506Ul5dj1KhReOqpp/DrX/8aAIzldXJywqFDhxAcHIxHHnkEcXFxePLJJ1FVVdXlFpJ58+ahuLgYlZWVrVYP/de//oWSkhIMGzYMixYtwjPPPAMfH59238vOzg4ff/wxkpOTMWTIELz22mv4wx/+0K366Kw1a9Zg0qRJePDBBzFr1izMnTsXERERJue88cYb2LNnD7RaLYYNG2Y8HhkZiUceeQSzZs3C9OnTMWjQILzzzjvtXuv3v/89Xn75ZWzYsAFxcXGYMWMGvvrqK4SFhXXrMxBRH2MhY0QU4m4GG/QSnU4HtVqNsrKyVl9Q1dXVyMzMRFhYWLszFfqzZcuWITk5GYcPH5a7KHd09OhRTJgwAenp6a2+XKl71q1bh507d7bqQupJ1v5vjajf+mwZcPETYPp6YFzHSxB0VUff37frX2NE+rHXX38d06ZNg7OzM7777ju8//77Hf7WK6cvvvgCLi4uiIqKQnp6OlatWoXx48czhBARWRLjGBHOmqFOOHnyJP70pz9Br9cjPDwcb731Fp566im5i9UmvV6PX/7yl8jOzoZGo8H999/faqwFERHJTFjGrBl2zRBRh/hvjaif+mQxkPQlMOt1YPSyHn3rrnTNcLAqERGRNbKQdUQYRIiIiKxRU4cIg0j3WHDPElG/wH9jRP2UhWx612eDSNMaG5WVlTKXhKh/a/o3dvsy9UTUx1nIgmZ9dtaMjY0N3N3djRuHOTk53XGPEyLqPCEEKisrUVBQAHd391YrvRJRH2ccI8Lpu3etaQnrntrFlIhac3d3b3O5eCLq49gi0n0KhQL+/v7w8fG5651Kiah9dnZ2bAkh6q8sZIxInw4iTWxsbPifJRERUVcYZ83IO6yhzw5WJSIiom6wkDEiDCJERETWyELGiDCIEBERWSMLGSPCIEJERGSN2CJCREREsuEYESIiIpINW0SIiIhILrqqGgDA3w9lyloOBhEiIiJrU1MOJ/11AMCNOhdZi8IgQkREZG1unoatoRo3hAY5qghZi8IgQkREZG2qSgEAN4UGSiXHiBAREVFvqtEBAMqFI5Qyb1zPIEJERGRtavQAAD0cYSNzEmEQISIisjaNQaRcOELJIEJERES9qimIwAk23H2XiIiIelXjGBG9YNcMERER9bbGWTPlcISSLSJERETUayqKgZTvAAAlwgU2MicBBhEiIiJrcnYrYKhDqUMQ9hpGsEWEiIiIelHWcQDAuYDHUQkHzpohIiKiXlSUCgAodJKWduesGSIiIuoddVVAibTZXZEqBAA4a4aIiIh6yY3TAATg4gu9rQcAcIwIERER9ZITf5fuQyfCIKQfOWuGiIiIzK84A0j+GoACuOfnMAgpibBFhIiIiMwvfa90HzYJCBqJBoP0sF/Pmjl06BB+9KMfISAgAAqFAjt37jTn5YiIiKg9t65K9/5DAMDYItKvZ81UVFRgyJAh+Otf/2rOyxAREdGdNM6WgYc0W6ahcZCI3C0ituZ88wceeAAPPPCAOS9BREREd2IwAAWXpZ/dQwEADRbSImLWINJVNTU1qKmpMT7W6XQyloaIiKgfKEwFvloFlGYBdk7NXTONLSKcNdPChg0boFarjTetVit3kYiIiPquiiLg48eArGOAjT0w9x3AxRtA8xgRubtmLCqIrFmzBmVlZcZbdna23EUiIiLqm+qqgPfulQapOqiB5UeBgQ8bnzbOmmHXTDOVSgWVSiV3MYiIiPq+hL8BpdcBRw9g8S7AO9rkaUuZNWNRQYSIiIi6SQhg/x+BQ3+SHk9dC/jHtzrNKmbNlJeXIz093fg4MzMTiYmJ8PT0RHBwsDkvTUREZH0MDcCel4GExmUzhj4BjPzfNk9tnjXTW4Vrm1mDyOnTpzF16lTj49WrVwMAlixZgq1bt5rz0kRERNbl1lVg5y+ArATp8ZjlwPQ/AMq2h4M2z5rpxy0iU6ZMgWhMXERERGQmF/8D7Pw50FAL2LsAs98Ahszv8CWWMmuGY0SIiIj6spPvAd8+L/0cMh6Yu8m4empHOGuGiIiI7p4+T1qoLPV76fGwRcCDGwGbzn21c9YMERERdZ0QwOl/Ad+vARpqAIUNMOZnwIw/Al0IFVYxa4aIiIh6iKFBav1I+Btw/ah0zGcg8ON/Aj5xXX87YRlLvDOIEBERWbqqEuDLlUDy19JjWwdg0vPAxOe71ArSkrFFhF0zRERE1CaDAbj4CfDdL4HqMunYPb+Qbu7d24+twRqm7xIREdFduvI18MMfgMIr0mPvWGDKGmDg3B55+6bVNdgiQkRERM2qdcC+3wGn3pMeq9yA8auA8c92ekZMZzStrMogQkRERNJ03KN/kWbE1FdLx8avAiasBhzde/xy7JohIiIiQJcDHPsrcPqfzQFEHQzM3ADEPWi2y3LWDBERkTW7cRo49DqQ+l3zsaBRwKingPjH7no2TGdx1gwREZG1aagDUv8LXNkl7Q8jGqTj3rHAxP8HDJ5n9gACSBveFehrALBrhoiIqP8rLwQufSYNQC1Obz4+8BFg6q8ATVSvFudsVgkK9TWwUSoQ4+vaq9e+HYMIERGROdRVSSuhnt8OpO8FDPXScScNEP8oED0TCJvUKy0gt2tqDRmqdYePm0OvX78lBhEiIqKe0lAHXD0InN0qdcE01DY/FzAMiJ8PDF0AOLjJVkQA0FfXAQDUjnaylgNgECEiIuq+W1eBU/8Ezn0IVJc2H1cHA/HzpADiHS1b8W6nr5ZaZ1wd5I8B8peAiIioL6ooBlK+BS79B7h6oPm4o4c06HTYIsBvsCxdL3eiYxAhIiLqg+prgPzLUsvH2febx31AAUTeB4xaBkRNA5Q2shbzTpq6Zlwd2DVDRERk2YQAshKkWS8XP23efA4AvKKAgQ8Dw54APELkK2MXsWuGiIjI0tXXAuf+DVz8DMg61nzcQQ0Ej5V2wJVp1kt3peTpAbBFhIiIyPKUXAMO/gm48AlgkLowoLSVxn0MngeETwWUMq+L3g1p+XpcvCm16si9hgjAIEJERAQYGoAzW4DEj4Gcc80rnjqogbErpYXHNJHylrGHHEorAgCMCPHA6DBPmUvDIEJERNasvBBI2gmc/TeQd6H5eMS9wKQXgMCRgK29bMUzh4QMKYhMH+Arc0kkDCJERGQ9hABKMoFrR4DLX0iLjzW1fqjcgMkvAgMeAtyD5S2nmdQ3GHDi6i0AwNgIL5lLI2EQISKi/q9pr5eTm6XFx1oKGAYM+rG07LqLjzzl6yWXc3TQ19TD1cEWAwPUchcHAIMIERH1R3XVwM0zQNp/gesJwI1TAKRt76G0AwKHA5HTgEGPAF4Rsha1N20/lQUAGBPmJfuuu00YRIiIqO8zGICCJCD7OHDxP1IIabnPCwD4DgJGPgnEPwaoXOQpp4wqa+ux41Q2AODRkUEyl6YZgwgREfU9DfVA3nng2lHg+jFpnY+WC40BgIsvEDIOiJoOhE7ot+M+OispRweDAHzdVJg+0E/u4hgxiBARkeWrrwGuHwVunJFWOc0+AdSWm55j7wL4xAGxDwJxPwI8w/vkYmPm0rR2yOBAyxgb0oRBhIiILEtDHVCcARRcBrJPAhn7gaJUGMd4NHFQA8HjgNDxQMh4wC8esOHXWnuagoilDFJtwj8xIiKSV41eGtORcw64cRpI2wM01LQ+z0EtDTDVjpa6XHwG9ukVTnuTrroOB1IKAbBFhIiIrF15AZB1XOpiuX5MWkhMGEzPsXcBvGOlqbWhE6S9XVx82NVyl3Yl5uBWRS00LircYyHrhzRhECEiIvOpqwZyz0vTZ2+cAnLOAqVZrc9zDwYCRwD+Q4HQidL0WoaOHlFWWYf/nLkBAFh0TwhcVJb11W9ZpSEior6rvlYa15F7ASi4Atw4Kf3ctHGckUJq7QgZJ92CxwLqQFmK3N9lFJbj8XePo0BfA6UCmGYhy7q3xCBCRERdV1spdalcPwrkXZJWKy1Iar12BwA4ewNBo4GgkVKrR8BQabwHmVVxeQ1WbDuLAn0NQr2csGZWHAYEuMldrFYYRIiIqGM15dJA0vzLQN5FqXulMLn1uA4AcHCXgoYmWgof2lGAewi7WXpZXlk1Fv3zBNIKyuHhZIdPfjYWPm4OcherTQwiRETUTAipWyX3PHArQ1qv43pCG90rkBYMCxoFaMcAHqGA32DpnqFDVldydXjor0dQ1yDg46rCR8vGWGwIARhEiIisV32ttD5HUQpQmArkX5JmsVTdan2uo6c0lsN3oDSTJWAo4BbQ60Wmjumr6/CTLadQ1yCtufKPJSMR6eMqc6k6xiBCRGQNhABKrgFXD0jdLIUpUstHTVnrc23spVYOr0hppdLI+61qY7i+at+VfLyy6zLydNUAgG1PjUF8kLu8heoEBhEiov7E0CAFjqK05paOpvu2QodKDXjHAN7R0riO4LGA/xDAVtXrRae7k15Qjj98k2RcsMxf7YANjwzG+EiNzCXrHAYRIqK+qK4KKE6XulZaho3i9LZXJQUAhY0UNLSjpS4WjzApdHBZ9D5JCIEdp7Lx26+SUFXXAAB4akIYVk+PhpN93/kz7TslJSKyRg31QFm2FDAKU6TAkXNOmjJ7+94rTWwdAK+oxlaOmOZ7z3DAznIHLVLn5euq8fKXl/Dfy/kAgPGRXvjVrDiL20emMxhEiIgsQdO6HLnnpTU5bmVK9yXX2p6xAkhTZb1jpC4V75jm0KHWAkqb3iw99RKDQeDbS7n4zc5LKKmsg52NAi/MiMFTE8KhVPbN2UoMIkREvam6DChIllo2bmUCJZnSeI6CK4BoaPs1NiqpNaOpZcMrAgibBLj6c6qsldBX1+FIWhHe/iEdSbk6AMAAfzf8eV58n2wFaYlBhIiopwkB6HMbp8amSS0dpVlAUTqgu9H+61x8gYDhgCYK8AyTxnB4RQJugdxl1oodSCnAczsSUVIptYy5qmwxf7QWz03rW2NB2tP3PwERkVzqa4DiDKA4rTl0NN3Xlrf/OrdAqTvFM1wKHJ4RgH+8dJwtHNQo+1YlXv0+Gd9cyAUA+Lk54KGhAVg+OQKezvYyl67nMIgQEXWkaf2N0uvSfUGyNHC0OF061tYy54A0Q8UzvHHF0UHNwcM7FnB0773yU59TVlWH9w5dxT+OXEV1nQEKBfDEmBCsnR0HB7v+N/aHQYSICJACR3m+tJdK3kXTqbFtrb/RROUmdZ94x0hdKprGcRweoYBt//mtlcyvsrYe/72ch7d/SMfVwgoAwOgwT6x5IBbDgj1kLp35MIgQkXVpqAcqi5uDRlEakH9Rmq1S3U7gsLGXgoU6CPCOkwKHV4QUOFx82J1C3VLXYMBb+9Kw9eg16GvqAQAaF3v85sEBeGhIABT9/O8XgwgR9U+1ldLMlIJkaafYwmRpZkppFtpdf0OhbB6v4R3boqUjGrCx69XiU/+XkqfH3iv52H4qC9m3qgAAWk9HzBuhxeKxIXB3so4WNQYRIurbaisbu1Aag0ZT6Ci5jnYDBxSAu7a5G8U7BggcLi0CxgW/yIyEEDh9vQR/2ZuGI+lFxuPO9jZ4blo0fjI+DDZ9dD2Qu8UgQkR9Q12VtLJoYQpQeKWxpeNKx4HDyUvqSvGJlVo4fOKksOHkxWXNqVfllVXju0u5+OT0DVxpXAfEzkaBkSGeeHh4IGYP9oezyjr/TlrnpyYiy1VXJbVwNAWNpl1iS67hjoHDO0YKG02hw7lvbPpF/Zeuug6bD2bgH4czUVMvzbBytLPBg/H+ePreKAR7OclcQvkxiBCRPOqqpIGirbpUrrU/JdbR0zRoeMdIAcTFu1eLTnQnyXk6/OtIJnYn5aO0cSGyoVp3zBrsh0dHaq1m/EdnMIgQkXkZGhpXF714W5fKtc4FDu/Yxq6VxhaOfj6DgPquzKIKfH0+B19fyEVKvt54PNLHBS/MiMH0Ab79fgbM3eiVIPLOO+/gz3/+M3JzczFw4EBs3LgREydO7I1LE1FvMRikTdoKk6X9UwpTgIIkqbWjrrLt1zh6tBjD0WIsh7M3Awf1CbrqOnx7IRffXsrDodRC43E7GwXujfXBwjEhGBfhBVsbLtHfHrMHkR07duDZZ5/FO++8g/Hjx2Pz5s144IEHkJSUhODgYHNfnoh6mqFB2qwt/xJw45R0X1EshZC6irZfY+cM+A0GfAe0aOWIY+CgPklfXYfvL+Vh1/kcJGQUo97QPHZpYpQGPxoSgBkD/KB24pTvzlAIIdqb39YjxowZg+HDh2PTpk3GY3FxcZg7dy42bNjQ4Wt1Oh3UajXKysrg5uZmzmISUVvKC5sHjOZdBPIvS60c7bVw2DpIIcMzXFqDw3cA4DNQWvyL29JTH1Zd14AjaUX48nwOvr+Ui7qG5q/OcI0z/mdEEKbEePf5nXB7Sle+v83aIlJbW4szZ87gpZdeMjk+ffp0HDt2rNX5NTU1qKmpMT7W6XTmLB4RNamtbBw0miSFjabAUVHY9vm2jlKLhu9AQDtG2o7eXSstBsZpsdRP1DUY8ENyAXadz8H+5AJU1jYYnwtQO2DBmGDMjg9AmMZZxlL2fWb9H6OoqAgNDQ3w9fU1Oe7r64u8vLxW52/YsAG//e1vzVkkIutWW9G4Q2xai0XAkqQdZNucGquQljb3jpFCh+8g6cYWDurHTl27hU9OZWN/SiGKypt/OfZ1U+H+OF/8eERQv977pbf1yq8ut48SFkK0OXJ4zZo1WL16tfGxTqeDVqs1e/mI+h1DgxQusk9IrRtFKVL4KMtu/zVOXlLY8BnY3KXiEwvY87c96v9q6htwLL0Ymw5m4GTmLeNxjYsKjwwPxIPx/hgcqOasFzMwaxDRaDSwsbFp1fpRUFDQqpUEAFQqFVQqlTmLRNT/VJdJ4zeK0oCs40BxGpCfBNRXtX2+k5e0rHnTTrFNoYObt5GVqa5rwL4rBfjqfA4Ophaiqk7qerG3UeJ/RgRi+kA/TIjUwI4zXszKrEHE3t4eI0aMwJ49e/Dwww8bj+/Zswdz5swx56WJ+h9DQ/P02KI0acbKzTPS1vVtsXMGvMKBsMlS4PCOkZY3d/bq3XITWZDymnqcyyrBd5fy8OW5m6hoMe7D21WFH8UH4KeTwuGn5p5DvcXsXTOrV6/GokWLMHLkSIwdOxbvvvsusrKysHz5cnNfmqhvEgKoKmlcBOyCND0272LH63G4BwMeYUDwPYDPAKmLxTOc4ziIANwsrcKx9CJ8mZiD41dNp9sGeTjiR0MCMHuwPwYGuLHrRQZmDyKPPfYYiouL8bvf/Q65ubkYNGgQvv32W4SEhJj70kSWzWAAdDeB4vTmW2EykJMIVJe2/Rpbx+Zt6X3igNCJUheLo3svFpzIsmXfqsTZrBIcTS9CwtViZN8y7aYM9nTCEK07FowOxj3hngwfMjP7OiLdwXVEqF8QAqi8Ja02euuq1MLRNJC0vemxgDQl1neQtBCY3yDAdzBnqxDdpsEgkF5QjlPXbuHCjVKcvlaCq0WmC+vZKBWID1JjUpQ3Hh4WiFBOtzU7i1lHhMjqGBqkPVTyLzXOVkkFsk9KLR9tUdo2L/7lFSGN4fCPlwaT2nNXTqLbpReUI7OoApdzynDmeglOXysxDjJtYqtUYECAG8aEeWJcpAajQj3houLXnaXinwzR3RACqChq3jk274IUPAqT2x/H4RoAeIY1L3HuFQGETgBsOVOMqC1F5TVIydPj0s0ynLpWgvM3SlGor2l1npO9DYZq3TEixAMDA9wwLlIDNwcur95XMIgQ3YnBAOhzpKCR/E1j98ploOpW2+e3XHW0KXSEjGMLB1EHhBC4WlSB41eLcTS9CEk5Olwrbh3qbZQKDPB3Q5jGGSNDPTA82ANx/m6wUXKcR1/FIELUUm2lNGi0KBXISgByG1s62tzMreWqo4Max3EM4mwVok4oKq/BD8kFSMrR4UquDin5epRW1pmco1BIA0sHBrhhSJA7RoZ6ItLHBWpHtnb0JwwiZL2EkMZu5CQCN08D1xOAnLNAQ23rc5V2gFsAED0DCBwpzVTxiQPsHHu92ER9TUlFLa7k6pCUq0NidikSs0txo6T1gnsqWyWGaN0xKUqDoVoPDAp0g7uTvQwlpt7EIELWpewmkHlIul3dD+hzW5/j6CFt3qYdDQSOaNxbJZKbuRF1ghACebpqnM8uw9msEuy7ko+MwrZaFIEB/m4YH+mFGD83xPq5IsrXBSpbtiZaG/7PSv1bXTWQ8YO0CmnSl8CtDNPnFTbSAmABQ4DgsYD/ECl4cF0Bok6pqW/AmWslOJhWiKQcHa4WVuBmaevWjmBPJ8T5uyI+yB1Dte4YHKTmgFICwCBC/U21Tpoum3MWSN8rdbs0tBhlr1ACAcOAsEnS0ufB97B7hagT6hsMuFZcgeQ8PVIab2kF5bheXAFDG6tRxfq5YqjWHWPCPXFvrC/HdVC7GESob9PnAwWXgYz9wLUjQG4iIAym57gFSiuQRkwFYh4AHNSyFJWoLxBCoEBf0xg4dMbgkVZQjtp6Q5uv8XS2x72xPhgZ4oEQL2cMDHRjawd1GoMI9S1lN6XZLGm7pfvSrNbneIRJi4JF3AeEjJfW62BXC1ErFTX1SMlvbuFIbgwet89eaeJoZ4NoP1fE+roi2s8VMb7SuA4fVxWXSae7xiBClq26DMg8DFw9IA0uLU43fV6hlDZ8Cx4HhE8BQscD6iA5SkpksW7vVrmSq0dKvq7VHixNlAogVOOMWD9XxPq5IcbPFbF+rtB6OEHJ9TqohzGIkGWpq5a6V9L2ACnfSSuVihbLNytspHU7omdI4zwCRwIO3IeICGijWyVXj+Q8PdIL2+9W8XZVNQYOV+PslUgfFzjYcfYK9Q4GEZJf1gkg5Vtp4bBrR4D6235L84yQxneETwXCJnKMBxGA8pp6pObrkZzbYixHG4uCNXGyt0G0b1PgcG1s5XCDpzPX6SB5MYhQ76qvBfIuAtcOS1Nqcy8AZbeN83BwlwJHzGypq8U9WJaiElmCpm6VK7lN4zju3K0SpnE2dqnE+Lkizs8NQR6O7FYhi8QgQubVUCdNoc08KIWPrBOtWzygAGJnS60e/kOlRcQ48I2sTINBIKe0ChmF5S0Gj3bcreLjqjKO32gKHuxWob6GQYR6lqFB2ok285A0yDQrAagtNz3HyQsIGA6ET5YWEPMbLK1mSmQFmlYebdrOPi2/HGkF+g5nqzR1q8T5SzNVmsZyeLBbhfoBBhHqHiGkrpas49KslutHpZkuLTl6SNNowyZLXS7esWzxoH5PCIF8XQ2uFpbjws0yXM7R4VpRBTKLKlBeU9/ma+xtldB6OCLW3w2xvs3jONitQv0Zgwh1TV2VtHJp5iEpfBRcBqpKTM9RuTUGj4nSQmK+gwClUp7yEvWC6roGJOXqkJ5fjks5ZUjMLkVGQTkqahvaPN9GqUColxNCvZwR5euKaF8XRPq4IM7fDXY2/LdC1oVBhDpWXwPcPCsFj2uHpRDScsl0ALBRASFjpXU8wiYBfkO4QRz1S3UNBmQUliM5V4/0gnJcLSrH1cIKpBeUo76Ndc5tlAoEezohyscFw0M8EK5xRpjGGcFeTtzcjagRvy3IlKEBKLgCJH8t7dWSex5oqDU9x9VfChyhE6TBpZpowM5BluISmYMQArll1cYBo8l5OqTk6ZFRWI66hjY2VgGgcbFHrJ8bIrydMSbcC9G+Lgj2dIa9LVs4iDrCIGLtaisbp9MeAq4fA26cBmp0puc4aZq7WcImc8l06leq6xqQXlCO1HxpHY7z2aW4eKOs3W4VF5UtYv2kJc7DNc6I8Ja6VYI8HLnMOdFdYBCxNg31QOl1afXSpC+lFUzrKk3PsXcBgkYCg34srePhEcbgQX1eTX0DrhVV4sKNUlzO0eF6cQWuFlUg61YlRBuNHLZKBcK9nY0zVJoWAgt0Z+Ag6kkMItZACGnF0qv7gfM7AN0N0+ddfKWl0iOmAsH3AD4DACX7r6lvqm8w4EZJFc7fKEVynh75umpkFFYgKaes3W4VT2d7RPm4INrXFQMC3DA82ANhGnarEPUGBpH+qmla7eUvpFtJZvNzNvaAVxQQPR0YMEca58Hf8KgPKqusQ3qhHpdu6pCQUYwbpZVIztW3OXAUAFxVtojxc8VQrTvCvV0QqnFCtK8rNC6qXi45ETVhEOlPDAYg/yJw9SBw4RPp5yZ2TtJGcTGzgNgHAXsn+cpJ1EW66jok5eiQlq9HWkE50vLLkV5YjkJ9TZvn29soMSDADfFBavirHRHo4YihQe7QerJbhcjSMIj0dRVF0gyXjB+klUyrbjU/Z6MCoqYBgx4BomcC9s7ylZOoE+oaDLhWJI3duFpYgUs5Zbh8swzXiivbfU2A2gERPi4YFuyBGF9XxAepEejOBcCI+goGkb5Inw+c+0CaXpt9EhAtRvfbu0iLiUXcC8Q/Cjh5yldOonbUNxhw/VYl0vL1SM0vR0q+Hmn5emQWVbQ7jiPQ3VHaot7XBVE+rojycUGEjwtcVPxvjKgv47/gvsJgALKOAZd3SiGkvrr5Of+h0qZxYZOBwOGAjZ1cpSQy0WAQyL5VidR8feNNmiZ7tbACtQ1tb+TmbG+DCB8XBHs6YWCAGoMC3TAwQM3t6on6KQYRSyYEcOOUNMX2/HagLKv5ucCRwLCFUsuHR6hsRSRqoquuw6UbZcgtq8bFm2W4kqvDpZvtr8fhaGeDqMbWjWhfacZKtJ8rAtQOHMdBZEUYRCyRwQAkfwUcel3aybaJSg1ETAFGPim1fvA/a5JBfYMB14orkJpfjqQcHS7llCE1T4+csuo2z1fZKhHZODU2ytcF0T7N63FwHAcRMYhYkuoyIPFjIPFDaeotIM12CZ8qDTiNnQ3YOcpbRrIqlbX1uHRThws3SpGYXYqkHB1ulFS1260S5OGIIA9HxPm7YVCAGjF+rojzd4MNAwcRtYNBxBIUpQMnNkkhpK5COmbvCtzzc+nGAadkZgaDwM3SKlzO0SE5Tyctd54nDR5ta0kOZ3sbaD2dMFTrLoWOQDeEaVw4joOIuoxBRC5CAJkHgRPvAqnfN8988Y4Fhi8GBj8KuHjLW0bql/TVdUjNL8ex9CJcuFmG68UVuF5ciZr6tls5fN1UGBLkjiFadwwOVCNM48xuFSLqMQwiva2qFDj1nrTUenFa8/HQicCkF6RdbTn2g3pAfYMBqfnlOJFZjCu5OlwrrkRmUUW7i4DZKhWI9pW6UmL9pPEccf5u8HXjzspEZD4MIr2hvkaa+XLqH8C1w4ChXjqutJNmvgxZAGhHM4DQXdNX1yG9oBxXCytw8WaZcWO39lo5/NwcEOvviolR3ojwdja2ctjacG8VIupdDCLmVFEMHHwNuPw5UFHYfFwdDEx4Fhj4MMd/UJcYDAJXiyqQlq/HuexS5OuqjVvYt7UQmKvKFoOD1BgZ6okQTydE+bogxNMZaieuNUNEloFBxByEABK3Abt/07zkuosfMOQxYPgSwCtC3vJRn1BT34C0ximySbnSINKMwva7VnxcVQjTOCPO3w1DtGrEB7kjzMuZYzmIyKIxiPS0shvAt78EUr6RHvsOAqb+CoiazhVPqU0NBoG0Aj0yCiqQWSR1r1zJk5Y8b2sXWQc7JUK9nDEs2B0hXs4I1zgjytcVoV5OXAiMiPocBpGeYjAAabuBr54ByvOl8R/3/hoYu4IBhIwMBoGsW5W43LgQ2KWbZTifXQpddX2b56sd7TDA3w0DA9wQ5++GEC8nDApUw8HOppdLTkRkHgwi3WUwSDvf7vtt8yqo3nHAI+8C/vHylo1kV1JRi8s5OhxKK8TJzFtIzdejso0lz53tbRDt54pwjQvCNE6I9nXFwEA1lzsnon6PQaQ7DAbg82XApf9Ij20dpeXX710L2DvLWzbqddV1DTiSVoS9V/Jxs7QKKXl6FLQxnsPeVok4P1cMaNzQbXCgGgP83ThjhYisEoPI3SrNlkJIVgIAhTQNd/yzgCZK7pJRL6ipb0BqXjkSb5TiVOYtXLhRiuu3KiHaWIU0yMMRo8M8MTnaGwMD1Aj1cmLoICJqxCByN/IuAtsXAKVZgJ0zMPsNYOjjcpeKzKS23oDUfL1xg7eTmbeQXlDe5kBSjYs9HowPQJy/KyIbN3dzUfGfGRFRe/g/ZFcIASR/DXyxHKgtB1x8gZ98x+m4/UyhvgZnrpfgzPVbOJRahIzCtkOHu5MdBgeqMTBAjfGRXoj1c4PGxZ5jOoiIuoBBpLPqqoEvfwFc+kx6HDoRmLcVcNbIWizqHiGkzd7OXC/BvisFyC2Tfr49d6gd7TAo0A0xvm4YHeaB+CB3+HMgKRFRtzGIdEZ9DbD9cWl2jNJOmpI7dS1gy51G+xp9dR1OXyvBqWu3cPGmNH22pLKu1Xkxvq4YHuKBcRFeGBHiwdBBRGQmDCJ3os8DPp4P5JyTxoMs2C5tTEd9QlVtAy7nlOGT09k4m1WKq4XlrVo7bJUKxPi5YkyYF4YGuyPOzxVRvq7yFJiIyMowiNzJN/9PCiEO7sCj/2YIsXBCCCTn6XEuqxQfnbyOpBxdq+AR4uWEMWGexm3tY/xcobLlAmFERHJgEOnIgdekwakKG+An3wK+A+UuEd1GCIG0gnIkZpXibFYJEq4W43pxpck5Ghd7jI3Q4JFhgYjzd4OfmtvaExFZCgaR9lxPAA78Ufp50vMMIRbCYBBIytXhSq4OV4sq8MOVAqTk603OUdkqMTzYAyNCPPD4mGAEujvKVFoiIroTBpG2CAF8tUr6WXsPMGWNvOWxcgW6aiRcLcbxq8U4kFKI3LJqk+ftbZQYHuKOIVp3jAj2wPhIDZy5dgcRUZ/A/63bkn0SKEqRumQefR/gbIleVVxeg+NXbyHhahESMoqRUVhh8ryLyhZDtGqEaZwRH+iOGQP9oHbixoJERH0Rg0hbTv9Luh/yOODqJ29ZrIC+ug6HUouw70o+Tl2/hexbVSbPKxTAwAA3jA33wtgIL4yL0HD3WSKifoJB5HaJHwMXtks/D18sb1n6KYNBmtmy41QWDqcV4VpxRauZLbF+rrinMXjcE+bFFg8ion7KrEFk/fr1+Oabb5CYmAh7e3uUlpaa83Ldl3se+OoZ6efB8wDtaHnL009U1zXgwo0ynL5+C6evleDM9RKUVZkuIhauccZ9cT6YHO2DAQFu8HTmYnFERNbArEGktrYW8+bNw9ixY/HPf/7TnJfqvmodsG0e0FALBI4AHn6XY0O64UZJJb44exMHUgtx8UYZahsMJs872tlgYpQGj43SYlCgGr5unFJLRGSNzBpEfvvb3wIAtm7das7L9IxDfwLK8wEXP+CxDwElt2nvrKb9Wo5lFONgaiHOXi9pNbNF46LCqFBpSu2oUE8MCHCDnQ3rmIjI2lnUGJGamhrU1NQYH+t0ut65cGk2cHyT9PNDbwNuAb1z3T6spKIWRzOK8M2FXJy5XoICfY3J80oFMDrME3OHBuKecC+EeDlxrxYiImrFooLIhg0bjK0overI/wGGeiBkPBA9vfev30dkFVfiqws5+Op8DpLzTBcRs1UqMDBQjcnR3hgX4YXBgWqu5UFERHfU5W+KdevW3TEsnDp1CiNHjuxyYdasWYPVq1cbH+t0Omi12i6/T5eUZjdP1x35pHmv1YcIIXAlV4+TmcU41TjANE9n2t0S6eOCqTHemDbAD4MD1XC055RaIiLqmi4HkZUrV2L+/PkdnhMaGnpXhVGpVFCpVHf12ruWex6AANyDgUH/07vXtkD5umocTCnEv45mtmr1UCiA8REazI73x31xPvBx5QBTIiLqni4HEY1GA41GY46yyCP/snQfMsEqZ8lU1tZjT1I+vr+Uh8TsUpNBpo52NhgV5onRoc0DTF0duJ4HERH1HLN24mdlZeHWrVvIyspCQ0MDEhMTAQCRkZFwcXEx56U7L/+SdG9Fm9qVVNTiYGohvr6QgwMphahvsZqYUgHE+LnhoSEBWDA6mAuJERGRWZk1iLz88st4//33jY+HDRsGANi/fz+mTJlizkt3TkMdcOUr6ed+HETqGww4kXkLe6/k42BKIa4Wme7dEujuiIeHBWJilAYDA9Vw4SBTIiLqJWb9xtm6datlryFyeguAxtYAv8GyFsUcLt0sw1cXcvDNhVzcKDHdvyXG1xX3D/DB3KGBiPJ1lamERERk7az7V9+889J96ETAuX+Me7lZWoUfruTji3M3cTar1Hjcw8kO98f54r44X9wT7gl3Jy6hTkRE8rPuIFKULt2PWCprMbqrqLwG/72ch09P38D5G6UQjY08SgUwY6Cf8cbptUREZGmsO4jcypDuvSLlLcddEELg1LUS/OtIJnYn5ZnsXjs82B1TYnwwf7SWU2yJiMiiWW8QaagHKoqkn90C5S1LF+iq67DjZDZ2Jt7E5ZzmJfBjfF3x6CgtZg/2h5+a4YOIiPoG6w0ilUUABKBQAk6ecpfmji7dLMM7B9KxJykfdQ1S84fKVolHhgdi6bgwxPhxwCkREfU91htEygukeycNoLTMsRPF5TU4mFqIfydcR2J2qfF4lI8LFo0NwYPxAfB05qBTIiLqu6w3iFQ0BhEXH3nL0YYruTr8O+EaPjtzE7UNBuPxe8I98fz0GIwI8eBOtkRE1C9YbxC5lSndW1AQScvXY+O+NHx3Mdc4+DTO3w1TY7zx+OhgaD2d5C0gERFRD7PeIJL8tXQfNlnecgDIKa3C3w9m4KMTWcbl1idGabBiaiTGhHmy9YOIiPot6w0i+jzpPmCYbEW4dLMM7x2+iu8u5hm7YGJ8XbHuoYG4J5wBhIiI+j/rDSKiceyFDANVs29V4s09qfji3E3jscGBaiwcE4z/GREEOxtlr5eJiIhIDgwiit750jcYBBJvlOKd/RnYeyXfePyBQX74+ZQIDA5UswWEiIisDoOIwvwtIkk5Orz0+QVcuFFmPDYixAOr7ovCxCgNAwgREVktBhEztYjUNRjwZWIO9iTl4b+XpRYQe1slZg/2x8p7IxHh7WKW6xIREfUlDCI9HESEEDiaXoyNe1Nx+nqJ8fjgQDXeenwYwjTOPXo9IiKivsyKg0jjQh092C2Sr6vGzz44Y1wF1d5WiWUTwzB9gB8GB6qhVLILhoiIqCUrDiI91yJSUlGLj05m4c//TQEAONrZ4NGRQXhqYjgXISMiIuoAg0g3g8jfD2bgzT2pqK2X3i9c44z3nxzNAEJERNQJDCJ3GUSEEHj7h3S8uScVABDr54qF94Rgwehg2LALhoiIqFMYRO4iiAgh8Mz2RHx1PgcAsGJqBJ6fHsNpuERERF3EINLFIFLfYMBbP6QbQ8gLM2LwiykRDCFERER3gUGkC0GktLIWL312Ed9flvapmRztjRVTI81ROiIiIqvAINLJIHL8ajFW70hETlk1bJUKPDpKyxBCRETUTdYbRAxNQeTOXSrHrxZjwXvHYRBAqJcT/vjIYIyL0Ji5gERERP2f9QaRTraIFJfX4KXPLsAggPtiffDmY0OhdrTrhQISERH1f9a733wng8jqT87jWnElvF1VeOPRIQwhREREPYhBpIMgcuraLRxMLYStUoFtT42Bu5N9LxWOiIjIOjCItBNEhBD47VeXAQD/MzwI0b6uvVUyIiIiq8Eg0k4Q2XulAJdu6qCyVeLFB2J7sWBERETWg0FEadPqqQMpBVj279MAgCXjQuHpzC4ZIiIic7DOICIEACH9fFuLiK66Dmu/uARA2j/m/02P7uXCERERWQ8rDiKNWgQRg0HgmY/P4WZpFTyc7PCfn4+DyrZ1iwkRERH1DCsNIobmn1ssaPavo5k4kFIIla0S/35yDFxU1rvMChERUW9gEGlsEUnJ0+PNPakAgJd/NACDg9RylIyIiMiqMIg0BpGXv7yEytoGjAr1wOOjgmUqGBERkXVhEFEoUV3XgHNZpQCA388dBKXyzvvPEBERUfcxiCiU+PD4ddQ2GODjqkIMFy4jIiLqNQwiCiV2J+UDAJaOD4WiE7vxEhERUc+w+iAioEBKnh4AMDnaW64SERERWSWrDyKF5XUoq6qDUgFEeLvIWCgiIiLrY6VBpHlBs+sl1QCAQA9HONhx8TIiIqLeZKVBpLlFJLukCgCg9XCSqzRERERWy0qDSIN0r1Ai61YlAAYRIiIiOVhpEGlsEVEojeuHRPpwfAgREVFvs+ogIhRKHL9aDACYHMMZM0RERL3NqoOIAQrU1BsQoHZAFFtEiIiIep1VB5EGIX38CVEaLmRGREQkAysPIlL4iOay7kRERLKw0iAirSPSFEQi2C1DREQkCysNIlKLSH3jumYcH0JERCQPqw4iBiigcbFHoLujzAUiIiKyTlYeRJSI83fjQFUiIiKZWHUQEVDAkfvLEBERycaqg4gBCtjZWmcVEBERWQLr/BZu0TVjb2OdVUBERGQJzPYtfO3aNfzv//4vwsLC4OjoiIiICLzyyiuora011yU7r2WLiA3HhxAREcnF1lxvnJycDIPBgM2bNyMyMhKXLl3CsmXLUFFRgddff91cl+2cFmNEbNkiQkREJBuzBZGZM2di5syZxsfh4eFISUnBpk2bLCCISAuIGISCXTNEREQyMlsQaUtZWRk8PT3bfb6mpgY1NTXGxzqdzjwFaTFGhF0zRERE8um15oCMjAy8/fbbWL58ebvnbNiwAWq12njTarXmKYzJGBG2iBAREcmly9/C69atg0Kh6PB2+vRpk9fk5ORg5syZmDdvHp566ql233vNmjUoKysz3rKzs7v+iTqjxRgRBhEiIiL5dLlrZuXKlZg/f36H54SGhhp/zsnJwdSpUzF27Fi8++67Hb5OpVJBpVJ1tUhd13L6LtcRISIikk2Xg4hGo4FGo+nUuTdv3sTUqVMxYsQIbNmyBUqlhXzpGxqkO07fJSIikpXZBqvm5ORgypQpCA4Oxuuvv47CwkLjc35+fua6bOdwjAgREZFFMFsQ2b17N9LT05Geno6goCCT50Tj9FnZmMyaYRAhIiKSi9m+hZcuXQohRJs32dnYoUzpgRLhwnVEiIiIZGSd38Jhk7AycAeeqFsLW44RISIiko11BhEAtfVS9wy7ZoiIiORjtd/CdQ0MIkRERHKz2m/heoM0VsXell0zREREcrHaINLUNWNrKWubEBERWSGr/RZuaGwR4WBVIiIi+VhtEDE0TiNWKhhEiIiI5GK1QaRpORMGESIiIvlYbRBpbhGRuSBERERWzIqDiHSvYIsIERGRbKw4iLBFhIiISG5WG0Q4RoSIiEh+VhxEOGuGiIhIblYbRJrHiMhbDiIiImtmxUFESiIMIkRERPKx4iAi3bNrhoiISD5WG0Q4RoSIiEh+VhtEOH2XiIhIflYcRKR7LmhGREQkHysOImwRISIikpvVBhEuaEZERCQ/qw0iBg5WJSIikp3VBxHmECIiIvlYcRCR7pUcJEJERCQbqw0igoNViYiIZGe1QYQrqxIREcnPioMIx4gQERHJzSqDiBCC03eJiIgsgJUGkeafGUSIiIjkY5VBxNAiiXCwKhERkXysNIg0/8y9ZoiIiORjpUGELSJERESWwCqDCMeIEBERWQarDCKmLSIMIkRERHKx+iDCHEJERCQfKw0izT+zRYSIiEg+VhlEBAerEhERWQSrDCJsESEiIrIMVhlEBMeIEBERWQSrDCJNLSIKBRc0IyIikpNVBpGmFhF2yxAREcnLKoOIwbjzrrzlICIisnZWGkSkJMJuGSIiInlZdRBhiwgREZG8rDKICGPXDJMIERGRnKwyiBg4WJWIiMgiWGkQke6ZQ4iIiORlpUGELSJERESWwCqDiOBgVSIiIotglUHEwMGqREREFsFKgwjXESEiIrIE1hlEDNI9u2aIiIjkZZ1BhINViYiILIJVBhHBvWaIiIgsglUGEY4RISIisgxmDSIPPfQQgoOD4eDgAH9/fyxatAg5OTnmvGSnGLtmrDKGERERWQ6zfhVPnToVn3zyCVJSUvDZZ58hIyMDP/7xj815yU7h9F0iIiLLYGvON3/uueeMP4eEhOCll17C3LlzUVdXBzs7O3NeukOCg1WJiIgsglmDSEu3bt3Ctm3bMG7cuHZDSE1NDWpqaoyPdTqdWcrCvWaIiIgsg9lHSbz44otwdnaGl5cXsrKy8OWXX7Z77oYNG6BWq403rVZrljJx+i4REZFl6HIQWbduHRQKRYe306dPG89/4YUXcO7cOezevRs2NjZYvHixsWvkdmvWrEFZWZnxlp2dffefrAMG7jVDRERkEbrcNbNy5UrMnz+/w3NCQ0ONP2s0Gmg0GkRHRyMuLg5arRbHjx/H2LFjW71OpVJBpVJ1tUhd1pSDFGASISIiklOXg0hTsLgbTS0hLceByKF5HRFZi0FERGT1zDZY9eTJkzh58iQmTJgADw8PXL16FS+//DIiIiLabA3pTZy+S0REZBnMFkQcHR3x+eef45VXXkFFRQX8/f0xc+ZMbN++vVe6XzoS7OmElVMj4e0qbzmIiIisnUK0N3LUAuh0OqjVapSVlcHNzU3u4hAREVEndOX7m4ucExERkWwYRIiIiEg2DCJEREQkGwYRIiIikg2DCBEREcmGQYSIiIhkwyBCREREsmEQISIiItkwiBAREZFsGESIiIhINgwiREREJBsGESIiIpINgwgRERHJxlbuAnSkaWNgnU4nc0mIiIios5q+t5u+xzti0UFEr9cDALRarcwlISIioq7S6/VQq9UdnqMQnYkrMjEYDMjJyYGrqysUCkWPvrdOp4NWq0V2djbc3Nx69L2pGeu5d7Ceew/runewnnuHuepZCAG9Xo+AgAAolR2PArHoFhGlUomgoCCzXsPNzY1/yXsB67l3sJ57D+u6d7Cee4c56vlOLSFNOFiViIiIZMMgQkRERLKx2iCiUqnwyiuvQKVSyV2Ufo313DtYz72Hdd07WM+9wxLq2aIHqxIREVH/ZrUtIkRERCQ/BhEiIiKSDYMIERERyYZBhIiIiGRjlUHknXfeQVhYGBwcHDBixAgcPnxY7iL1KRs2bMCoUaPg6uoKHx8fzJ07FykpKSbnCCGwbt06BAQEwNHREVOmTMHly5dNzqmpqcHTTz8NjUYDZ2dnPPTQQ7hx40ZvfpQ+ZcOGDVAoFHj22WeNx1jPPePmzZt44okn4OXlBScnJwwdOhRnzpwxPs967hn19fX49a9/jbCwMDg6OiI8PBy/+93vYDAYjOewrrvu0KFD+NGPfoSAgAAoFArs3LnT5PmeqtOSkhIsWrQIarUaarUaixYtQmlpafc/gLAy27dvF3Z2duK9994TSUlJYtWqVcLZ2Vlcv35d7qL1GTNmzBBbtmwRly5dEomJiWL27NkiODhYlJeXG8959dVXhaurq/jss8/ExYsXxWOPPSb8/f2FTqcznrN8+XIRGBgo9uzZI86ePSumTp0qhgwZIurr6+X4WBbt5MmTIjQ0VMTHx4tVq1YZj7Oeu+/WrVsiJCRELF26VJw4cUJkZmaKvXv3ivT0dOM5rOee8Yc//EF4eXmJr7/+WmRmZopPP/1UuLi4iI0bNxrPYV133bfffivWrl0rPvvsMwFAfPHFFybP91Sdzpw5UwwaNEgcO3ZMHDt2TAwaNEg8+OCD3S6/1QWR0aNHi+XLl5sci42NFS+99JJMJer7CgoKBABx8OBBIYQQBoNB+Pn5iVdffdV4TnV1tVCr1eLvf/+7EEKI0tJSYWdnJ7Zv32485+bNm0KpVIrvv/++dz+AhdPr9SIqKkrs2bNHTJ482RhEWM8948UXXxQTJkxo93nWc8+ZPXu2ePLJJ02OPfLII+KJJ54QQrCue8LtQaSn6jQpKUkAEMePHzeek5CQIACI5OTkbpXZqrpmamtrcebMGUyfPt3k+PTp03Hs2DGZStX3lZWVAQA8PT0BAJmZmcjLyzOpZ5VKhcmTJxvr+cyZM6irqzM5JyAgAIMGDeKfxW1WrFiB2bNn4/777zc5znruGbt27cLIkSMxb948+Pj4YNiwYXjvvfeMz7Oee86ECROwb98+pKamAgDOnz+PI0eOYNasWQBY1+bQU3WakJAAtVqNMWPGGM+55557oFaru13vFr3pXU8rKipCQ0MDfH19TY77+voiLy9PplL1bUIIrF69GhMmTMCgQYMAwFiXbdXz9evXjefY29vDw8Oj1Tn8s2i2fft2nD17FqdOnWr1HOu5Z1y9ehWbNm3C6tWr8atf/QonT57EM888A5VKhcWLF7Oee9CLL76IsrIyxMbGwsbGBg0NDVi/fj0ef/xxAPw7bQ49Vad5eXnw8fFp9f4+Pj7drnerCiJNFAqFyWMhRKtj1DkrV67EhQsXcOTIkVbP3U0988+iWXZ2NlatWoXdu3fDwcGh3fNYz91jMBgwcuRI/PGPfwQADBs2DJcvX8amTZuwePFi43ms5+7bsWMHPvzwQ3z00UcYOHAgEhMT8eyzzyIgIABLliwxnse67nk9Uadtnd8T9W5VXTMajQY2Njat0ltBQUGrtEh39vTTT2PXrl3Yv38/goKCjMf9/PwAoMN69vPzQ21tLUpKSto9x9qdOXMGBQUFGDFiBGxtbWFra4uDBw/irbfegq2trbGeWM/d4+/vjwEDBpgci4uLQ1ZWFgD+fe5JL7zwAl566SXMnz8fgwcPxqJFi/Dcc89hw4YNAFjX5tBTdern54f8/PxW719YWNjtereqIGJvb48RI0Zgz549Jsf37NmDcePGyVSqvkcIgZUrV+Lzzz/HDz/8gLCwMJPnw8LC4OfnZ1LPtbW1OHjwoLGeR4wYATs7O5NzcnNzcenSJf5ZNLrvvvtw8eJFJCYmGm8jR47EwoULkZiYiPDwcNZzDxg/fnyr6eepqakICQkBwL/PPamyshJKpenXjo2NjXH6Luu65/VUnY4dOxZlZWU4efKk8ZwTJ06grKys+/XeraGufVDT9N1//vOfIikpSTz77LPC2dlZXLt2Te6i9Rk///nPhVqtFgcOHBC5ubnGW2VlpfGcV199VajVavH555+Lixcviscff7zN6WJBQUFi79694uzZs+Lee++16il4ndFy1owQrOeecPLkSWFrayvWr18v0tLSxLZt24STk5P48MMPjeewnnvGkiVLRGBgoHH67ueffy40Go345S9/aTyHdd11er1enDt3Tpw7d04AEG+++aY4d+6ccVmKnqrTmTNnivj4eJGQkCASEhLE4MGDOX33bv3tb38TISEhwt7eXgwfPtw47ZQ6B0Cbty1bthjPMRgM4pVXXhF+fn5CpVKJSZMmiYsXL5q8T1VVlVi5cqXw9PQUjo6O4sEHHxRZWVm9/Gn6ltuDCOu5Z3z11Vdi0KBBQqVSidjYWPHuu++aPM967hk6nU6sWrVKBAcHCwcHBxEeHi7Wrl0rampqjOewrrtu//79bf6fvGTJEiFEz9VpcXGxWLhwoXB1dRWurq5i4cKFoqSkpNvlVwghRPfaVIiIiIjujlWNESEiIiLLwiBCREREsmEQISIiItkwiBAREZFsGESIiIhINgwiREREJBsGESIiIpINgwgRERHJhkGEiIiIZMMgQkRERLJhECEiIiLZMIgQERGRbP4/BviVp84DkI4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log10(sigma_stan), label=\"log10 eigenvalue stan\")\n",
    "plt.plot(np.log10(sigma_nutpie), label=\"log10 eigenvalue nutpie\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0204068e-b465-48c6-a7bd-1d318cca92ef",
   "metadata": {},
   "source": [
    "## Estimate mass matrix with a finite sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3b9b8242-aef6-4ab6-9c6d-28e3ad87c432",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_draws = 200\n",
    "draws = stats.multivariate_normal(np.zeros(N), cov=Sigma).rvs(num_draws)\n",
    "grads = - linalg.solve(Sigma, draws.T, assume_a=\"pos\").T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3ced6aac-727c-40dd-81ca-57e30d781d43",
   "metadata": {},
   "outputs": [],
   "source": [
    "estimate_no_decay = np.sqrt(draws.var(axis=0) / grads.var(axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "caf47f45-5b92-49af-a5a6-96a0657f5dc7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fad7f4c4830>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKbRJREFUeJzt3X901PWd7/HXJJIMP5LBIQ0zlABZwF1iFA2KRqoLFLhBL4rdYr17cdWr3oUSbpVz91a6pxvoLpvd1V13j9QUeyz2Xg62e3rqD6qbs7FasCo3i0g1pmrBIFyTlB+RmRhMAjPf+0c6KUlmwiSZ73y+M9/n4xzOaYav+b570javfj7vz/vjsSzLEgAAgAE5pgsAAADuRRABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYMwlpgsYTjQaVWtrqwoKCuTxeEyXAwAAkmBZljo7OzVt2jTl5Ay/5uHoINLa2qqSkhLTZQAAgFE4fvy4pk+fPuwzjg4iBQUFkvr+jRQWFhquBgAAJCMcDqukpKT/9/hwHB1EYtsxhYWFBBEAADJMMm0VNKsCAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjHH0QDMAAGCPSNRSY0uHTnR2q7jAq4WlfuXmpP9eN4IIAAAuU9/Upq17mtUW6u7/LOjzqmZVmarKg2mtha0ZAABcpL6pTet3HRwQQiSpPdSt9bsOqr6pLa31EEQAAHCJSNTS1j3NsuL8XeyzrXuaFYnGe8IeBBEAAFyisaVjyErIhSxJbaFuNbZ0pK0mgggAAC5xojNxCBnNc6lAEAEAwCWKC7wpfS4VCCIAALjEwlK/gj6vEh3S9ajv9MzCUn/aaiKIAADgErk5HtWsKpOkIWEk9nXNqrK0zhMhiAAA4CJV5UHVra1QwDdw+yXg86pubUXa54gw0AwAAJepKg9qeVmAyaoAAMCM3ByPKmdPMV0GWzMAAMAcgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGFuDSG1tra699loVFBSouLhYq1ev1gcffGDnKwEAQAaxNYjs3btXGzZs0P79+9XQ0KDz589rxYoV6urqsvO1AAAgQ3gsy7LS9bKTJ0+quLhYe/fu1U033XTR58PhsHw+n0KhkAoLC9NQIQAAGKuR/P5Oa49IKBSSJPn9/nS+FgAAOFTabt+1LEubNm3Sl770JZWXl8d9pqenRz09Pf1fh8PhdJUHAAAMSNuKSHV1td555x0988wzCZ+pra2Vz+fr/1NSUpKu8gAAgAFp6RHZuHGjnnvuOe3bt0+lpaUJn4u3IlJSUkKPCAAAGWQkPSK2bs1YlqWNGzfq2Wef1S9+8YthQ4gk5efnKz8/386SAACAg9gaRDZs2KDdu3fr+eefV0FBgdrb2yVJPp9P48ePt/PVAAAgA9i6NePxeOJ+vnPnTt1zzz0X/ec5vgsAQOZx1NYMAABAItw1AwAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjLjFdAAAAIxGJWmps6dCJzm4VF3i1sNSv3ByP6bIwSgQRAEDGqG9q09Y9zWoLdfd/FvR5VbOqTFXlQYOVYbTYmgEAZIT6pjat33VwQAiRpPZQt9bvOqj6pjZDlWEsCCIAAMeLRC1t3dMsK87fxT7buqdZkWi8J+BkBBEAgOM1tnQMWQm5kCWpLdStxpaO9BWFlCCIAAAc70Rn4hAymufgHAQRAIDjFRd4U/ocnIMgAgBwvIWlfgV9XiU6pOtR3+mZhaX+dJaFFCCIAAAcLzfHo5pVZZI0JIzEvq5ZVcY8kQxEEAEAZISq8qDq1lYo4Bu4/RLweVW3toI5IhmKgWYAgIxRVR7U8rIAk1WzCEEEAGDcSMa25+Z4VDl7SporhF0IIgAAo0Y6tp27ZrILQQQAYExsbPvgeaixse2Dez+4ayb70KwKADBipGPbuWsmOxFEAABGJDu2ff+R09w1k8UIIgAAIxqa25N6bsPug9r+ymHumslSBBEAQNrVN7XpB68fTerZM5+f02Mvf5jUs9w1k3kIIgCAtIpts9iBu2YyD6dmAABJiR2bbQ99ro6uXvkn5StQOPLjsxfrDRkNj/omrHLXTOYhiAAALiresdkY/8Q8/c1t5br5yuSOz451+8QjDWha5a6ZzMbWDABgWImOzcZ0dPXq67sPqval5LZbxrJ98tCyudw1k2VYEQEAJDTcsdnBduxr0fzpk3XzldOGfW5hqV9Bn3dE2zOxrZfqpXNVvXQuk1WzCEEEAJDQSPs5/uIn7+iTM5+raFK+Ar7xcUNCbo5HNavKtG7XwRHVcuHWC3fNZA+CCAAgoZH2c3T1RrTtpff7v/ZPHKfbr/qilpUFBoSSqvKgnvjTClU/c1AXm0HGCPfsRhABACQ01uOwHV3n9NTrR/XU60eHBIqbrwxqu67W13e/nfCff2hZ31YMWy/Zi2ZVAEBCsX6OVIh3J8zNV07T99ZWDHlH0OfV99ZW6BvLLiOEZDmPZVmOHcwfDofl8/kUCoVUWFhouhwAcKVEN+SORqzp9JffXDogYMRmlNCAmh1G8vubFREAwLCWlwX04LK5mjx+3Ji/V6I7YXJzPKqcPUW3XfVFVc6eQghxEXpEAAAJxRtkNn5cjixJ3eeio/6+3AmDGFZEAABxJRpk1n0uqu5zUVVdPlUT83NH9b2LJuanokRkAVZEAABDDDfILPZZ/Xu/Hf0L2HnB79i6IrJv3z6tWrVK06ZNk8fj0XPPPWfn6wAAKbL/yOmUX0x3oVd+PYYQg6xiaxDp6urS/PnztX37djtfAwBIgUjU0ptHTus7e97Tn+86YOu7nj30iSIXm2QGV7B1a2blypVauXKlna8AAKTAcLfr2qGj65waWzoY1Q5n9Yj09PSop6en/+twOGywGgDIfpGope2v/EaPvfybtL+bkzOQHBZEamtrtXXrVtNlAIAr1De1acsLzWoPmwkEYx0fj+zgqOO7mzdvVigU6v9z/Phx0yUBQFaKHc01EUI86hvhvrDUn/Z3w3kctSKSn5+v/HzOlgOAnSJRS1teeC8lI9tHKnZqt2ZVGdNTIclhQQQAYL/trxxWe7jn4g/aIDDoBl7A1iDy2Wef6fDhw/1ft7S06NChQ/L7/ZoxY4adrwYAxFHf1KbHXv4wLe8K+rz69i3zdOnEfC6zQ0K2BpEDBw5oyZIl/V9v2rRJknT33Xfr6aeftvPVAIBBYtNS7XbfollaVhYgdCAptgaRxYsXy7IYWAMATtDY0pHSOSET83LV1Rvp/zrItgtGgR4RAMgikailxpaOAVshkrT/o9P64RtHU/quv7n9CgUKvWy7YEwIIgCQBWKDyXa+flRnPj/X//mEvFydj0TVG0n96nSg0MtkVIwZQQQAMtCFKx9HT3XpB6+3KPT5+SHPnb1g6yRVPOo7/cIcEKQCQQQAMky674W5EHNAkGoEEQDIILGJqKaOATAHBKlGEAEAB0rUdLp1T7OxECJJ375lHiEEKUUQAQCHibf14p84TteX+o1sx8R4JP31i7/WfyoPsi2DlHHUpXcA4HaxrZfBgaOj65xeavqtoar6WJLaQt1qbOkwWgeyC0EEABwiNvnU6WMgT3SaW5VB9iGIAIBDpHryqV2KC7ymS0AWIYgAgEM4YaXBP3GcEnV/eNQ3xp35IUglgggAOITplYagz6u/ua1ckoaEEeaHwC4EEQBwiIWlfgV93oQrEnbyqC9k3HzlNNWtrVDANzAUBXxe1a2t4OguUs5jOfh63HA4LJ/Pp1AopMLCQtPlAIDt6pvatG7XwbS+M96tufHmmLASgmSN5Pc3c0QAwJALf9kXTcxX1LLU9ElY86f79Kv/F7L9/ffeMFMrLg/GDRm5OR4utENaEEQAwID6pjZteaFZ7eH0N6gGCvO15dbL2WaBIxBEAMBmg1c+/uNoh/75579Jaw3fvmWeigry2WaB4xBEAMBGJm/KlfqaUAM+r+5ZVEr4gCMRRADAJqZvyo3hyC2cjCACADaIRC09/NN3jYaQeKdhAKchiACADf6l4UOdOXvOyLurLp+qu28opRcEGYEgAgAptu3FZn3/tRYj7w76vPruf11AAEHGIIgAQApte/E9ff+1o8beTz8IMg0j3gEgRfb8qtVoCLlv0Sz6QZBxWBEB4HqpGGde39Smjc+8bVOFyVlWFjD6fmA0CCIAXC3enI+LnTaJBZf20Ofq6OrV5Al52vbSr9NV8hCxWSELS/3GagBGiyACwLUSzfloD3Vr/a6DcW+bNT2gbLDYug29IchU9IgAcKVI1NLWPc1x53zEPtu6p1mR6O+fiAUXp4QQqW8lJF5gAjIFKyIAXKmxpWPYQGFJagt1q7GlQ5WzpwwbXNIl6PPq27fM06UT88fUzwI4CUEEgCud6ExuVSP23MWCi52ql8zWojlfIHQgKxFEALhScYE3qedOhHv02ocn9b/fPGpvQXHEmlAfWv6HBBBkLYIIAFdaWOpX0OdVe6h72O0WU6dhaEKFW9CsCsCVcnM8qllVJun3v/SdhCZUuAUrIgBcq6o8qLq1FY46jnvvDTO14vIg/SBwDYIIAFerKg9qeVlAP/hli9GhZJL0xJ9W6OYrWQGBuxBEALheQ3O7vvvqYWPvv3TCJar9ypVsw8CVCCIAXCsStbT9ld/osZd/Y6yGKRPz9ObmLyvvElr24E4EEQCuMPhiu0+7evWdnzWrPWyuN8Qjadvt5YQQuBpBBEDWc9r9MJKU45G2/5er2Y6B6xFEAGSFwSsesVMniS62My1qSZdOzDddBmAcQQRAxou34hEo9OrOa0v09BtHHRdCYpIdMw9kM4IIgIyWaMWjPdytf/65uSbUZCQ7Zh7IZgQRABnLCTfijkbsDpmFpX7TpQDG0aoNIGOZvBF3tLhDBhiIFREAGSsTeywCPq9qVpVxWgb4HYIIgIyVST0Wf1Y5UyvLuUMGGIwgAiBjfdrVK4+UET0iK8uDqpw9xXQZgOMQRABkHCeMZk8WjanA8AgiADJCbGBZQ3O7njvUqo6uXtMlDTF4dYbGVODiCCIAHM+JI9oHe2jZZfrRfxwbOFSNxlTgoggiABzNqSPaY2JbL9VL56h66Zy4Y+YBJEYQAeBYn/dG9Bc/+ZWjQ4g0cOuFhlRgZAgiAByp9qVmPflaiyynphCx9QKkAkEEgKNEopYe/NFB7Xmn3XQpA8S2YB796nyd6uph6wVIEYIIAEfoO5J7WD94vUWhz8+ZLmeAC7dgFs0tMloLkG0IIgCMq29q08M/fVdnzjorgMSwBQPYhyACIK1i80BiJ0s+7erV13cfNF3WEA8tm6tZRRPZggFsRhABkDbx5oF4HPj7/d4bZuobyy4zXQbgCgQRAGmRaB6IE0/FrLicLRggXXJMFwAg+0WilrbuaXbsPJALBbkXBkgrgggA2zW2dDh6PHuMR9wLA6QbWzMAbHei0/khJMjJGMCItKyIPPHEEyotLZXX69WCBQv02muvpeO1AByiuMBruoS4vlrxRT12x3w988D1+uU3lxJCAANsXxH58Y9/rAcffFBPPPGEFi1apB07dmjlypVqbm7WjBkz7H49gDQYfCR38HHXq0ommysujskTxunvvnIFwQNwAI9l2duzft1116miokJ1dXX9n82bN0+rV69WbW3tsP9sOByWz+dTKBRSYWGhnWUCuIhEYSPekdygz6tv3TxPJ8Ld+rjjrNrOfK6GX58wWP3vfbXii/r7r86nDwSw0Uh+f9u6ItLb26u33npLDz/88IDPV6xYoTfeeGPI8z09Perp6en/OhwO21kegCQlChu3zg/qyX0tQ07DtIW6tfGZt9NbZBIChfmEEMBhbO0ROXXqlCKRiKZOnTrg86lTp6q9feiFVrW1tfL5fP1/SkpK7CwPQBJi8z8Gn3ppC3VrR5wQ4mRbbr2cEAI4TFqaVT2DRidaljXkM0navHmzQqFQ/5/jx4+nozwACWTS/I+LeWjZXHpCAAeydWumqKhIubm5Q1Y/Tpw4MWSVRJLy8/OVn59vZ0kARiBT5n9czOTx41S9dK7pMgDEYeuKSF5enhYsWKCGhoYBnzc0NOiGG26w89UAUiAT5n8k495FpWzJAA5l+/HdTZs26a677tI111yjyspKPfnkkzp27JjWrVtn96sBjJFT53+MxOQJ41S9dI7pMgAkYHsQ+drXvqbTp0/rO9/5jtra2lReXq6XXnpJM2fOtPvVAMZoYalfQZ9X7aHujO0T+buvXMFqCOBgts8RGQvmiADmxU7NSHJsGLnlioDe+vhTtYd/f/yfke2AOY6ZIwIg8y0vC+jBZXO1Y99HOtsbMV3OEA8tm6tvLLvsotNdATgTQQRAQvEGmTlJ0OftPw2Tm+NR5ewphisCMFIEEQCSho5w/7SrVxt2H3TkdkxsnaNmVRmrHkCGI4gAiLvykeMx3xNy7w0z1dl9Xg2/PqHQ5+f6Pw/Q/wFkDYII4FKxFZCXm9v11OtHh/x91HQKkVT/3m/1y28u1d9L9H8AWYogAriQ03s/YtpC3Wps6VDl7Cn0fwBZiiACuEzsOK4DFjySki3TXQHEl5ZL7wA4QyRqacsL72VMCJGyY7orgMRYEQFcZPsrhwcM/XIyj/qaUheW+k2XAsBGrIgALlHf1KbHXv7QdBlJs8TxXMANCCKAC0SilrbuaTZdxoj8t0WzOJ4LuABBBHCBxpYOx5+QGWx5WcB0CQDSgCACuICTTp585eppChTmK9GGi0d9o9vpDQHcgWZVIIvEhpS1hz5XR1ev/JPyVTwpXycc0qDqG3+JHllzlRqa27V+10F5NHB6K6PbAfchiABZIhOGlP39n1yp3ByPqsqDqltbMaReRrcD7kMQAbKA04eUTcjL0T/dcdWAgFFVHtTysgCj2wGXI4gAGejCm3KLJuU7dkiZR9ItVwb1L3deHTdg5OZ4GN0OuBxBBMgwmbAFc+PcIi2+7Au6q3KW8i6hJx5AYgQRIIM4fQsm5uuL57DSASAp/F8VIEPEhpI5PYRw9BbASBBEgAyRKUPJOHoLYCQIIkCGaGhuN13CsCZPGKfvra3g6C2AEaFHBHCoSNTS/iOn9eZHp3T45Geqb/qt6ZLimjx+nO5dVKrqpXNYCQEwYgQRwIHqm9r08E/f1Zmz50yXMqw/uXqa/mHNVQQQAKNGEAEcpr6pTet2HTRdxkVNzMslhAAYM4II4ACxAWWtZz5XzQtNpstJyj/eMZ8QAmDMCCKAYU4eUHbphHHKy/Xot529/Z8FCvO15dbLaUoFkBIEEcAgJw8o80iq/coV3AcDwFYEEcAQJw8oCw66BZcpqQDsQhABDHHigLJCb67q1l6j6/9gCqseANKCIAIYcqLTOSEkFjn+4avztWhOkdFaALgLQQQw5OipLtMl9AsM2ooBgHQhiAAG1De16bGXf2O6DEmSf+I47f2LJcq7hBsfAKQf/8sDpFkkamnLC++ZLqNfR9c5vfXxp6bLAOBSBBEgzRpbOtQe7jFdxgBO6lcB4C5szQA26z0f1f9586g+7jirmf4JKvQ67792xQVe0yUAcCnn/S8ikEVqX2rW919rUdSJw0LUd1om4PNqwcxL9eaR0wwtA5B2BBHAJrUvNWvHvhbTZSQUixm3zg/qjx95dcBMk8EDzQDALvSIADboPR/V919zRgjxSPrzm0oV9A3cfgn4vPrvN5XqyX0tQwartYe6tX7XQdU3taWxUgBuxIoIYINv/fRdR2zHTJmYp223l6uqPKj/VTVvwJ0xC2Zeqj9+5NW4I+Yt9QWYrXuatbwswDYNANsQRIAUiUQtNbZ0qKG5XT85+P9MlyP/xHF6c/OX++eD5OZ4BtwZ8+aR08OOmLcktYW61djSwV0zAGxDEAFSoL6pTVv3NDvm7hiPpL+9/Yphh5Qle2SXo70A7EQQAcbopXda9fXdb5suo59/4jj97e1XXLTRNNkjuxztBWAnmlWBMXjpnTZVP+OcECJJ3/7Plyd12mVhqV9Bn1eJuj886js9s7DUn9L6AOBCBBFglOqb2vT13Qcd0ZR6oUBhcisYuTke1awqk6QhYST2dc2qMhpVAdiKIAKMgtPui4kZ6QpGVXlQdWsrFIhztLdubQVzRADYjh4RYIQiUUvf/MmvHHdfjEejW8GoKg9qeVlgwNFeJqsCSBeCCDBI7BhuvF/K9U1tevin7+rM2XOGqxxorJNQBx/tBYB0IYgAF4h3DNc/MU+rr5om3/g8Pfbyhwar+72gz6s7r52hWUUTWMEAkNEIIsDv1De1af2ug0MmjXZ09eoHrx81UVK/QGG+/vGOq3Tqsx6CB4CsQhAB1Lcds3VPc9xx56Z5JG259XItmlNkuhQASDlOzQCSGls6jE9FnZiXqwl5uQM+C3J6BUCWY0UEkNTQ3G70/Q8tu0zVS+dIEqdXALgKQQSuF4la+tcDZi6pm5iXq3+8Y/6AFQ9OrwBwE4IIXG/7K4f1Wc95I+9+8q5rtGguvR8A3IsgAleLRC3t2Hck7e/1qG966fWsfgBwOZpV4WqP//xDne2NpPWd3OMCAL/HighcJxK1tP/Iae36v0f1b02/Tfv7A2OcggoA2YQgAlcxOaL99qum6Y5rZ3ASBgAuQBCBa9Q3tWndroPG3l8wfhwnYgBgEIIIsl4kamn/R6e16V8PGa1jpn+C0fcDgBPZGkS2bdumF198UYcOHVJeXp7OnDlj5+uAIeJdYmdCjke6q3KW0RoAwIlsDSK9vb1as2aNKisr9dRTT9n5KqBfbAVk1/6P9W9NZiemxjxwY6nyLuGQGgAMZmsQ2bp1qyTp6aeftvM1QD+Tzajx5Hj6Qsjmm8tMlwIAjkSPCLKG6WbUmMpSv+YGCjTTP0F3Vc5iJQQAhuGoINLT06Oenp7+r8PhsMFqkEkiUUtbXmg2XYYk6WjHWe164HqO6AJAEkb8f9W2bNkij8cz7J8DBw6Mqpja2lr5fL7+PyUlJaP6PnCPSNTSm0dO67GGD9QeNtuQGtMW6lZjS4fpMgAgI4x4RaS6ulp33nnnsM/MmjVrVMVs3rxZmzZt6v86HA4TRpCQU07ExHOi03k1AYATjTiIFBUVqajInttC8/PzlZ+fb8v3Rnapb2rT+l0HZZkuJIHiAq/pEgAgI9jaI3Ls2DF1dHTo2LFjikQiOnTokCRpzpw5mjRpkp2vRhaLRC1t3dNsLITkeKRogpfHbtVdWOpPa00AkKlsDSJ/9Vd/pR/+8If9X1999dWSpFdffVWLFy+289XIUr3no/rrn71nZDvmzypnamV5UJ929WjD7rclaUAY4lZdABg5j2VZTl3dVjgcls/nUygUUmFhoelyYFjtS836/mstCVcj7PbMA9f33xUTrz8lyK26ACBpZL+/HXV8F0hk24t9IcSU4KDtlqryoJaXBdTY0qETnd0qLvByqy4AjAJBBI73s0OtRkOIR/G3W3JzPNymCwBjRBCBo0Si1oBVhk+7elT9o7eN1cN2CwDYiyACx4jXd5HujY71i/9ARRPz5Z+Yp4BvPNstAGAzgggcIdFckHT2pQZ9Xv3PFX9E8ACANOI2Lhhnei5IDMduASD9CCIwrrGlw/iY9oeWzaUPBAAMYGsGaTe4IdX0ZXVBn1fVS+carQEA3IoggrSK15B66YRxRmphEioAmMfWDNIm1pA6eBvm07PnbH1vjkd64MZSBX0DL6IL+LyqW1vBlgwAGMSKCNLCZENq1JKW/tFUPbxyHpNQAcBhCCJIC9MNqSc6u5mECgAOxNYM0uJEp9mG1OIC78UfAgCkHSsisM2Fp2NOdfYYqcGjvl6QCy+sAwA4B0EEtkg0rt2OHpHJE8bpzNlzQ74/p2IAwPkIIki5dIxrf/DLc1X6hYn9TacNze1Dgk+AC+sAwPEIIkgpu0/HXDphnGq/csWQcFFVHtTysgCnYgAgwxBEkFJ2nY6ZPH6c7l00S9VL5yYMF5yKAYDMQxBBSjU0t6f0+335j76g+2+czeoGAGQpgghSpr6pTT94/WhKvtek/Ev0D39ypW6+kv4OAMhmBBGkRKw3JFV2rF2gRXOLUvb9AADORBBBUi6cCVI0MV/ySKc+6+lvCk1lb0igMF/X0+sBAK5AEMFFxZsJcqGgz6uiSXkpe9+WWy+nHwQAXIIggmElmglyobZQd8pWQx5adhlzPwDARbhrBgn1no/qW8++m7YbcwOF+apeOidNbwMAOAFBBHHVN7Xp+tqfq6PrXMq/9+BNF8/v/rAlAwDuQxDBELHtmI6u3pR/74eWXaaAb+BNuAGfV3VrK9iSAQAXokcEA9g5on3COI+ql85R9dI5jGIHAEgiiGAQu0a0S9IDN83pDxyMYgcASGzNYJATnfaEkAl5ufofX55ry/cGAGQugggGKC7wXvyhUfinO+az/QIAGIKtGQyZmhoo9Ko9nLopqVtuvZxGVABAXAQRl4s3NXVSfu6Yv+/qq6bpa9fOoBEVADAsgoiLJZqa+llPZMzf+4//sJiGVADARdEj4lJ2HtOVpEChPb0mAIDsQhBxKTuP6QZ9fbNBAAC4GLZmstiFTaiDB4fZcUw31glSs6qMvhAAQFIIIlkqXhNq0OdVzaoyVZUHbTmmG7jg+wMAkAyCSBZK1ITaHurW+l0HVbe2QsvLAgr6vGoPdaekT6R6yWw9tPwPWQkBAIwIPSJZZrgm1NhnW/c0S+rbQkmVRXO+QAgBAIwYQSSDRaKW3jxyWs8f+kRvHjnd3xMyXBOqJakt1K3Glg5VlQdVt7ZCgcL8UdfgEc2pAIDRY2smQ9U3tWnLC80DJqAGCr265YpAUv98rFm1qjyo5WUBfeNHb+tn77SNqAaaUwEAY0UQyUD1TW1at+vgkM/bw9166vWjSX2Pwc2qb3386YjroDkVADBWBJEME4laevin7w77jMcjyVLcPhGP+gLEhVspyc4U8U/M092VszSraMKQ48AAAIwGQSTD7P/otM6cPTfsM9bvEsjv8ki/RFspyc4U+fYt83R7xfTkiwUA4CJoVs0wbx45ndRzK8sDCvgGbr8EfF7Vra0YspWS7EyRgG98ckUCAJAkVkQyTnJTP2Z/YaK2/2lFwsmqF1pY6h92pki87RwAAFKBFZEMU/kHRUk/l5vjUeXsKbrtqi+qcvaUhP0cuTme/pkig5/gZAwAwE4EkQxzbam/rxl1GB5P33Mj0T9TJMntHAAAUoGtmQzz1sef9jejJmJZfc9Vzp4you8dmymSzHYOAACpQBDJMMmecBnt7bqx7RwAANKBrZkMk+wJFztu1wUAINUIIhlmYalfkyeMG/aZSyeM44QLACAjEESyUHIHfAEAMI8gkmEaWzouOln1zNlzamzpSFNFAACMHkEkw9jdrAoAQDoRRDIMzaoAgGxCEMkwsXHsiSZ7eCQFGccOAMgQBJEMwzh2AEA2IYhkIMaxAwCyBZNVMxTj2AEA2YAgksEYxw4AyHS2bc0cPXpU9913n0pLSzV+/HjNnj1bNTU16u3tteuVAAAgw9i2IvL+++8rGo1qx44dmjNnjpqamvTAAw+oq6tLjz76qF2vBQAAGcRjWRe7VD51HnnkEdXV1emjjz5K6vlwOCyfz6dQKKTCwkKbqwMAAKkwkt/fae0RCYVC8vsTz7fo6elRT09P/9fhcDgdZQEAAEPSdnz3yJEjevzxx7Vu3bqEz9TW1srn8/X/KSkpSVd5AADAgBEHkS1btsjj8Qz758CBAwP+mdbWVlVVVWnNmjW6//77E37vzZs3KxQK9f85fvz4yP8dAQCAjDHiHpFTp07p1KlTwz4za9Yseb19w7ZaW1u1ZMkSXXfddXr66aeVk5N89qFHBACAzGNrj0hRUZGKioqSevaTTz7RkiVLtGDBAu3cuXNEIQQAAGQ/25pVW1tbtXjxYs2YMUOPPvqoTp482f93gUDArtcCAIAMYlsQ+fd//3cdPnxYhw8f1vTp0wf8XRpPDAMAAAezba/knnvukWVZcf8AAABI3L4LAAAMIogAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMucR0ASZEopYaWzp0orNbxQVeLSz1KzfHY7osAABcx3VBpL6pTVv3NKst1N3/WdDnVc2qMlWVBw1WBgCA+7hqa6a+qU3rdx0cEEIkqT3UrfW7Dqq+qc1QZQAAuJNrgkgkamnrnmZZcf4u9tnWPc2KROM9AQAA7OCaINLY0jFkJeRClqS2ULcaWzrSVxQAAC7nmiByojNxCBnNcwAAYOxcE0SKC7wpfQ4AAIyda4LIwlK/gj6vEh3S9ajv9MzCUn86ywIAwNVcE0RyczyqWVUmSUPCSOzrmlVlzBMBACCNXBNEJKmqPKi6tRUK+AZuvwR8XtWtrWCOCAAAaea6gWZV5UEtLwswWRUAAAdwXRCR+rZpKmdPMV0GAACu56qtGQAA4CwEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxjp6salmWJCkcDhuuBAAAJCv2ezv2e3w4jg4inZ2dkqSSkhLDlQAAgJHq7OyUz+cb9hmPlUxcMSQajaq1tVUFBQXyeNxxKV04HFZJSYmOHz+uwsJC0+W4Gj8LZ+Dn4Az8HJwhU34OlmWps7NT06ZNU07O8F0gjl4RycnJ0fTp002XYURhYaGj/0PmJvwsnIGfgzPwc3CGTPg5XGwlJIZmVQAAYAxBBAAAGEMQcZj8/HzV1NQoPz/fdCmux8/CGfg5OAM/B2fIxp+Do5tVAQBAdmNFBAAAGEMQAQAAxhBEAACAMQQRAABgDEHEwY4ePar77rtPpaWlGj9+vGbPnq2amhr19vaaLs11tm3bphtuuEETJkzQ5MmTTZfjGk888YRKS0vl9Xq1YMECvfbaa6ZLcp19+/Zp1apVmjZtmjwej5577jnTJblSbW2trr32WhUUFKi4uFirV6/WBx98YLqslCCIONj777+vaDSqHTt26L333tNjjz2m733ve/rWt75lujTX6e3t1Zo1a7R+/XrTpbjGj3/8Yz344IP6y7/8S7399tu68cYbtXLlSh07dsx0aa7S1dWl+fPna/v27aZLcbW9e/dqw4YN2r9/vxoaGnT+/HmtWLFCXV1dpksbM47vZphHHnlEdXV1+uijj0yX4kpPP/20HnzwQZ05c8Z0KVnvuuuuU0VFherq6vo/mzdvnlavXq3a2lqDlbmXx+PRs88+q9WrV5suxfVOnjyp4uJi7d27VzfddJPpcsaEFZEMEwqF5Pf7TZcB2Kq3t1dvvfWWVqxYMeDzFStW6I033jBUFeAcoVBIkrLi9wFBJIMcOXJEjz/+uNatW2e6FMBWp06dUiQS0dSpUwd8PnXqVLW3txuqCnAGy7K0adMmfelLX1J5ebnpcsaMIGLAli1b5PF4hv1z4MCBAf9Ma2urqqqqtGbNGt1///2GKs8uo/k5IL08Hs+Ary3LGvIZ4DbV1dV655139Mwzz5guJSUuMV2AG1VXV+vOO+8c9plZs2b1/+vW1lYtWbJElZWVevLJJ22uzj1G+nNA+hQVFSk3N3fI6seJEyeGrJIAbrJx40a98MIL2rdvn6ZPn266nJQgiBhQVFSkoqKipJ795JNPtGTJEi1YsEA7d+5UTg6LWKkykp8D0isvL08LFixQQ0ODbr/99v7PGxoadNtttxmsDDDDsixt3LhRzz77rH7xi1+otLTUdEkpQxBxsNbWVi1evFgzZszQo48+qpMnT/b/XSAQMFiZ+xw7dkwdHR06duyYIpGIDh06JEmaM2eOJk2aZLa4LLVp0ybddddduuaaa/pXA48dO0aPVJp99tlnOnz4cP/XLS0tOnTokPx+v2bMmGGwMnfZsGGDdu/ereeff14FBQX9q4U+n0/jx483XN0YWXCsnTt3WpLi/kF63X333XF/Dq+++qrp0rLad7/7XWvmzJlWXl6eVVFRYe3du9d0Sa7z6quvxv3P/t133226NFdJ9Ltg586dpksbM+aIAAAAY2g4AAAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGPP/Aag+sZRfw8x4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log(nutpie_target), np.log(estimate_no_decay))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "97dfcd78-fa1b-4b51-8ba4-7377405e8c0a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(13596.147354690482)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lower, upper = linalg.eigvalsh(Sigma, np.diag(estimate_no_decay))[[0, -1]]\n",
    "np.sqrt(upper / lower)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40e0be7e-941c-40a9-8fc0-fc61be2622f3",
   "metadata": {},
   "source": [
    "### Use exponentially weighted windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5c16a91f-8e5a-47a7-8819-23c0d010ffcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "draws_pd = pd.DataFrame(draws)\n",
    "grads_pd = pd.DataFrame(grads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "c4db9922-d208-495e-9c15-f2d2f2f2bd7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "adjust = True\n",
    "ewm_draw_var = draws_pd.ewm(alpha=0.05, adjust=adjust).var().iloc[-1]\n",
    "ewm_grad_var = grads_pd.ewm(alpha=0.05, adjust=adjust).var().iloc[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b5c14240-b090-4b7e-af1e-26ddea2f84b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "estimate_decay = np.sqrt(ewm_draw_var / ewm_grad_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6c9abd28-c3a5-4c18-95e4-6e1d1c7c3c2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fad7e6a7530>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMBVJREFUeJzt3X90VPWd//HXJOYnJIMhhRlKkAh2SxqrBkXxNxTY2C5V2/Vbu8Vqj+0WBE+VPdtW3W7IWpfu0dbuqRW13UL3y+rXfrfrD9Z+ORu3Coq4VCOlafpDMBRqkgoEZjCQBDL3+0ecmB/z485k7nzuzH0+zsk5ncmd3I+Eel9+Pu/P++OzLMsSAACAAQWmBwAAALyLIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAmDNMDyCRSCSizs5OVVRUyOfzmR4OAACwwbIsHT9+XDNmzFBBQeI5D1cHkc7OTtXU1JgeBgAASMPBgwc1c+bMhNe4OohUVFRIGvoHqaysNDwaAABgRzgcVk1NzfBzPBFXB5HockxlZSVBBACAHGOnrIJiVQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxrm5oBgAAnDEYsbSro0fvHO/TtIpSLaitUmFB9s91I4gAAOAxW9u61LylXV2hvuH3gv5SNS2vU2N9MKtjYWkGAAAP2drWpVWbW0eFEEnqDvVp1eZWbW3ryup4CCIAAHjEYMRS85Z2WTG+F32veUu7BiOxrnAGQQQAAI/Y1dEzbiZkJEtSV6hPuzp6sjYmgggAAB7xzvH4ISSd6zKBIAIAgEdMqyjN6HWZQBABAMAjFtRWKegvVbxNuj4N7Z5ZUFuVtTERRAAA8IjCAp+altdJ0rgwEn3dtLwuq/1ECCIAAHhIY31QG1Y0KOAfvfwS8Jdqw4qGrPcRoaEZAAAe01gf1NK6AJ1VAQCAGYUFPi2cM9X0MFiaAQAA5hBEAACAMY4GkfXr1+uiiy5SRUWFpk2bpuuuu06/+93vnLwlAADIIY4GkW3btmn16tV69dVX1dLSotOnT2vZsmXq7e118rYAACBH+CzLytrJNocOHdK0adO0bds2XXnllUmvD4fD8vv9CoVCqqyszMIIAQDARKXy/M5qjUgoFJIkVVVlr2MbAABwr6xt37UsS2vXrtXll1+u+vr6mNf09/erv79/+HU4HM7W8AAAgAFZmxFZs2aN9uzZoyeeeCLuNevXr5ff7x/+qqmpydbwAACAAVmpEbn99tv19NNPa/v27aqtrY17XawZkZqaGmpEAADIIanUiDi6NGNZlm6//XY99dRTevHFFxOGEEkqKSlRSUmJk0MCAAAu4mgQWb16tR5//HE988wzqqioUHd3tyTJ7/errKzMyVsDAIAc4OjSjM8X+/CcjRs36pZbbkn6ebbvAgCQe1y1NAMAABAPZ80AAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAY84wPQAAAFIxGLG0q6NH7xzv07SKUi2orVJhgc/0sJAmgggAIGdsbetS85Z2dYX6ht8L+kvVtLxOjfVBgyNDuliaAQDkhK1tXVq1uXVUCJGk7lCfVm1u1da2LkMjw0QQRAAArjcYsdS8pV1WjO9F32ve0q7BSKwr4GYEEQCA6+3q6Bk3EzKSJakr1KddHT3ZGxQygiACAHC9d47HDyHpXAf3IIgAAFxvWkVpRq+DexBEAACut6C2SkF/qeJt0vVpaPfMgtqqbA4LGUAQAQC4XmGBT03L6yRpXBiJvm5aXkc/kRxEEAEA5ITG+qA2rGhQwD96+SXgL9WGFQ30EclRNDQDAOSMxvqgltYF6KyaRwgiAICUmG6xXljg08I5U7N2PziLIAIAsI0W68g0akQAALY42WJ9MGJp574jemb329q57wgdUj2EGREAgKTESy7JWqz7NNRifWldIOVlGmZZvI0gAgBIGgZSabGeSv1GdJZlbMCJzrKwGyb/sTQDAB5nZ8nFiRbrHGQHiSACAJ5mNwxUTy6x9fNSabHOQXaQCCIA4Gl2w4AsJWyxLklnlhel1GKdg+wgEUQAwNPsPuQP9/araXldzJmTqKMnTqmlvdv2vTnIDhJBBAA8LZUwsLQuoCnlRXGvie6csVvTke5Bdmz1zS/smgEAD4uGge5QX8zZDp+GznJZUFulV986omMnTsX9WanunIkeZLdqc6t80qj7xzvIjq2++YcZEQDwMLun2ra0d2v1v7Xa+pmp1HSkcpCdkw3VYA4zIgDgcdEwMHamIfDeTIOkmL0+4oku99g9k6axPqjFH56u/71zv/7Qc0JnVZXrpoWzVXzG+/+t7GRDNZhFEAEAxD3VVpIu/6ef2w4hknS0tz+lJZRY137/xb26/vwPakldQAtqqxxrqAbzCCIA4GHJZi127juSMADE8rc/3aPe/sFx78fqlhqvs2pP7yn9y479+pcd+xX0l+qa+oCte7PVN/cQRADAo+zMWqTzYI8VQqTxSyiStO7ZXyedbekO9elHO/bbujdbfXOPo8Wq27dv1/LlyzVjxgz5fD49/fTTTt4OAGCT3cLPTD/YRy6hPPTzveoO99v6jE9SotKPeFt94X6OBpHe3l6dd955euihh5y8DQAgBckKPy293w8kWa+PdLW0d+vB539v+3pLUrRdSKLdPRSq5h5Hg8g111yjb37zm/rUpz7l5G0AAClIVvgpvT9rkWh770T85LU/pvW5Wy+bbWurL3IHNSIA4DHdYXt1H9HrltYFdMeSc7Rxx34dOxm/oVkq3u0/ndbnltQFdPcn6mxtC0ZucFUQ6e/vV3//++uF4XDY4GgAID/1vJu8LiN6XayCVn9ZkeoCldrz9jH1DsQuTM20kR1eCwt8bNHNI67qrLp+/Xr5/f7hr5qaGtNDAoC8UzWp2NZ1fzx6ImZBa+jkKe3sOJLVECJRA5KvXBVE7rrrLoVCoeGvgwcPmh4SALjSRA5+C/jLbF33761vp9TIzCnUgOQ3Vy3NlJSUqKSkxPQwAMDVtrZ1ad2z7aNqPQKVpVr3SXsHv0V3wiQqWPVJOt6XXh1HMj5J/vIihd47QC/WYXd3LPmQZleXUwPiAY4GkXfffVd79+4dft3R0aHdu3erqqpKs2bNcvLWAJCXtrZ1aeXm8YfPdYf7tHJzqx6xMXNQWODTJ88L6tHtHXGvcXImxJL0rU+dK0lxz7dh9sM7fJZlOfb37cUXX9SiRYvGvX/zzTdr06ZNST8fDofl9/sVCoVUWVnpwAgBIHcMRizN/2aLjp2Iv3PlzPIivfZ3SxPOIPxsT6due/wNJ4Zoy8gx2j0YD7kllee3ozMiV199tRzMOQDgKa++dSRhCJGkoydO6dW3juiyudXD74182Hcc6tU///ebTg81oaMnTg0fTscOGLiqRgQAEN/OfUdsXffKvsPDQSTW9ls34HA6RBFEACBn2Jth3rRjv879oF+SYp5s6wYcTocogggA5IiFZ1froRf2Jb2ud2BQKze3akp5ketCyMjGZIDksj4iAID4LpkzVVPKi2xfn6yeJNtoTIZYCCIAkCMKC3zD215zEY3JEAtBBABySGN9UI+saJC/NLdW1r/xiXl6+WuLCSEYhyACwPMm0i492wYjlvxlxfrcJbnVFLK6ooTlGMSUW5EaADIs1vbWoAu6e8Zq9NXS3m18K+6Xr6zVs7/sSnkM7JJBPAQRAJ61ta0r5vbW7lCfVm1uNVbPECscTSkvckXx6QWzztRXG+cNh6TqySX6m5/s1p/C/TF36LBLBsmwNAPAkwYjlpq3tMd8eEbfa97SnvVlmmg4Gjvj4IYQ4tPQn4kkLZwzVdee/0FdNrda6z75keHvj71eYpcMEiOIAPCkXR09CZcXLEldoT7t6ujJ2pgShSM3iPdn0lgf1IYVDQr4Ry+/sEsGdrA0A8CT7LYYH3tdskPa0jnELfqZHXsPua4Veyyx/uwa64NaWhfgADukjCACwJPsFk+OvC5ZYWs6ha9uPQsmkXh/dhxgh3SwNAPAkxbUVinoLx1X1xDl01CIiBZZxqvd6Ar1aeXmVt33XHvC7/9sT9e4e8T7mW419s8EyASCCABPKizwqWl5naTkRZZ2ajd+8FJHwu+veaJVP9vTOfzarfUgiz/8AUkUniJ7CCIAPMtukWWywlY7IpZ02+NvaGtbV8Z+ZqZ9+cpa/eiWBXqEwlNkETUiADzNTpGl3cJWO5q3tGtpXSCjPzNTnv1ll77aOI/CU2QVQQSA5yUrssxkV9CuUJ827ejQhwOVGfuZmRLdmrtwzlQKT5E1BBEASCJa2JqppZR7n/uNpleUaEp5kUInTrmqTsSNMzXIb9SIAEASIwtbM+VPx/t1zGUhROJMGGQfQQSAJ6V64m5jfVDf++wFGR/HpOJC+UvdMTnN1lyY4I6//QCQRek2HvvG021Jf/aSeR/Q8785ZHssvQODtq91kk9szYUZzIgA8JR4TcSiJ+5Gt9eO/czKza06djL5wXO/7jyuh//qAgX9ubPEEWRrLgxiRgSAZyQ7cTd6uuzSusDwzMBgxNLXfrrH9j26Qn3qCvXpq40f1k9fP6iX9x7JyNgzaUpZkVYvmqPqySUK+MvYmgujCCIAPCOVE3ejW1dffeuIQidPp3Sfe5/7zUSG6Zho1PjWp89l9gOuwdIMAM9I58TdnfvcN6Nh19hJDrqjwo2YEQHgGemcuCvXbbC1584l52jV1XP1+h+O0h0VrkYQAeAZ0cZk3aG+uPFiSlmRIpalwYg11F307Go99MK+rI5zIgKVJVr3yY8Mz3rQHRVux9IMAM9IdOJu1LGTp/S5H/6P5t/bon9+/k1dVFulKeVF2RvkBH37f53P0gtyCkEEgKfEO3F3rGMnT+nB53+vBf/4vD5z4cwsjW7iDr/bb3oIQEoIIgA8YWQnVX9Zsbb97SL9260Xa0pZ4tmOYydO6dHtHZpUXJilkU4MLdqRa6gRAZD3YnVSrZpUpM9fcpatJmWSezqgJkKLduQiggiAvBbtpDq2OLWn95S++997jYzJKbRoRy5iaQZA3krUSTXf/GXDTIpUkZOYEQFgzGDE0q6OHkf6XAxGLG3a0ZGwk2o+uWwu23SRmwgiAIxI5wRcOwYjlh76+ZvauGO/7fqPfBDwl5keApAWggiArItXtxE9ATfdNuRb27r09f/4lY6d8E4AkShSRW6jRgRAViU7AVcaOgF3MJJaZUc03ORLCFky7wO2r6VIFbmMIAIgq1I5AdeufCxKDdpcarlzyTkUqSKnEUQAZJXdE3C7Qydt/8xk4SYXnT9zioL+0rit6KWhJZk1i8/J2pgAJxBEAGSV3c6f9z73G21t67J1rd1wk0tmnFke91wc33tfLMkgHxBEAGTVApuHyB3tHdCqza22wkj15JJMDM01ApUlWlBbFfdcnIC/NO2CXsBt2DUDwJUsDf1Xf/OWdi2tC4z7L/+RPUjeCefXjEjf6Yha2rvVWB9UY31QS+sCjvVbAUwjiADIql0dPbZ3tkQLV19964gKfL7hB/HR3n7d+9xv8q4uJCp04tSobcyFBT4tnEPDMuQnggiArEqnnuMLG3dpYDCf9sQklmw2CMgn1IgAyKp0jqn3UgiJSmcbM5CLCCIAsmpBbVXSbal4Xz7uCAJGIogAyKrCAl/cbaleMrnE3r9+05lBAnIJQQRA1kW3pU6rKDY9FCOmTirWL+5ZlnBmyCfOkIE3EEQAGOO12o9oI7L7rq9XWXFhwoZlEg3L4A0EEQATMhixtHPfET2z+23t3HfE1mF10QPqjubJAXV2jW1ERsMygO27ACZga1uXmre0j+rnEfSXqml5XdyHaD4eUJfIVz42V2d/YHLcRmQ0LIPXEUQApCU6qzE2UHSH+kY14xorHw+oi8cn6Sev/VEvf21xwmBBwzJ4GUszAFKWaFYj+l7zlvaYyzRe2o5KLxAgOYIIgJQMRixt2tGRcFYj1gM4Wkvy5p+OZ2GU7uKl8AWkiqUZALbFqglJJPoATvVz+YZeIEB8BBEAtsSrCUlkWkVpWp/LFz4N7YChFwgQH0EEQFLp7HQpLy7U/LPO1FX3v+DZECLRCwRIhiACIKl0drqcGBjU2p/s9uxyTCDJNmYAQwgiAJJKt9jyuT1dGR5JbvjGJ+bplstqmQkBbMjKrpmHH35YtbW1Ki0t1fz58/XSSy9l47YAMiTdYkuvLclEz4chhAD2OR5EnnzySd1xxx2655579MYbb+iKK67QNddcowMHDjh9awAZsqC2KuEBbYkUF3rjgUxNCJAex4PId77zHd1666364he/qHnz5um73/2uampqtGHDBqdvDSBDCgt8cQ9oSyZfD7abUl406jXnwwDpcbRGZGBgQK+//rq+/vWvj3p/2bJleuWVV8Zd39/fr/7+/uHX4XDYyeEBSEH0gDYv9wORpCllRfrWp8/lfBggQxwNIocPH9bg4KCmT58+6v3p06eru7t73PXr169Xc3Ozk0MCMMJgxErpYTr2gLaW9j/pPz1WkPr9zzXosrnVksT5MEAGZGXXjM83+l9slmWNe0+S7rrrLq1du3b4dTgcVk1NjePjA7wonZNzpdEHtJWcUeCZIBJtTnbJ2YQPIJMcrRGprq5WYWHhuNmPd955Z9wsiSSVlJSosrJy1BeAzIt2Ox27xBI9OXdrW+JwMRix9NLvD+lv/u8vnRyma1CICjjH0SBSXFys+fPnq6WlZdT7LS0tuvTSS528NYA4kp2ca0n6+k9/pR17D8c8PXdrW5fmf7NFN/1ol3r7B50erhElY3b6UIgKOMfxpZm1a9fqpptu0oUXXqiFCxfqscce04EDB7Ry5Uqnbw0gBjtdUo+dPKXP/fB/xi3VbG3r0srNrdkYplH9g5buXHKOZldPohAVcJjjQeQzn/mMjhw5on/4h39QV1eX6uvr9bOf/UxnnXWW07cGEEMqXVKjSzUbVjRoaV1A655td3Bk7vJ/fnFQL39tMQEEcFhWilVvu+023Xbbbdm4FYAkUumSammoPqJ5S7sqSovUHfbOtt2uUJ92dfSwMwZwWFZavANwj1S7pFoaeijv3HfEyWG5Urpn7ACwjyACeEy6XVL/2HPCmQG5WLpn7ACwjyACeFC0S2rAb/9B+/QvOx0ckbtED69bUFtleihA3stKjQgA91laF1BFaZFe2XtYm17Zr96B/NyKmy56hgDZQRABPCTa0r2lvVtP7+5UT++A6SEZ49P7xbgju6UU+KQvXVFLzxAgSwgigEfEaunuZVPKi3T0xKlxjd0sS3pse4cumHUmYQTIAmpEAA+I19Ldq+75+DyVnBH7X3/RYNK8pT1mZ1kAmUUQAXLAYMTSzn1H9Mzut7Vz35GUHpCJWrp7TbQItS5Yqe5wf9zroluWd3X0ZG1sgFexNAO4ULSW453jfdp/uFdP7Dow6sFp55TcKDst3XPR5JJCvZvGWTdNy+t0uDd+CBmJPiKA8wgigMvYqeUY2Xo9WRjJx4fplPIi7bp7iV7/w9H3wtoJfff53yec9RkZ3uw2Z6OPCOA8ggjgItFajmTLKCNbry+tCwxvMx05kxI9rC0fH6ZfuLRWxWcUjGq//meByeMC3NRJxbr2/BlaWhcYdXBdtLtsd6gv5p+1T0Mn7tJHBHAeQQRwiVRrOUbWMSycMzXmTErQX6pvfGJewodurplSXqQ1i+eOe7+xPqildYFxQSxWL5Bod9lVm1vHbd+NXk0fESA7CCKAS6Rby/H/2rq0q6Mn5tJEV6hPtz3+hr50Ra1++FJHZgZq2Lc+dW7cgFBY4LN1SN1gxJK/rFhfuGz2uH4qgRTqbwBMHEEEcIl0azn+decfkl7zLy936ItX1Oonrx1U6OTptO5j2pTyIn3rU+dOOCDEmjmqmlSk68//oJaMWcIB4DyCCOASTtZyRCzpBy91aFJJoWP3cMqU8iJ94dJarVk8d8IBIV4NztHeU/rRjv26iBACZB19RACXiBZQOvkY7E1ju6tJ3/jEPO26e4kW1FbpP/d0ptxDZaRENTg0MQPMYUYEcIlEBZReFPSXKugv1VX3vzCuADedGo5kNThji38BZAczIoCLNNYHtWFFgwL+0cs0U8qK9JcNM3XlOdWGRpZ9f/HRoFY//sa48BDtobK1rSuln2e3Bicf+64AbkYQAQyK1bq9sT6ol7+2WHcu+ZCmlBVJko6dPKV/b/2jXtp72PCIs+enrW9ndBnFbg1OPvZdAdyMpRnAkHh9P5qW10lSzO24lofWa0ZuqR0rnWUUmpgB7kQQAQyIt3ujO9SnlZtbNam40PM1InaksoxCEzPAnViaAbLMzu6N3oHc2t2SST4N9fWwI9VllHg1OAF/qa1zewBkHjMiQJbl62m4mfTNa+t173O/cWQZJZVW8ACcRxABsoxdGfGN3JpbUOBzbBnFbit4AM4jiAAZEuvk21gPSnZljBbvhNzoMsrYgl7OggHyC0EEyIBEO2DGPjCT7d7IZ9HTgM+cVGJrWYRlFCD/+SzLvRsCw+Gw/H6/QqGQKisrTQ8HiCneDpjoo3LDioZxD9Ojvf1a/fgbkrzTQfUbn5inWy6rJUQAHpDK85sgAkzAYMTS5f/087jFpz4NHdpWckaBusP9w+8H/aX65HlBPfnaH3XsxKlxn3Ht/ynTEC0sfflriwkhgEek8vxmaQaYADvnlxwdEzSkoX4hj27viPuZfEF/DgDJEESACUh3B0w+hY1EKCwFkAxBBJgAdsDE9vmFZ+ma+iCFpQCSorMqMAHRHTA8ake7pj6ohXOmEkIAJEUQASYgen6JJMLIewp80tEEB9YBwEgEEWCC4p1fEvSXakp5kecCSsSSVj/eqq1tXaaHAiAHsH0XyJBYnVVb2rtj9hjJd2zZBbwtlec3MyJAhkTPL7n2/A9qQW2VdnX0qP90RHcs+ZD8ZfZOk80XlqSuUJ92dfSYHgoAl2PXDJBhsdq9TyouNDgiczjgD0AyBBEgg+K1e+8dGDQyHtPY3gwgGZZmgAwZjFhq3tLumXqQRKUfPg0V6y6orcraeADkJoIIkCHJ2r3nmy9dUSufxm9bpq07gFQQRIAM8Uo9xJnlRXpkRYPu+nhdzG3LAX+pNqxooK07AFuoEQEyJF/rIT53cY2mlBXL5xvaFXTJ2e93TG2sD2ppXWDctmVmQgDYRRABMmRBbZUClSXqDvebHkpGLagd2pIcT3TbMgCkg6UZIEMKC3z67IJZpoeRcfk60wPAHZgRAUaI1R010TLDYMTSq/uOaOdbhyX5ZOXRnplod1R2vgBwEkEEeE+sRmRBf6maltepsT44LnScUeDTplc6dOzkaXODToNPShqX2PkCIFsIIoDiNyLrCvVp5eZWfezDH9D/dBzVu/25FTrGWv7RgF77w7FRYevM8iJZko6dODX8XmBEAAMAJxFE4Hl2GpH9928PZW08TvrF/qPa/tXFev0PR0ctP0li5wsAIwgi8DwvNSLrDvfr9T8cjbnLhZ0vAExg1ww8zyuNyKK89s8LwN0IIvC8/YdPmB5CVrEdF4CbsDQDTxs4HdGPd+43PYysYDsuADdiRgSetbWtS5esf149vQOmh+I4tuMCcCtmROBJ8bbr5iu24wJwK4IIPMfOdt1cFp3vuGPJhzS7upztuABcjSACz8n37brMfgDIJQQReE4+bl+tmlSkb/zFRxSoZPYDQG4hiMBz8mn7ajRu/OP15zIDAiAnEUSQd5KdoHu0t18FPimSg0UiVZOK1NPLmTAA8gdBBHkl1gm6VZOKdd35M7S0LqCjvf1a/fgbOVmoOnVSsXbe9bFx58SwDAMgl/ksy3Ls38n33XefnnvuOe3evVvFxcU6duxYSp8Ph8Py+/0KhUKqrKx0ZpDIG3a25ObqTIgkPfxXDfr4R5n5AOB+qTy/HW1oNjAwoBtuuEGrVq1y8jbwkMGIpZ37juiZ3W9r574jGnwvVdjdkpurIeTLV9YSQgDkJUeXZpqbmyVJmzZtcvI28IhYyy6BylJ9dsEsnRoczMstuZNKCnX/pz+qj390humhAIAjXFUj0t/fr/7+/uHX4XDY4GjgJvGWXbrDfXrw+d8bGZOTyosL9eUr52jN4rnUgADIa646a2b9+vXy+/3DXzU1NaaHBBfI906osfzg8xfqK0vOIYQAyHspB5F169bJ5/Ml/HrttdfSGsxdd92lUCg0/HXw4MG0fg7yS753Qh3JJynoL9UlZ081PRQAyIqUl2bWrFmjG2+8MeE1s2fPTmswJSUlKikpSeuzyF/52Ak1EU7IBeAlKQeR6upqVVdXOzEWIKZ86oSaSJDmZAA8yNFi1QMHDqinp0cHDhzQ4OCgdu/eLUmaO3euJk+e7OStkUcW1FYp6C9Vd6gvb+pEpk4q1t3XfFjHTp5S1eQSzogB4FmOBpG///u/149//OPh1xdccIEk6YUXXtDVV1/t5K2RRwoLfGpaXqdVm1vlk3I6jERjxn3X1zPzAQByuLPqRNFZ1X2SnePipFh9RNxubCdXll8AeEEqz29X9RGBu8UKAtl8sDbWB7W0LqBdHT3asfewHnphr+P3TFc0mj302Qt05qQSzoYBgDiYEYEt8RqKRR+pG1Y0OBJG4s3ADEYsXf5PP3ft7AgzHwC8jBkRZFSihmKWhsJI85Z2La0LZPS/9pPNwHzyvKAe3d6RsftNlM8nPfCX52nGlDJmPgDAJld1VoU7JWsoZknqCvVpV0dPxu4ZnYEZe9/uUJ9WbW7V+p+16zEXhRBJsizp2IkBLZwzlRACADYRRJCU3YZimWo8lmwGRpJ+8FKHK3fP7D9ywvQQACCnsDSDpOw2FMtU4zE7MzDurWxy7cAAwJUIIkgqWUMxn6SAf6iQNF0ji1Lf/NO7af8c01iRAYDUEESQVKKGYtHn7kTOR8nF/iDx/HjnAS2cU81uGQCwiRoR2NJYH9SGFQ0K+EcvvwT8pRPauhuvKDWXNW9p12CEJRoAsIMZEdg2sqFYJhp0DUYsrXs2dlFqLusK9enVfUd02TkcDgkAyRBEkJLCAp8WzpmakZ/10M/fVHc4f2ZCRlr9eKu+9elzWaIBgCRYmoERW9u69ODzb5oehmOOnTylVZtbtbWty/RQAMDVCCLIumifEC+gXgQAEiOIIOuS9QnJF050nAWAfEMQQda1tHebHkJWZarjLADkI4IIsmprW5d+tGO/6WGMMrnkDFWWFCa8xqf0m5VlquMsAOQjggiyxq21IZ+7uEZvNP25vvKxuTG/H80fX7qiVr4Rr5Pxaei04Il0nAWAfEcQQda4tTbk2V8O7Wy5c+mf6ZEVDQrGadp218frYjZ1iyUTHWcBwAvoI4KscWutRLSgdOGcqUmbtsX6/tHeAd373OgW9QF/qZqW19FHBACSIIggbSMPqrPTZXX/4d4sji41I0NSsqZtsb7/5/WZ6zgLAF5CEEFaYh1UF0wwCzAYsfSjHR3ZHOK4A/oSmWhBaSY7zgKAl1AjgpTFO6iuO9Q3rpvoYMTSzn1H9J2W3yl08rRjY/pIsEJVk4pHvTelvCjp5ygoBQCzmBFBSoYOqvt1zJmG6Ht3P/UrnTwV0f7DvfrXnft19MQpx8f1d3/xES2orRpeHqmeXKK/+cluW5+loBQAzCGIICUP/XyvusP9Ca/p6T2lO5/cnZ0BSZo6qXi4JiO6PLJz35Gk45SkO5Z8iIJSADCIpRnYNnRQ3e9ND2Oca8+fMW5Gw+4OndnV5U4MCQBgE0EEtgycjujup9pMDyOmpXWBce/ZLT6l6ykAmEUQQVJb27p0yfrn1dM7YHoo48QrNF1QW6WgvzRuF1SKVAHAHQgiGN7Z8szut7Vz35FRx9ZHd8j09DpfcDpW0F+qL18Zu6169L14haaFBT41La8bvnbsZ5XgswCA7KFY1eMS9QNZWhdQ85Z22704MuWa+un6/MLa4QLUC2adOW6MdjqXNtYHtWFFQ1qfBQBkh8+yrGw/Z2wLh8Py+/0KhUKqrKw0PZy8E53tGPsXIDpHcMeSc/Tg829mdUyTSwr1y6Y/HzdTkWoX10x9FgCQulSe38yIeFT0JNx4/UB8kn70cnY7oUpSxLLU0t49brZiIp1L6XoKAO5FjYhHJTsJ15IU6nOuE2o8JwYiWjmmOysAIH8RRDzKxEm4n274oO1rm7e0jyqaBQDkJ4KIR2W7f4bPJ106p9r29V2hPu3q6HFwRAAANyCIeFSyPhuZZlnSsROp9SExMWsDAMgugohHJeqz4ZSqScUK+u3PxND1FADyH0HEw6J9NgIphIOJCPjL1LS8LmnwoespAHgH23c9rrE+qKV1geE+G1XlxVr5b6+rt38w4/c62jugj390KPx8/T9+pWMnxndrpespAHgLMyJ5LFHr9pGifTZKzijQV3+6x5EQIkn3Pje0E6axPqjX/26p7lxyjqaUFY26JuAv1YYVDXQ9BQCPYEYkTyVq3R7rIR+vy6odiz78AV0+p1r3PvebhNdFd8IsnDNVhQU+fWXJh7Rm8Tl0PQUADyOI5KF4oaI71KdVm1vHzTgk6rJqx19fMcf2Dpex19H1FAC8jaWZHBZr6SVZ63ZpfLOwZF1Wk4nOZtjBThgAwEjMiOSoeEsvN15Uk7R1+8glEmni/TqiSypBf2nce/s0VP/BThgAwEjMiOSg6NLL2Id+d6jP9mm5I8NHurMUI7fZFhb49MnzEheYshMGADAWQSTH2Fl6sWNk+JhIl9VouNja1qXHtsc/rfevr6xlJwwAYByCSI6ZaD1HrGZhI7us2lU1qWi46DVZsatP0rO/7OIQOwDAOASRHJNKPcfYGY5EzcKiXVanVxQn/blTJxXr1buWDM9wJAtHI+tSAAAYiSCSY+zWc9y55EPjWrfbaRbm88X/K+F77+u+6+tVfMb716W7dRcAAHbN5JhoPUd3qC/mUkh0d8qaxXO1ZvFc283C7DQ0C/hLdeNFs9R/OqKd+44M/zy27gIA0kUQyTHReo5Vm1vjXjNy6cVOszA7Dc0mlxQqEonowed/P/xetFPr0rqArXDE1l0AwFgszeSgxvqg/vrKWo2d3Cjwpbc7xU4B7Lv9g/rT8YFR70U7tba0dw8Xu6ZSlwIAAEEkB0W3yo7dhGJZ0mPbO7S1rSuln5du7cbITq1L6wLasKIhrboUAIB3sTSTY5L1EfHp/WBgdwZiIrUbI3fENNYHtbQuwCF2AADbCCI5JpWtsnYPk0tWAGtHdFaFQ+wAAKlgaSbHOLFVdmRDs3TnLtgRAwBIB0Ekxzi1VTba0GxsjUfQX6op5UVxA0qsTq0AANjF0kyOmX/WmfIp8bkyvveuS1W8Go+W9m6t2tw67r7siAEATBQzIjnmFx09Ses4rPeuS0e0xuPa8z+ohXOmqrDAF3e2hB0xAICJYkYkx+x867Dt6y47pzpj92VHDADACQSRnGP3wZ/5gMCOGABApjm2NLN//37deuutqq2tVVlZmebMmaOmpiYNDAwk/zDishsECAwAgFzg2IzIb3/7W0UiET366KOaO3eu2tra9KUvfUm9vb164IEHnLpt3rvk7KmaUl6kYydOxb3mzPIiXXI2QQQA4H4+y7LS7WGVsvvvv18bNmzQW2+9Zev6cDgsv9+vUCikyspKh0eXO7a2dWllgkPvHqGAFABgUCrP76zumgmFQqqqit9vor+/X+FweNQXxmusD+qRFQ0KVI7v+UEIAQDkkqwVq+7bt0/f+9739O1vfzvuNevXr1dzc3O2hpTT2MUCAMgHKS/NrFu3LmlY+MUvfqELL7xw+HVnZ6euuuoqXXXVVfrhD38Y93P9/f3q7+8ffh0Oh1VTU8PSDAAAOSSVpZmUg8jhw4d1+HDiXhazZ89WaenQskFnZ6cWLVqkiy++WJs2bVJBgf3VIGpEAADIPak8v1NemqmurlZ1tb1GWW+//bYWLVqk+fPna+PGjSmFEAAAkP8cqxHp7OzU1VdfrVmzZumBBx7QoUOHhr8XCAScui0AAMghjgWR//qv/9LevXu1d+9ezZw5c9T3srhjGAAAuJhjayW33HKLLMuK+QUAACBx+i4AADCIIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGPOMD0AEwYjlnZ19Oid432aVlGqBbVVKizwmR4WAACe47kgsrWtS81b2tUV6ht+L+gvVdPyOjXWBw2ODAAA7/HU0szWti6t2tw6KoRIUneoT6s2t2prW5ehkQEA4E2eCSKDEUvNW9plxfhe9L3mLe0ajMS6AgAAOMEzQWRXR8+4mZCRLEldoT7t6ujJ3qAAAPA4zwSRd47HDyHpXAcAACbOM0FkWkVpRq8DAAAT55kgsqC2SkF/qeJt0vVpaPfMgtqqbA4LAABP80wQKSzwqWl5nSSNCyPR103L6+gnAgBAFnkmiEhSY31QG1Y0KOAfvfwS8Jdqw4oG+ogAAJBlnmto1lgf1NK6AJ1VAQBwAc8FEWlomWbhnKmmhwEAgOd5amkGAAC4C0EEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYIyrO6taliVJCofDhkcCAADsij63o8/xRFwdRI4fPy5JqqmpMTwSAACQquPHj8vv9ye8xmfZiSuGRCIRdXZ2qqKiQj6fNw6lC4fDqqmp0cGDB1VZWWl6OJ7G78Id+D24A78Hd8iV34NlWTp+/LhmzJihgoLEVSCunhEpKCjQzJkzTQ/DiMrKSlf/JfMSfhfuwO/BHfg9uEMu/B6SzYREUawKAACMIYgAAABjCCIuU1JSoqamJpWUlJgeiufxu3AHfg/uwO/BHfLx9+DqYlUAAJDfmBEBAADGEEQAAIAxBBEAAGAMQQQAABhDEHGx/fv369Zbb1Vtba3Kyso0Z84cNTU1aWBgwPTQPOe+++7TpZdeqvLyck2ZMsX0cDzj4YcfVm1trUpLSzV//ny99NJLpofkOdu3b9fy5cs1Y8YM+Xw+Pf3006aH5Enr16/XRRddpIqKCk2bNk3XXXedfve735keVkYQRFzst7/9rSKRiB599FH9+te/1oMPPqhHHnlEd999t+mhec7AwIBuuOEGrVq1yvRQPOPJJ5/UHXfcoXvuuUdvvPGGrrjiCl1zzTU6cOCA6aF5Sm9vr8477zw99NBDpofiadu2bdPq1av16quvqqWlRadPn9ayZcvU29tremgTxvbdHHP//fdrw4YNeuutt0wPxZM2bdqkO+64Q8eOHTM9lLx38cUXq6GhQRs2bBh+b968ebruuuu0fv16gyPzLp/Pp6eeekrXXXed6aF43qFDhzRt2jRt27ZNV155penhTAgzIjkmFAqpqqrK9DAARw0MDOj111/XsmXLRr2/bNkyvfLKK4ZGBbhHKBSSpLx4HhBEcsi+ffv0ve99TytXrjQ9FMBRhw8f1uDgoKZPnz7q/enTp6u7u9vQqAB3sCxLa9eu1eWXX676+nrTw5kwgogB69atk8/nS/j12muvjfpMZ2enGhsbdcMNN+iLX/yioZHnl3R+D8gun8836rVlWePeA7xmzZo12rNnj5544gnTQ8mIM0wPwIvWrFmjG2+8MeE1s2fPHv7fnZ2dWrRokRYuXKjHHnvM4dF5R6q/B2RPdXW1CgsLx81+vPPOO+NmSQAvuf322/Xss89q+/btmjlzpunhZARBxIDq6mpVV1fbuvbtt9/WokWLNH/+fG3cuFEFBUxiZUoqvwdkV3FxsebPn6+WlhZdf/31w++3tLTo2muvNTgywAzLsnT77bfrqaee0osvvqja2lrTQ8oYgoiLdXZ26uqrr9asWbP0wAMP6NChQ8PfCwQCBkfmPQcOHFBPT48OHDigwcFB7d69W5I0d+5cTZ482ezg8tTatWt100036cILLxyeDTxw4AA1Uln27rvvau/evcOvOzo6tHv3blVVVWnWrFkGR+Ytq1ev1uOPP65nnnlGFRUVw7OFfr9fZWVlhkc3QRZca+PGjZakmF/Irptvvjnm7+GFF14wPbS89v3vf98666yzrOLiYquhocHatm2b6SF5zgsvvBDz7/7NN99semieEu9ZsHHjRtNDmzD6iAAAAGMoOAAAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABjz/wHrvm3JSgb1QwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log(nutpie_target), np.log(estimate_decay))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0307b19d-0ee3-4062-a5d4-c19f9def0380",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(13704.157180258602)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lower, upper = linalg.eigvalsh(Sigma, np.diag(estimate_decay))[[0, -1]]\n",
    "np.sqrt(upper / lower)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c69d3f0-c48c-4962-befd-29034b3a3580",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e0a19d8-8050-4e68-a22b-c3c728802376",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (Pixi)",
   "language": "python",
   "name": "pixi-kernel-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.12.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}