Simple examples of simulation based inference

simple-sbi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a7a70d0f-b789-4d98-a5b5-ef8030860ec0",
   "metadata": {},
   "source": [
    "# simple simulation based inference\n",
    "\n",
    "this notebooks explains simulation based inference with simpler methods than the more well known packages like mackelab's. \n",
    "## what is sbi actually doing?\n",
    "The approach is very general:  the simulator generates many samples of pairs of (parameter values, data features), and a Bayesian regression of parameters on data features is done.  In a very simple case this could be just least squares (for point esimates) or a simple linear Bayesian model. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b96bac7-6ebc-4dbf-b1af-fce605b4cefd",
   "metadata": {},
   "source": [
    "## setup\n",
    "\n",
    "in case of least squares we don't need anything other than numpy,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fb5d1b3-af2a-4b6f-9c84-be01af89f3ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d2b13a5-145a-4494-9b85-a82686885117",
   "metadata": {},
   "source": [
    "some dimensions of the problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b21c3749-4fff-48d4-bab7-6f8c466aaffc",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_params = 10\n",
    "num_samples = 1000\n",
    "num_features = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4b11ea9-d7d2-4fd5-9e1e-c951cee20aea",
   "metadata": {},
   "source": [
    "a simple simulator of a latent linear system where we are going to estimate its initial conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "20f3bb71-738c-45ed-9964-894cc17f0ce0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -85435.66261261,  308974.00453229,  243658.29502642,\n",
       "       -215082.5979462 , -388077.38372195])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.randn(num_params, num_params)\n",
    "G = np.random.randn(num_features, num_params)\n",
    "\n",
    "def simulator(p):  # linear system\n",
    "    x = p\n",
    "    for i in range(10):\n",
    "        x = A @ x\n",
    "    return G @ x\n",
    "\n",
    "simulator(np.random.randn(num_params))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79b8e9ad-0410-4911-bae0-fffab1ec7d38",
   "metadata": {},
   "source": [
    "sample parameters and run simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "6403a706-e988-44b8-bcdf-c23a6a69e6ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = np.random.randn(num_samples, num_params)\n",
    "features = np.array([simulator(p) for p in params])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f93c4df2-da45-41b1-a9ac-a6e37957b7cd",
   "metadata": {},
   "source": [
    "do least squares and compute losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "e155f089-2f5a-460b-a096-5aa2d5781ec2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5077.507087203544"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predictor, *_ = np.linalg.lstsq(features, params, rcond=None)\n",
    "assert predictor.shape == (features.shape[1], num_params)\n",
    "params_est = features @ predictor\n",
    "loss = np.sum(np.square(  params - params_est ))\n",
    "loss"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "812ea6cd-3caa-4960-9468-160ad81fdd0d",
   "metadata": {},
   "source": [
    "compare the true and estimated, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "371cdefc-f0ed-4163-8612-ac62f675814d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAADZCAYAAAA34mgzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABKfklEQVR4nO3deVxU1f8/8Bcgq2zhAposLuQW7olYuZRCaqbZYpt79tFwRTMpFZE+aWqaHzP9aCl9M9NMLTNSCLdU1EQpNSVF3EUrP4IrIJzfH/7mNjPcOzMXZ+f1fDzmwcy95577PnfGO2/PnHuuixBCgIiIiIjIgbjaOgAiIiIiIrWYxBIRERGRw2ESS0REREQOh0ksERERETkcJrFERERE5HCYxBIRERGRw2ESS0REREQOh0ksERERETmcarYOwJrKy8tx8eJF+Pn5wcXFxdbhEJGTEULg+vXrqFu3LlxdnbOPgOdRIrI0U8+lVSqJvXjxIkJDQ20dBhE5uXPnzqFevXq2DsMieB4lImsxdi6tUkmsn58fgHsHxd/f32r7LS0tRXp6OmJjY+Hu7m61/VqDs7bNWdsFOG/b7KFdRUVFCA0Nlc41zshW59HKsofPRVXBY20dVeE4m3ourVJJrOanL39/f6snsT4+PvD393e6D5yzts1Z2wU4b9vsqV3O/DO7rc6jlWVPnwtnx2NtHVXpOBs7lzrnoC0iIiIicmpMYomIiIjI4TCJJSIiIiKHwySWyMlFTP7B1iEQERGZXZW6sIuoKjo9q5dN9x8x+Qebx0BE8sz5n1xPN4HZ7YGHp29BcZnuBTk8B5AlsCeWiCyKX15ERGQJTGKJqAIOQSAiInvHJJaIKmDvKRER2TsmsURENrZ48WK0aNFCuoFATEwMfvzxR2n9nTt3EB8fjxo1asDX1xfPPfccLl++rFPH2bNn0atXL/j4+KB27dp46623cPfuXZ0y27dvR5s2beDp6YlGjRohNTXVGs0jIrIIJrFERDZWr149zJo1C9nZ2Thw4ACeeOIJ9OnTB0ePHgUAjB8/Ht9//z3Wrl2LHTt24OLFi+jXr5+0fVlZGXr16oWSkhLs2bMHn3/+OVJTUzFt2jSpTH5+Pnr16oWuXbsiJycH48aNw+uvv44tW7ZYvb1ERObAJJaITMJxspbTu3dv9OzZE5GRkXjooYfw73//G76+vti7dy8KCwvx2WefYd68eXjiiSfQtm1brFixAnv27MHevXsBAOnp6fj999+xcuVKtGrVCj169EBKSgoWLVqEkpISAMCSJUtQv359fPjhh2jatClGjRqF559/HvPnz7dl04mIKu2+p9gqKyvD4cOHER4ejgceeMAcMRGRHeI4WesoKyvD2rVrcfPmTcTExCA7OxulpaXo1q2bVKZJkyYICwtDVlYWOnTogKysLERFRSE4OFgqExcXh5EjR+Lo0aNo3bo1srKydOrQlBk3bpzBeIqLi1FcXCy9LioqAnDv/u2lpaVmaLFlaWJ0hFhtwdNNmK8uV6HzVxuPv/lUhc+0qW1TncSOGzcOUVFRGDZsGMrKytC5c2fs2bMHPj4+2LRpE7p06aK2SiKiKu/w4cOIiYnBnTt34Ovriw0bNqBZs2bIycmBh4cHAgMDdcoHBwejoKAAAFBQUKCTwGrWa9YZKlNUVITbt2/D29tbNq6ZM2ciOTm5wvL09HT4+PhUqq22kJGRYesQ7NLs9uavM6VdeYVlaWlp5t9RFefMn+lbt26ZVE51EvvNN9/gtddeAwB8//33yM/Px/Hjx/HFF1/g3Xffxe7du9VWSURU5TVu3Bg5OTkoLCzEN998g0GDBmHHjh22DguJiYlISEiQXhcVFSE0NBSxsbHw9/e3YWSmKS0tRUZGBrp37w53d3dbh2N3Hp5uvjHRnq4CKe3KMfWAK4rLdW92cGR6nNn2U9VVhc+05hcfY1QnsX/99RdCQkIA3Puf1QsvvICHHnoIQ4cOxYIFC9RWR0READw8PNCoUSMAQNu2bfHLL79gwYIF6N+/P0pKSnDt2jWd3tjLly9L5+KQkBDs379fpz7N7AXaZfRnNLh8+TL8/f0Ve2EBwNPTE56enhWWu7u7O9QXqKPFay36d9YyS53lLhXq5bE3P2f+TJvaLtUXdgUHB+P3339HWVkZNm/ejO7duwO41/Xr5uamtjoiIpJRXl6O4uJitG3bFu7u7sjMzJTW5ebm4uzZs4iJiQEAxMTE4PDhw7hy5YpUJiMjA/7+/mjWrJlURrsOTRlNHUREjkZ1T+yQIUPw4osvok6dOnBxcZEuFNi3bx+aNGli9gCJyLIiJv/Ai7ZsLDExET169EBYWBiuX7+OVatWYfv27diyZQsCAgIwbNgwJCQkICgoCP7+/hg9ejRiYmLQoUMHAEBsbCyaNWuGAQMGYPbs2SgoKMCUKVMQHx8v9aKOGDECH3/8MSZNmoShQ4di69at+Prrr/HDD5x1gogck+okdvr06Xj44Ydx7tw5vPDCC9IJ0s3NDZMnTzZ7gERkWfoJLJNa67ty5QoGDhyIS5cuISAgAC1atMCWLVukX7rmz58PV1dXPPfccyguLkZcXBw++eQTaXs3Nzds2rQJI0eORExMDKpXr45BgwZhxowZUpn69evjhx9+wPjx47FgwQLUq1cPn376KeLiOFaRiBxTpabYev755yssGzRo0H0HQ0TmVZmElAms9X322WcG13t5eWHRokVYtGiRYpnw8HCjV4B36dIFhw4dqlSMRET2plJJ7C+//IJt27bhypUrKC/XnUpj3rx5ZgmMiHQxISUiIvqH6iT2/fffx5QpU9C4cWMEBwfDxeWfKxC1nxOReSklpGqS24jJP+BESqw5wyIiIrIJ1UnsggULsHz5cgwePNgC4RCRWsYSWO0k9/SsXjp3QuH4VyIiclSqp9hydXXFo48+aolYiMiKIqemM4ElIiKHpTqJHT9+vMGLC4jIviglqhxWQEREjkz1cIKJEyeiV69eaNiwIZo1a1bhrgrr1683W3BEpJ7+EAH915FT07GA89sTEZGDU90TO2bMGGzbtg0PPfQQatSogYCAAJ2HpcycOROPPPII/Pz8ULt2bfTt2xe5ubkW2x+RLUVMNjwBvWa9sXLAPz2xmrLsgSUiImeguif2888/x7p169Crl3XH0u3YsQPx8fF45JFHcPfuXbzzzjuIjY3F77//jurVq1s1FiJLMzZWVftCLaV1+j2wase/GrroS+2MCBx7S0RE5qa6JzYoKAgNGza0RCwGbd68GYMHD0bz5s3RsmVLpKam4uzZs8jOzrZ6LES2ZErvK2A8aY2cml7p7dUkpUxgiYjIEip129mkpCSsWLECPj4+lojJJIWFhQDuJdVKiouLUVxcLL0uKioCAJSWlupMM2Rpmn1Zc5/W4qxts4d2RU5NV/zpXzsu/XKGttPf1lA9lWWuetSyh/fM2f4dEBHZMxchhFCzQevWrZGXlwchBCIiIipc2HXw4EGzBiinvLwczzzzDK5du4Zdu3Yplps+fTqSk5MrLF+1apVNE3AitcZmVcOCmLsGl2ley5Ul67h16xZeeeUVFBYWwt/f39bhWERRURECAgIcpo2lpaVIS0tDz549K3xfkem/7JjC001gdvsyTNrvhuIy69/8qKr86lMVPtOmnmdU98T27dv3fuIyi/j4eBw5csRgAgsAiYmJSEhIkF4XFRUhNDQUsbGxVj35lpaWIiMjA927d3e6D5yzts0W7VLqUVX6OzYrHT179pTKa55q/ir1iN5P22zVy2oKe/gsan7tISIiy1OdxCYlJVkiDpONGjUKmzZtws6dO1GvXj2DZT09PeHp6Vlhubu7u02+5Gy1X2tw1rZZs12nZ/WqcHct7Tjk1mt6UfTXKV1MpX3b2WYztlXYjykx2jtbfhad8d8AEZG9Un1hl60IITBq1Chs2LABW7duRf369W0dEpFE6Sc5Q8u112kno/oMJab6r439NKhd/kRKrE4iTERE5EhUJ7FlZWWYO3cu2rdvj5CQEAQFBek8LCU+Ph4rV67EqlWr4Ofnh4KCAhQUFOD27dsW2yeRqZR6KA0t108g1cz9KldGU6fSfpXqdYTeVSIiIn2qk9jk5GTMmzcP/fv3R2FhIRISEtCvXz+4urpi+vTpFgjxnsWLF6OwsBBdunRBnTp1pMeaNWsstk+i+yWXpOq/1vSw6ieg+tsq9brK9aZq9/RqnjNZJSIiZ6I6if3yyy+xbNkyTJgwAdWqVcPLL7+MTz/9FNOmTcPevXstESOAe8MJ5B6DBw+22D6JTKU0PEDutXbyqd0ja2jogf5tZOXqkqtbO0HmsAEiInImqpPYgoICREVFAQB8fX2l+Vqffvpp/PADvySpatFODOXujqX919D4VWO9pPpl9evT73mV21bNEAMiIiJ7p3p2gnr16uHSpUsICwtDw4YNkZ6ejjZt2uCXX36RnQmAyJkZGoOqoZ9EaveQ6vecKpWXW65dh7ExuUpl9GMxx21miYiIrEF1T+yzzz6LzMxMAMDo0aMxdepUREZGYuDAgRg6dKjZAySyJ0o/+5syO4FmW0PjWwH5XlvtdZrycsmwElNuIavmNrPswSUiIltT3RM7a9Ys6Xn//v0RHh6OPXv2IDIyEr179zZrcET2RqlnUy7JM3ShlvZ2+hd46c/3KldWfyiBKUmmZvt7N0qoBq37JCgy1INrSjkiIiJLUdUTW1paiqFDhyI/P19a1qFDByQkJDCBpSpFaUYA7fX66+SGDMiVNTR8QHvf+nXKxaj5q1/G1NvSmpqYMoElIiJrU5XEuru7Y926dZaKhchuRE5NV1xnrDfU2Hyv+jMSGEuG5WYc0F4ntz9DSaVc2zg8gIiIHI3qMbF9+/bFt99+a4FQiCzP1GRNc2tWOab8vK69XE2Ca6w3V27uWP0hCYbqPz2rlzScQDuZZU8qERE5GtVjYiMjIzFjxgzs3r0bbdu2RfXq1XXWjxkzxmzBEZmD3PhUtdtrtjVljKr2PuXGzhqa99XQMrn9Kg0pMNbeBTF3MTZL958/x7USEZEjUZ3EfvbZZwgMDER2djays7N11rm4uDCJJbtjSmKmn8BFTk3HgpiK2xsav6o/hZZc/cZmE5Cb8srQ1Fqa13LbG2ofERGRo1M9nCA/P1/xcerUKUvESGRx+j/3yw0n0B/Dqv1cqWfVWO+tXC+tpg65oQimXBAmlyTLJbAnUmINjseVaz8REZG9UJ3EEjkTpWQvcmq69HO7/kwA+lf9y9UhNwuB3ByxcmNc9depHQ5hqLfXUHJraKYDU/dNRERkLaqHEwDA+fPnsXHjRpw9exYlJSU66+bNm2eWwIgqy9hP54aSQs26EymxSEtLQ+TUdNm6lHpTtdfpk5tNQH/YgClzyRqq29DQB82y0tJSk+aJ5RAEIiKyZ6p7YjMzM9G4cWMsXrwYH374IbZt24YVK1Zg+fLlyMnJsUCIROoYS7zU9DYqDSvQlNUeYiDX06pdp9xQgMrMF6vUe2rsYjFtmnli1cx1S0REZE9UJ7GJiYmYOHEiDh8+DC8vL6xbtw7nzp1D586d8cILL1giRiKTqb3SX398q35Z/TlV5abLkqtT7ewDhuaLNTQ9l9zUWobikxtOwMSViIgckeok9tixYxg4cCAAoFq1arh9+zZ8fX0xY8YMfPDBB2YPkEgNtUmY/vhWpZ5YucRRuw5DQwruJyYlSkmvsfG6wD+JueYvE1ciInJEqpPY6tWrS+Ng69Spg7y8PGndX3/9Zb7IiFQyNF7UUBm5mwgAkG4IYGw8qmad3GwDctsZm7LL1ARYf8otY8MWtMf7Aso3dOAsBERE5AhUJ7EdOnTArl27AAA9e/bEhAkT8O9//xtDhw5Fhw4dzB4gkSH601wpXRylKauf1Mklr5oyC2LuKiZ6ctNZKdUnF4dSGbnZDIzRT9LlknZjMw8o1WcIk10iIrIl1bMTzJs3Dzdu3AAAJCcn48aNG1izZg0iIyM5MwFZndLFUNrUXNGvPU/s2KxqGJtV8Sd3Q8ml0sVYhuKoDP0hB0o9zHJDCgzdUrcyMRAREdmC6iS2QYMG0vPq1atjyZIlZg2I6H7oJ25K41VNTS5PpMTC3d29wk/0SkMM1Nw8QL8uJcYu0DKll1V7fWlpqWx9apNSTsFlPjNnzsT69etx/PhxeHt7o2PHjvjggw/QuHFjqcydO3cwYcIErF69GsXFxYiLi8Mnn3yC4OBgqczZs2cxcuRIbNu2Db6+vhg0aBBmzpyJatX+OdVv374dCQkJOHr0KEJDQzFlyhQMHjzYms0lIjKLSt/s4MCBA/jiiy/wxRdfVLj9LJE5qenBVOpdNTbsQD8Z056VQG1iqoapdd3PRV5KtNtYmWSUCaz57NixA/Hx8di7dy8yMjJQWlqK2NhY3Lx5Uyozfvx4fP/991i7di127NiBixcvol+/ftL6srIy9OrVCyUlJdizZw8+//xzpKamYtq0aVKZ/Px89OrVC127dkVOTg7GjRuH119/HVu2bLFqe4mIzEF1Env+/Hk8/vjjaN++PcaOHYuxY8fikUcewWOPPYbz589bIkaq4tQmcEpX9yv10MpdTKX5yV1/ii1rMrW3WGkMrPZ6/e3HZlWTve2sPo57tY7Nmzdj8ODBaN68OVq2bInU1FScPXtW6iAoLCzEZ599hnnz5uGJJ55A27ZtsWLFCuzZswd79+4FAKSnp+P333/HypUr0apVK/To0QMpKSlYtGiRdDHukiVLUL9+fXz44Ydo2rQpRo0aheeffx7z58+3WduJiCpL9XCC119/HaWlpTh27Jj0U1dubi6GDBmC119/HZs3bzZ7kFQ1GOs91L8wS/9ndLntleo0tC9zJW6mDBUwdVtD43eVepQNtVFzswP2ptqnwsJCAEBQUBAAIDs7G6WlpejWrZtUpkmTJggLC0NWVhY6dOiArKwsREVF6QwviIuLw8iRI3H06FG0bt0aWVlZOnVoyowbN04xluLiYhQXF0uvi4qKANwbliI3NMXeaGJ0hFhtwdNNmK8uV6Hz19qqyntcFT7TprZNdRK7Y8cO7NmzR2esVuPGjbFw4UI8/vjjaqsjkqgZ1yl3AZVcoqcpI5fcmTJ11v2oTL1KMywozaJgqC2GpvJSGw+TXespLy/HuHHj8Oijj+Lhhx8GABQUFMDDwwOBgYE6ZYODg1FQUCCV0U5gNes16wyVKSoqwu3bt+Ht7V0hnpkzZyI5ObnC8vT0dPj4+FSukTaQkZFh6xDs0uz25q8zpV25+Ss1QVpamk32ayvO/Jm+deuWSeVUJ7GhoaGyGXJZWRnq1q2rtjoiiSnJklwCp79Of9J/7XXa2+o/twdKMwooldXQHxphKGlXOzsBE1jrio+Px5EjR6SpDG0tMTERCQkJ0uuioiKEhoYiNjYW/v7+NozMNKWlpcjIyED37t3h7u5u63DszsPTzTce2tNVIKVdOaYecEVxuYvZ6jXVkelxVt+nLVSFz7TmFx9jVCexc+bMwejRo7Fo0SK0a9cOwL2LvMaOHYu5c+eqrY5IYkqypOZGAYDhK/ttkcBWdjYCcySSmjqM/UzDnlfbGTVqFDZt2oSdO3eiXr160vKQkBCUlJTg2rVrOr2xly9fRkhIiFRm//79OvVdvnxZWqf5q1mmXcbf31+2FxYAPD094enpWWG5u7u7Q32BOlq81lJcZv5ks7jcxSL1GlPV3l9n/kyb2i7VF3YNHjwYOTk5iI6Olk5u0dHROHjwIIYOHYqgoCDpQWSMsemjTNnW0B2qtJcZ2481Ejdj88saisNYvHI3WzB1ZgVDMzWQ5QkhMGrUKGzYsAFbt25F/fr1dda3bdsW7u7uyMzMlJbl5ubi7NmziImJAQDExMTg8OHDuHLlilQmIyMD/v7+aNasmVRGuw5NGU0dRESORHVP7EcffWSBMMhZaf+ULZcc6a9T+ilc7rl+/fr1ysWgv1/9WG1Fbsyrhtzxk2uffjml4y3XE6t0/Mk64uPjsWrVKnz33Xfw8/OTxrAGBATA29sbAQEBGDZsGBISEhAUFAR/f3+MHj0aMTEx0p0SY2Nj0axZMwwYMACzZ89GQUEBpkyZgvj4eKkndcSIEfj4448xadIkDB06FFu3bsXXX3+NH36wnyE1RESmUp3EDho0yBJxkJMylFDpl9GnlJgaGzeqf8GTvSWsSpSScbmxvmp6TvUT08ip6Vig0PHGBNY2Fi9eDADo0qWLzvIVK1ZINyKYP38+XF1d8dxzz+nc7EDDzc0NmzZtwsiRIxETE4Pq1atj0KBBmDFjhlSmfv36+OGHHzB+/HgsWLAA9erVw6effoq4uKoxlpCInIvqJJbIWgzNf2ookbXnBNVQ8q00LEKujNIMDIb+Q6BZdyIlFmlpaYicmm7ylFxkWUIYn5LIy8sLixYtwqJFixTLhIeHG71Cu0uXLjh06JDqGImI7A2TWLILpiRQckMPtP/K1WlovTVVJhalC7w0Y11NvUhM6bjqz1LABJbIvtjDuYvInlX6trNEahjrMTU0blaph1J/qIKhnkhHYWhsqqHeV6V2V+Z48IuTiIgcAZNYMitTrr6XGxqgVI/chV6G5od19ARM6QIvuYRWu1dWqe1ywxGM3UrX0f8jQEREVQOTWDIrtXO9ar/W/3nc0MVZpkyh5eiM9a7q98QqHT/9bdXe7ICIiMgemTQmtl+/fiZXuH79+koHY4pFixZhzpw5KCgoQMuWLbFw4UK0b2+B++bRfTM2zlXupgNyCa3c3beU5oN19J5YbfrJu37PrH5Ptv6wAmdN7omIiAATe2IDAgKkh7+/PzIzM3HgwAFpfXZ2NjIzMxEQEGCxQAFgzZo1SEhIQFJSEg4ePIiWLVsiLi5OZ3Jvsh9ywwD018uNZ5XrYdX/yVzNRP72Sm76Mf1EVC5xNXSclDjasSEiIjLGpCR2xYoV0iM4OBgvvvgi8vPzsX79eqxfvx6nTp3CSy+9hJo1a1o02Hnz5mH48OEYMmQImjVrhiVLlsDHxwfLly+36H5JHbnkUv/nbrmhAYYm/NfmLD2MSsdCm7H5dZXGF+sfT2c5ZkRERBqqp9havnw5du3aBTc3N2mZm5sbEhIS0LFjR8yZM8esAWqUlJQgOzsbiYmJ0jJXV1d069YNWVlZstsUFxejuLhYel1UVATg3r3jjd0/3pw0+7LmPi0tcmo6TqTEShPna7ftREqs4usTKbGyyZumLu3t9BM6TRln7lXUHB/t46F9zLTLaS7QkjvWcsdT85454+cRsI9/Z852TImI7JnqJPbu3bs4fvw4GjdurLP8+PHjKC8vN1tg+v766y+UlZUhODhYZ3lwcDCOHz8uu83MmTORnJxcYXl6ejp8fHwsEqchGRkZVt+npSyIAdLS0qQ7Pxlr29isalgQc1faVrPs3uu7/z8Zvvv/J2qX/1jqX1WvqU9TjyNZEHNXOiba8WuOqe7xqHjMtI+90uT2mvdI6bUzfR612bJdt27dstm+iYiqGtXf/kOGDMGwYcOQl5cnXVC1b98+zJo1C0OGDDF7gPcjMTERCQkJ0uuioiKEhoYiNjYW/v7+VoujtLQUGRkZ6N69O9zd3a22X2vQtG1sVjWDV7337Cm/LHJqus62xqZ/AqCqrL3SJKOatmvaonmtOTb6x1TuOGpo6jA2+4D257HZjG1OM1uBPfw70/zaQ0RElqc6iZ07dy5CQkLw4Ycf4tKlSwCAOnXq4K233sKECRPMHqBGzZo14ebmhsuXL+ssv3z5MkJCQmS38fT0hKenZ4Xl7u7uNvmSs9V+reFESmyFthmanUDpoiT9KbXkZhxwd3d3uCEFhsav6ifjmuNo6AYQ2uRu/GAKd3d3pxwra8t/Z87675uIyB6pnifW1dUVkyZNwoULF3Dt2jVcu3YNFy5cwKRJk3TGyZqbh4cH2rZti8zMTGlZeXk5MjMzERMTY7H9kjJNEqb0c76hKa+U5jXVT8jkkjhHS2AB+QvZDJWVu9DN2I0QnDEhJSIiUlKpmx3cvXsXP/30E7766iu4uLgAAC5evIgbN26YNTh9CQkJWLZsGT7//HMcO3YMI0eOxM2bN+1uGIOzMJYsapImzdhUue1MmSdW7rlcQuuMlNqnP5UW54AlIiLSpTqJPXPmDKKiotCnTx/Ex8fjzz//BAB88MEHmDhxotkD1Na/f3/MnTsX06ZNQ6tWrZCTk4PNmzdXuNiLzKOyiZKanlPt3la5/Rm6paq90k9M5dqln6AaOgb6t5fVMHZMHOmYERERqaU6iR07dizatWuH//3vf/D29paWP/vsszo/9VvKqFGjcObMGRQXF2Pfvn2Ijo62+D6rOnMkQ3J1KN2Jy1DSZ8/k5mpVGlJh6Paw2kmr9jKl2/XeLya7RETkiFQnsT///DOmTJkCDw8PneURERG4cOGC2QIj+2Hs4iwl+rdD1b/jluav0gVP+nXZeyJrLFmV66HVPzb6x0wpcTUl8TT1eNn7cSUiIpKjOoktLy9HWVlZheXnz5+Hn5+fWYIi+6afdBq6sEu/vNJP7IaGDahJ3CzNlIRPLgFVmm1BLonXX27omBAREVVVqpPY2NhYfPTRR9JrFxcX3LhxA0lJSehpaBJLchr6CZT+hV2AclKq3wNr7MItQzMc2IKh28JqJ9v6PbBKCal+4qrfS6tfr7F4iIiIqgrV88R++OGHiIuLQ7NmzXDnzh288sorOHHiBGrWrImvvvrKEjGSAzElMTUlGdP+id3eEllt2r2o+s/lxsbKJbr387M/e2SJiKiqUt0TW69ePfz666949913MX78eLRu3RqzZs3CoUOHULt2bUvESHZGKaE0lsDKDStQunjJlKvw7SmB009glWZo0E9c7akNREREjkR1T+zOnTvRsWNHvPrqq3j11Vel5Xfv3sXOnTvRqVMnswZIjkNu2iilMtqvK9vLaq3eWWMxyg0XMDRXrtwsDI5w4RoREZE9Ud0T27VrV1y9erXC8sLCQnTt2tUsQZF9M5Zs6f+0rn8HLkA+2bPX+WANjWPVXq80PZh+HUpJLREREZlOdRIrhJDu0qXt77//RvXq1c0SFDkGTWKmPTuBsdupmpIA2zulIQFKc7vqj4VVqkd/PRERESkzOYnt168f+vXrBxcXFwwePFh63a9fP/Tp0wdxcXHo2LGjJWMlG1PqkZSbncDQnacMlTG23hzUXiAlN35VacosuXG9hmZeUBMHERER/cPkJDYgIAABAQEQQsDPz096HRAQgJCQELzxxhtYuXKlJWMlG1P6GV1unlili5tMHSurXb8pZdVQmgnBUHm5+V21/8rNC8sklYiIyHJMvrBrxYoVAO7dmWvixIkcOlCFKA0DOD2rF0pLS3V6YpWSUf1pp4z1fBpLKu+HUkKq9iIz/WEEcr2yHO9KRERkGarHxCYlJTGBrWKUElNTEj6lJE7/ZgdK+1GK534SQ0NThBnaRmlaLP2L15TuTKZ/AwhTphEjIiIieaqTWAD45ptv8OKLL6JDhw5o06aNzoOqBu1ETTOcwFAvpLEkTc2tVc2R8MklonJOpMTKXqClNjbtabWUxtgSERGR6VQnsf/5z38wZMgQBAcH49ChQ2jfvj1q1KiBU6dOoUePHpaIkeycZjiB3FX2SlNt6ZeT2+Z+hxQYSgxN+dl/QcxdRE5NN7itfu+q9rAJa1ykRkREVFWpTmI/+eQTLF26FAsXLoSHhwcmTZqEjIwMjBkzBoWFhZaIkSzgfnozTU0g9ceMmpKcmoux8ajGLiQ7kRKLsVnVcCIlVrZO/eRcrleViSsREZHlqE5iz549K02l5e3tjevXrwMABgwYgK+++sq80ZHF3E+CpX+RFgCpx1Lpgin9MaNyMRhKOE2JV+4OWHIXksnN76q/PnJqus4Fa8ZuF2vshghK5YiIiKhyVN92NiQkBFevXkV4eDjCwsKwd+9etGzZEvn5+RBCWCJGskOaJK20tBQApB5LQ3Ok6ieNhqbQMuWiJ7lySvXKzZCgX4f+tmOzqmFsVrrs8AdTx8jKxUxEROZnrU4Cnsfth+qe2CeeeAIbN24EAAwZMgTjx49H9+7d0b9/fzz77LNmD5AsR82YU1OGA+iPDZXrrZSrw1ivpn45/WTY0Bhb/WENcj3Bcj22J1JisSDmboXk3FBPsiltICIiIvNQncQuXboU7777LgAgPj4ey5cvR9OmTTFjxgwsXrzY7AGS5RhKHg1N1G/qHK/af+USTqWZC5SGD2jXpz1jQGUSYu31+rMoREz+QRoeYezCLiVq/6fOpJeIiEgd1cMJXF1d4er6T+770ksv4aWXXjJrUGQbpiSncmXGZlVDz54VyyrVqzQWVamn19BYV31y2yrdYEGpPdrbaF/Ypb/enD8p8ecpIiIidSo1T+ydO3ewf/9+bNq0CRs3btR5kPNS+tneWHn9Hlm5ZUqJplwZNWNR5eaqlZsOSz9mueRVrk4ic9m5cyd69+6NunXrwsXFBd9++63OeiEEpk2bhjp16sDb2xvdunXDiRMndMpcvXoVr776Kvz9/REYGIhhw4bhxo0bOmV+++03PP744/Dy8kJoaChmz55t6aYREVmE6iR28+bNCAsLQ4cOHfDMM8+gb9++0oNjYh1LZW5AoFmuncTJzROrdBGXoWEK2vs0lKTKjWWVW2doVgS5WI21W6k8kTncvHkTLVu2xKJFi2TXz549G//5z3+wZMkS7Nu3D9WrV0dcXBzu3LkjlXn11Vdx9OhRZGRkYNOmTdi5cyfeeOMNaX1RURFiY2MRHh6O7OxszJkzB9OnT8fSpUst3j4iInNTncSOHj0aL7zwAi5duoTy8nKdR1lZmSViJAsxlowpJXemJHFK41wNzQ6gXbdcb6tST7DccAO5eOR6f40l1UTW0qNHD7z33nuynQFCCHz00UeYMmUK+vTpgxYtWuD//u//cPHiRanH9tixY9i8eTM+/fRTREdH47HHHsPChQuxevVqXLx4EQDw5ZdfoqSkBMuXL0fz5s3x0ksvYcyYMZg3b541m0pEZBaqx8RevnwZCQkJCA4OtkQ8ZCfkErzKJnfGxrFq6pabEUGpHkNzympva2xOWKVlQMWxvsZwmAFZSn5+PgoKCtCtWzdpWUBAAKKjo5GVlYWXXnoJWVlZCAwMRLt27aQy3bp1g6urK/bt24dnn30WWVlZ6NSpEzw8PKQycXFx+OCDD/C///0PDzzwQIV9FxcXo7i4WHpdVFQE4N70epop9uyZJkZHiFWfp5tjTVvp6Sp0/jorW3+WHPkzbSpT26Y6iX3++eexfft2NGzYUHVQ5DiUejMNJYNy25uSCCv1pOonopW9QMvQ0AGlZQB0bnZgCiawZCkFBQUAUKHzIDg4WFpXUFCA2rVr66yvVq0agoKCdMrUr1+/Qh2adXJJ7MyZM5GcnFxheXp6Onx8fCrZIuvLyMiwdQiqzW5v6wgqJ6Vdua1DsKi0tDRbhwDAMT/Tprp165ZJ5VQnsR9//DFeeOEF/Pzzz4iKioK7u7vO+jFjxqitkuycoemxlMrKkUuA5YYW6E+hpV3O2JhZNVNVseeUyLjExEQkJCRIr4uKihAaGorY2Fj4+/vbMDLTlJaWIiMjA927d6/wfWXvHp6+xdYhqOLpKpDSrhxTD7iiuNzF1uFYzJHpcTbdvyN/pk2l+cXHGNVJ7FdffYX09HR4eXlh+/btcHH554Pq4uLCJNYJGesJlSurT/8mCEq9rNrl5W4sYGw4gJqk1FDZe7edNbkqIosKCQkBcG84V506daTlly9fRqtWraQyV65c0dnu7t27uHr1qrR9SEgILl++rFNG81pTRp+npyc8PT0rLHd3d3eoL1BHixcAisscMxEsLndx2NhNYS+fI0f8TJvK1HapvrDr3XffRXJyMgoLC3H69Gnk5+dLj1OnTqkOlOyPoTtm6Sd+kVPTMTarmsHtNNtqlzF2owO59UpjayuTVGvXSWTv6tevj5CQEGRmZkrLioqKsG/fPsTE3PvfVkxMDK5du4bs7GypzNatW1FeXo7o6GipzM6dO3XGm2VkZKBx48ayQwmIiOyZ6iS2pKQE/fv317nhATkXU3szIyb/IN2eVX+59l/t5/pDEvTHrcrNG2vsrlz3Q6kuQ/PEElnCjRs3kJOTg5ycHAD3LubKycnB2bNn4eLignHjxuG9997Dxo0bcfjwYQwcOBB169ZF3759AQBNmzbFU089heHDh2P//v3YvXs3Ro0ahZdeegl169YFALzyyivw8PDAsGHDcPToUaxZswYLFizQGS5AROQoVGeigwYNwpo1aywRCzkQTeKpfVtWuemy5Kaykuv9NDabgP4wBLm6lO74ReQIDhw4gNatW6N169YAgISEBLRu3RrTpk0DAEyaNAmjR4/GG2+8gUceeQQ3btzA5s2b4eXlJdXx5ZdfokmTJnjyySfRs2dPPPbYYzpzwAYEBCA9PR35+flo27YtJkyYgGnTpunMJUtE5ChUj4ktKyvD7NmzsWXLFrRo0aLCuAXON1g1aN/ZKi0tDZFT02XHpSpdwKVUn+a5qVNhaZc1ZfysPl7cRfaiS5cuEEJ5aiIXFxfMmDEDM2bMUCwTFBSEVatWGdxPixYt8PPPP1c6TiIie6E6iT18+LDUU3DkyBGdddoXeZH9M2cCp/n5XSmBVJqyS2kbuWEGcheEKdVtaLnaMkRERGR/VCex27Zts0QcBp0+fRopKSnYunUrCgoKULduXbz22mt49913dSbtJnXMncAp3YxAf3yr3P6N3VyhMrMPEBERkfNyiKuzjh8/jvLycvz3v//F0aNHMX/+fCxZsgTvvPOOrUOrcpTGneoPJ1BKVuVmONCvX2ncK8e6EhERkYZJPbH9+vVDamoq/P390a9fP4Nl169fb5bAtD311FN46qmnpNcNGjRAbm4uFi9ejLlz55p9f6RMabyqZjiB/vhUuRslaOrRXq//XK6smrGuRERE5NxM6okNCAiQxrv6+/sjICBA8WEthYWFCAoKstr+6B/6CaRmnljNOrleWLneVTXjZi01xRYRERE5JpN6YlesWCE9T01NtVQsJjt58iQWLlxotBe2uLgYxcXF0mvNbcxKS0t1Jvu2NM2+rLlPa4icmo7fp3XFgpi7KC0tReTUdJxIicWJlFjpNQCdZZrtND23mmXa67WdSImV5qOV2177uTk563sGOG/b7KFdznZMiYjsmYswNKeLjCeeeALr169HYGCgzvKioiL07dsXW7duNbmuyZMn44MPPjBY5tixY2jSpIn0+sKFC+jcuTO6dOmCTz/91OC206dPR3JycoXlq1atgo+Pj8lxknpjs6pJN0HQfi73mshZ3Lp1C6+88goKCwvh7+9v63AsoqioCAEBAQ7TxtLSUqSlpaFnz54Od4tOR7sOwNNNYHb7Mkza7+bUt5219a+BjvyZNpWp5xnVsxNs374dJSUlFZbfuXNH9dyDEyZMwODBgw2WadCggfT84sWL6Nq1Kzp27KgzgbeSxMREnTvRFBUVITQ0FLGxsVY9+ZaWliIjIwPdu3d3ug9c5NR0LIi5W6FtPXtC57l2j+nYrHT01Cpgqd7U++HM75mzts0e2qX5tYeIiCzP5CT2t99+k57//vvvKCgokF6XlZVh8+bNePDBB1XtvFatWqhVq5ZJZS9cuICuXbuibdu2WLFihUm3vfX09ISnp2eF5e7u7jb5krPVfi1Jc7MDY20zdDMDY1Nt2ZIzvmcazto2W7bLGY8nEZG9MjmJbdWqFVxcXODi4oInnniiwnpvb28sXLjQrMFpXLhwAV26dEF4eDjmzp2LP//8U1oXEhJikX1WVcaSSHMlmUoXcZlSt70lukRERGR9Jiex+fn5EEKgQYMG2L9/v04PqoeHB2rXrg03NzeLBJmRkYGTJ0/i5MmTqFevns46lUN6yQhjyaH++nvDCYzXa0riaWpiygSWiIiITL7ZQXh4OCIiIlBeXo527dohPDxcetSpU8diCSwADB48GEII2QdZh9IFBqaOZWXiSUREROak+sKuzz//HDVr1kSvXveSkkmTJmHp0qVo1qwZvvrqK4SHh5s9SLI+U+ZwJSKqihxt1gAiZ6X6trPvv/8+vL29AQBZWVn4+OOPMXv2bNSsWRPjx483e4BkG/pJK0/aREREZE9UJ7Hnzp1Do0aNAADffvstnn/+ebzxxhuYOXOm6im2yHGo6YllwktERESWpno4ga+vL/7++2+EhYUhPT1dmofVy8sLt2/fNnuA5Dg0ySuHHhAREZGlqU5iu3fvjtdffx2tW7fGH3/8IU1af/ToUURERJg7PrIDpianTF6JiIjIWlQPJ1i0aBFiYmLw559/Yt26dahRowYAIDs7Gy+//LLZAyTbOz2rl8kJKocSEBERkTWo7okNDAzExx9/XGF5cnKyWQIix8beWCIiIrIG1T2xAPDzzz/jtddeQ8eOHXHhwgUAwBdffIFdu3aZNTiyb+x1JSIiIltRncSuW7cOcXFx8Pb2xsGDB1FcXAwAKCwsxPvvv2/2AMl+sdeViIiIbEV1Evvee+9hyZIlWLZsGdzd3aXljz76KA4ePGjW4Mj+RU5Nt3UIREREVAWpTmJzc3PRqVOnCssDAgJw7do1c8REDsTU284SERERmZPqJDYkJAQnT56ssHzXrl1o0KCBWYIiIiIiIjJE9ewEw4cPx9ixY7F8+XK4uLjg4sWLyMrKwsSJEzF16lRLxEhERERkF6x1UTOvOzFOdRI7efJklJeX48knn8StW7fQqVMneHp6YuLEiRg9erQlYiQiIiIi0qF6OIGLiwveffddXL16FUeOHMHevXvx559/IiUlxRLxkQPjFFxERERkKZWaJxYAPDw80KxZM7Rv3x6+vr7mjInMyJaJJH8KISIiIkupdBJLjoGJJBERETkjJrFERERE5HCYxBIRERGRw2ESS0REREQOh0lsFcRZA4iIiMjRqZ4nlhwfL/YiImck9x90TzeB2e2Bh6dvQXGZiw2iIiJLYU8sERERETkcJrFERERE5HCYxBIRVTGLFi1CREQEvLy8EB0djf3799s6JCIi1TgmloioClmzZg0SEhKwZMkSREdH46OPPkJcXBxyc3NRu3ZtW4dHRP+f0kXY5hzn7ejXyLAnloioCpk3bx6GDx+OIUOGoFmzZliyZAl8fHywfPlyW4dGRKRKleqJFUIAAIqKiqy639LSUty6dQtFRUVwd3e36r4tzVnb5qztApy3bfbQLs25RXOusTclJSXIzs5GYmKitMzV1RXdunVDVlaW7DbFxcUoLi6WXhcWFgIArl69itLSUpP3HT0zs5JRm07uC61aucCtW+WoVuqKsnLOTmBJPNbWYc7j3Gji12aKyrB9iU+qKn/9+nUAxs+lVSqJ1RyU0NBQG0dCRM7s+vXrCAgIsHUYFfz1118oKytDcHCwzvLg4GAcP35cdpuZM2ciOTm5wvL69etbJEZLeMXWAVQhPNbW4WjHueaHldvO2Lm0SiWxdevWxblz5+Dn5wcXF+v9L7GoqAihoaE4d+4c/P39rbZfa3DWtjlruwDnbZs9tEsIgevXr6Nu3bo22b8lJCYmIiEhQXpdXl6Oq1evokaNGlY9j1aWPXwuqgoea+uoCsfZ1HNplUpiXV1dUa9ePZvt39/f32k/cM7aNmdtF+C8bbN1u+yxB1ajZs2acHNzw+XLl3WWX758GSEhIbLbeHp6wtPTU2dZYGCgpUK0GFt/LqoSHmvrcPbjbMq5lBd2ERFVER4eHmjbti0yM/8Zn1peXo7MzEzExMTYMDIiIvWqVE8sEVFVl5CQgEGDBqFdu3Zo3749PvroI9y8eRNDhgyxdWhERKowibUCT09PJCUlVfhJzhk4a9uctV2A87bNWdtlbv3798eff/6JadOmoaCgAK1atcLmzZsrXOzlLPi5sB4ea+vgcf6Hi7DXuWCIiIiIiBRwTCwRERERORwmsURERETkcJjEEhEREZHDYRJLRERERA6HSayFREREwMXFRecxa9Ysg9vcuXMH8fHxqFGjBnx9ffHcc89VmJTclk6fPo1hw4ahfv368Pb2RsOGDZGUlISSkhKD23Xp0qXCsRgxYoSVola2aNEiREREwMvLC9HR0di/f7/B8mvXrkWTJk3g5eWFqKgopKWlWSlS082cOROPPPII/Pz8ULt2bfTt2xe5ubkGt0lNTa3w/nh5eVkpYtNNnz69QpxNmjQxuI0jvGdkPZU9h5Fxas+npF5lzu/OjkmsBc2YMQOXLl2SHqNHjzZYfvz48fj++++xdu1a7NixAxcvXkS/fv2sFK1xx48fR3l5Of773//i6NGjmD9/PpYsWYJ33nnH6LbDhw/XORazZ8+2QsTK1qxZg4SEBCQlJeHgwYNo2bIl4uLicOXKFdnye/bswcsvv4xhw4bh0KFD6Nu3L/r27YsjR45YOXLDduzYgfj4eOzduxcZGRkoLS1FbGwsbt68aXA7f39/nffnzJkzVopYnebNm+vEuWvXLsWyjvKekfXczzmMlKk9n1LlVPb87tQEWUR4eLiYP3++yeWvXbsm3N3dxdq1a6Vlx44dEwBEVlaWBSI0j9mzZ4v69esbLNO5c2cxduxY6wRkovbt24v4+HjpdVlZmahbt66YOXOmbPkXX3xR9OrVS2dZdHS0+Ne//mXROO/XlStXBACxY8cOxTIrVqwQAQEB1guqkpKSkkTLli1NLu+o7xlZlynnMDJM7fmUzMOU87uzY0+sBc2aNQs1atRA69atMWfOHNy9e1exbHZ2NkpLS9GtWzdpWZMmTRAWFoasrCxrhFsphYWFCAoKMlruyy+/RM2aNfHwww8jMTERt27dskJ08kpKSpCdna1zrF1dXdGtWzfFY52VlaVTHgDi4uLs+r0B7r0/AIy+Rzdu3EB4eDhCQ0PRp08fHD161BrhqXbixAnUrVsXDRo0wKuvvoqzZ88qlnXU94ysy9RzGMmrzPmUzMPU87sz4x27LGTMmDFo06YNgoKCsGfPHiQmJuLSpUuYN2+ebPmCggJ4eHggMDBQZ3lwcDAKCgqsELF6J0+exMKFCzF37lyD5V555RWEh4ejbt26+O233/D2228jNzcX69evt1Kkuv766y+UlZVVuENRcHAwjh8/LrtNQUGBbHl7fW8AoLy8HOPGjcOjjz6Khx9+WLFc48aNsXz5crRo0QKFhYWYO3cuOnbsiKNHj6JevXpWjNiw6OhopKamonHjxrh06RKSk5Px+OOP48iRI/Dz86tQ3hHfM7IuU89hpKwy51O6f6ae350dk1gVJk+ejA8++MBgmWPHjqFJkyZISEiQlrVo0QIeHh7417/+hZkzZ9rdreLUtEvjwoULeOqpp/DCCy9g+PDhBrd94403pOdRUVGoU6cOnnzySeTl5aFhw4b3Fzwpio+Px5EjRwyOGwWAmJgYxMTESK87duyIpk2b4r///S9SUlIsHabJevToIT1v0aIFoqOjER4ejq+//hrDhg2zYWRka5Y+hxHZG1PP786OSawKEyZMwODBgw2WadCggezy6Oho3L17F6dPn0bjxo0rrA8JCUFJSQmuXbum0xt7+fJlhISE3E/YRqlt18WLF9G1a1d07NgRS5cuVb2/6OhoAPd6QWyRxNasWRNubm4VZn4wdKxDQkJUlbe1UaNGYdOmTdi5c6fq3lR3d3e0bt0aJ0+etFB05hEYGIiHHnpIMU5He8+o8qx9DqN/VOZ8Svfnfs7vzoZJrAq1atVCrVq1KrVtTk4OXF1dUbt2bdn1bdu2hbu7OzIzM/Hcc88BAHJzc3H27FmdXjJLUNOuCxcuoGvXrmjbti1WrFgBV1f1w6pzcnIAAHXq1FG9rTl4eHigbdu2yMzMRN++fQHc+2kmMzMTo0aNkt0mJiYGmZmZGDdunLQsIyPD4u+NWkIIjB49Ghs2bMD27dtRv3591XWUlZXh8OHD6NmzpwUiNJ8bN24gLy8PAwYMkF3vKO8Z3T9rn8PoH5U5n1LlmOP87nRsfWWZM9qzZ4+YP3++yMnJEXl5eWLlypWiVq1aYuDAgVKZ8+fPi8aNG4t9+/ZJy0aMGCHCwsLE1q1bxYEDB0RMTIyIiYmxRRNknT9/XjRq1Eg8+eST4vz58+LSpUvSQ7uMdrtOnjwpZsyYIQ4cOCDy8/PFd999Jxo0aCA6depkq2YIIYRYvXq18PT0FKmpqeL3338Xb7zxhggMDBQFBQVCCCEGDBggJk+eLJXfvXu3qFatmpg7d644duyYSEpKEu7u7uLw4cO2aoKskSNHioCAALF9+3ad9+fWrVtSGf22JScniy1btoi8vDyRnZ0tXnrpJeHl5SWOHj1qiyYomjBhgti+fbvIz88Xu3fvFt26dRM1a9YUV65cEUI47ntG1mPKOYzUM3Y+JfMw5fxe1TCJtYDs7GwRHR0tAgIChJeXl2jatKl4//33xZ07d6Qy+fn5AoDYtm2btOz27dvizTffFA888IDw8fERzz77rF2dXFesWCEAyD409Nt19uxZ0alTJxEUFCQ8PT1Fo0aNxFtvvSUKCwtt1Ip/LFy4UISFhQkPDw/Rvn17sXfvXmld586dxaBBg3TKf/311+Khhx4SHh4eonnz5uKHH36wcsTGKb0/K1askMrot23cuHHScQgODhY9e/YUBw8etH7wRvTv31/UqVNHeHh4iAcffFD0799fnDx5UlrvqO8ZWY8p5zCqHEPnUzIPU87vVY2LEEJYq9eXiIiIiMgcOBiIiIiIiBwOk1giIiIicjhMYomIiIjI4TCJJSIiIiKHwySWiIiIiBwOk1giIiIicjhMYomIiIjI4TCJJSIiogq6dOmic9tmInvDJJaIiIhUmz59Olq1amXrMKgKYxJLpEUIgbt379o6DCIisygpKbF1CCgtLbV1CKo5YsxVEZNYcmhdunTBqFGjMGrUKAQEBKBmzZqYOnUqNHdT/uKLL9CuXTv4+fkhJCQEr7zyCq5cuSJtv337dri4uODHH39E27Zt4enpiV27diEvLw99+vRBcHAwfH198cgjj+Cnn37S2XdERATee+89DBw4EL6+vggPD8fGjRvx559/ok+fPvD19UWLFi1w4MABaZszZ86gd+/eeOCBB1C9enU0b94caWlp1jlYRGR2p0+fhouLS4VHly5dZMsLITB9+nSEhYXB09MTdevWxZgxY6T1xcXFePvttxEaGgpPT080atQIn332mbR+x44daN++PTw9PVGnTh1MnjxZ5z/emnPiuHHjULNmTcTFxQEAjhw5gh49esDX1xfBwcEYMGAA/vrrL2m7mzdvSueyOnXq4MMPPzTY7tTUVCQnJ+PXX3+V2pyamgoAcHFxweLFi/HMM8+gevXq+Pe//43U1FQEBgbq1PHtt9/CxcVFZ9l3332HNm3awMvLCw0aNEBycrLBjoXBgwejb9++SE5ORq1ateDv748RI0boJO+bN2/GY489hsDAQNSoUQNPP/008vLypPWa93DNmjXo3LkzvLy88OWXX+Lvv//Gyy+/jAcffBA+Pj6IiorCV199pbP/Ll26YPTo0Rg3bhweeOABBAcHY9myZbh58yaGDBkCPz8/NGrUCD/++KO0zf/+9z+8+uqrqFWrFry9vREZGYkVK1YYPN6kQBA5sM6dOwtfX18xduxYcfz4cbFy5Urh4+Mjli5dKoQQ4rPPPhNpaWkiLy9PZGVliZiYGNGjRw9p+23btgkAokWLFiI9PV2cPHlS/P333yInJ0csWbJEHD58WPzxxx9iypQpwsvLS5w5c0baNjw8XAQFBYklS5aIP/74Q4wcOVL4+/uLp556Snz99dciNzdX9O3bVzRt2lSUl5cLIYTo1auX6N69u/jtt99EXl6e+P7778WOHTuse9CIyGzu3r0rLl26JD0OHTokatSoIaZOnSpbfu3atcLf31+kpaWJM2fOiH379knnKyGEePHFF0VoaKhYv369yMvLEz/99JNYvXq1EEKI8+fPCx8fH/Hmm2+KY8eOiQ0bNoiaNWuKpKQkaXvNOfGtt94Sx48fF8ePHxf/+9//RK1atURiYqI4duyYOHjwoOjevbvo2rWrtN3IkSNFWFiY+Omnn8Rvv/0mnn76aeHn5yfGjh0r245bt26JCRMmiObNm0ttv3XrlhBCCACidu3aYvny5SIvL0+cOXNGrFixQgQEBOjUsWHDBqGdhuzcuVP4+/uL1NRUkZeXJ9LT00VERISYPn264vEfNGiQ8PX1Ff379xdHjhwRmzZtErVq1RLvvPOOVOabb74R69atEydOnBCHDh0SvXv3FlFRUaKsrEwIIUR+fr4AICIiIsS6devEqVOnxMWLF8X58+fFnDlzxKFDh0ReXp74z3/+I9zc3MS+fft0jrefn59ISUkRf/zxh0hJSRFubm6iR48eYunSpdJ3Q40aNcTNmzeFEELEx8eLVq1aiV9++UXk5+eLjIwMsXHjRsU2kjImseTQOnfurJMkCiHE22+/LZo2bSpb/pdffhEAxPXr14UQ/ySx3377rdF9NW/eXCxcuFB6HR4eLl577TXp9aVLlwQAnS+vrKwsAUBcunRJCCFEVFSUwRMyETmu27dvi+joaPH0009LCZK+Dz/8UDz00EOipKSkwrrc3FwBQGRkZMhu+84774jGjRvrnO8WLVokfH19pf117txZtG7dWme7lJQUERsbq7Ps3LlzAoDIzc0V169fFx4eHuLrr7+W1v/999/C29tbMYkVQoikpCTRsmXLCssBiHHjxuksMyWJffLJJ8X777+vU+aLL74QderUUYxh0KBBIigoSEoQhRBi8eLFOsdE359//ikAiMOHDwsh/kliP/roI8X9aPTq1UtMmDBBet25c2fx2GOPSa/v3r0rqlevLgYMGCAt03w3ZGVlCSGE6N27txgyZIjRfZFxHE5ADq9Dhw46P0nFxMTgxIkTKCsrQ3Z2Nnr37o2wsDD4+fmhc+fOAICzZ8/q1NGuXTud1zdu3MDEiRPRtGlTBAYGwtfXF8eOHauwXYsWLaTnwcHBAICoqKgKyzRDGMaMGYP33nsPjz76KJKSkvDbb7/db/OJyE4MHToU169fx6pVq+DqKv/1+sILL+D27dto0KABhg8fjg0bNkg/l+fk5MDNzU06T+k7duwYYmJidM53jz76KG7cuIHz589Ly9q2bauz3a+//opt27bB19dXejRp0gQAkJeXh7y8PJSUlCA6OlraJigoCI0bN67cgUDFc6opfv31V8yYMUMnzuHDh+PSpUu4deuW4nYtW7aEj4+P9DomJgY3btzAuXPnAAAnTpzAyy+/jAYNGsDf3x8REREAjH8PlJWVISUlBVFRUQgKCoKvry+2bNli8HvAzc0NNWrUMPg9MHLkSKxevRqtWrXCpEmTsGfPHlMPEelhEktO686dO4iLi4O/vz++/PJL/PLLL9iwYQOAihc7VK9eXef1xIkTsWHDBrz//vv4+eefkZOTg6ioqArbubu7S881Xyxyy8rLywEAr7/+Ok6dOoUBAwbg8OHDaNeuHRYuXGimFhORrbz33nvYsmULNm7cCD8/P8VyoaGhyM3NxSeffAJvb2+8+eab6NSpE0pLS+Ht7W2WWPTPZzdu3EDv3r2Rk5Oj8zhx4gQ6depkln0ai8HV1VW6VkFD/+KpGzduIDk5WSfGw4cP48SJE/Dy8qp0LL1798bVq1exbNky7Nu3D/v27QNg/Htgzpw5WLBgAd5++21s27YNOTk5iIuLM/g9ANw77xv6HujRowfOnDmD8ePH4+LFi3jyyScxceLESrevKqtm6wCI7pfmhKSxd+9eREZG4vjx4/j7778xa9YshIaGAoDORVaG7N69G4MHD8azzz4L4N7J9fTp02aJNzQ0FCNGjMCIESOQmJiIZcuWYfTo0Wapm4isb926dZgxYwZ+/PFHNGzY0Gh5b29v9O7dG71790Z8fDyaNGmCw4cPIyoqCuXl5dixYwe6detWYbumTZti3bp1EEJIidHu3bvh5+eHevXqKe6vTZs2WLduHSIiIlCtWsWv/YYNG8Ld3R379u1DWFgYgHsXH/3xxx+KvcIA4OHhgbKyMqPtBYBatWrh+vXruHnzppQs5uTkVIgzNzcXjRo1MqlOjV9//RW3b9+W/hOwd+9e+Pr6IjQ0FH///Tdyc3OxbNkyPP744wCAXbt2mVTv7t270adPH7z22msA7iWhf/zxB5o1a6YqPjm1atXCoEGDMGjQIDz++ON46623MHfu3Puut6phTyw5vLNnzyIhIQG5ubn46quvsHDhQowdOxZhYWHw8PDAwoULcerUKWzcuBEpKSkm1RkZGYn169cjJycHv/76K1555RXpf9H3Y9y4cdiyZQvy8/Nx8OBBbNu2DU2bNr3veonINo4cOYKBAwfi7bffRvPmzVFQUICCggJcvXpVtnxqaio+++wzHDlyBKdOncLKlSvh7e2N8PBwREREYNCgQRg6dCi+/fZb5OfnY/v27fj6668BAG+++SbOnTuH0aNH4/jx4/juu++QlJSEhIQExeELABAfH4+rV6/i5Zdfxi+//IK8vDxs2bIFQ4YMQVlZGXx9fTFs2DC89dZb2Lp1K44cOYLBgwcbrBO4N0NLfn4+cnJy8Ndff6G4uFixbHR0NHx8fPDOO+8gLy8Pq1atkmYz0Jg2bRr+7//+D8nJyTh69CiOHTuG1atXY8qUKQbjKCkpwbBhw/D7778jLS0NSUlJGDVqFFxdXfHAAw+gRo0aWLp0KU6ePImtW7ciISHBYH0akZGRyMjIwJ49e3Ds2DH861//wuXLl03a1pBp06bhu+++w8mTJ3H06FFs2rSJ3wOVxCSWHN7AgQNx+/ZttG/fHvHx8Rg7dizeeOMN1KpVC6mpqVi7di2aNWuGWbNmmfw/3Xnz5uGBBx5Ax44d0bt3b8TFxaFNmzb3HWtZWRni4+PRtGlTPPXUU3jooYfwySef3He9RGQbBw4cwK1bt/Dee++hTp060qNfv36y5QMDA7Fs2TI8+uijaNGiBX766Sd8//33qFGjBgBg8eLFeP755/Hmm2+iSZMmGD58OG7evAkAePDBB5GWlob9+/ejZcuWGDFiBIYNG2Y0yatbty52796NsrIyxMbGIioqCuPGjUNgYKCUqM6ZMwePP/44evfujW7duuGxxx6rMLZW33PPPYennnoKXbt2Ra1atSpMP6UtKCgIK1euRFpamjRV1fTp03XKxMXFYdOmTUhPT8cjjzyCDh06YP78+QgPDzcYx5NPPonIyEh06tQJ/fv3xzPPPCPV7erqitWrVyM7OxsPP/wwxo8fjzlz5hisT2PKlClo06YN4uLi0KVLF4SEhKBv374mbWuIh4cHEhMT0aJFC3Tq1Alubm5YvXr1fddbFbkI/UEqRA6kS5cuaNWqFT766CNbh0JERFY2ePBgXLt2Dd9++62tQyEbYE8sERERETkcJrFERERE5HA4nICIiIiIHA57YomIiIjI4TCJJSIiIiKHwySWiIiIiBwOk1giIiIicjhMYomIiIjI4TCJJSIiIiKHwySWiIiIiBwOk1giIiIicjhMYomIiIjI4fw/t62edvdk9vAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(8,2)); subplot(121)\n",
    "plot(params.ravel(), params_est.ravel(), ',')\n",
    "grid(1), axis('equal'), xlabel('params'), ylabel('estimated params');\n",
    "subplot(122), hist(params.ravel() - params_est.ravel()), grid(1), xlabel('z scored true params');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0da43fb-e0c2-4adb-8315-0f89085f3b67",
   "metadata": {},
   "source": [
    "given how simple least squares is, this is already quite ok.\n",
    "\n",
    "## w/ Bayesian MLE instead\n",
    "\n",
    "The full SBI techniques use various (deep neural) predictors of parameters of the approximate posterior, but it doesn't work with a plain least squares, since some derivatives are required, so we can try with Jax,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "8943d961-adc2-4fe6-8817-ccdaf6ae8594",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax, jax.numpy as jp, jax.scipy as js"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "4be32af2-32b0-4911-b5eb-32e5ae3c2eac",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Array(-1.7370857, dtype=float32, weak_type=True)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "js.stats.norm.logpdf(0,1,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f4e7443-791a-4987-b5ff-3a72e28c25d6",
   "metadata": {},
   "source": [
    "instead of solving a least squares problem directly, we write a loss function with a linear Bayesian model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "c0cdd2e9-6b45-475f-8237-46f2351c2379",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(127184.44, dtype=float32)"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j_params = jp.array(params)      # num_samples, num_params\n",
    "j_features = jp.array(features)  # num_samples, num_features\n",
    "\n",
    "# predict a mean & std per parameter\n",
    "j_predictor = jp.ones((2, num_features, num_params))\n",
    "\n",
    "def loss(j_predictor):\n",
    "    mu = j_features @ j_predictor[0]\n",
    "    sd = j_features @ j_predictor[1]\n",
    "    lp_per = js.stats.norm.logpdf(j_params, mu, sd)\n",
    "    return -jp.sum(lp_per)\n",
    "\n",
    "loss(j_predictor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3670d68-d696-400e-8d8f-c331b071adcb",
   "metadata": {},
   "source": [
    "with a loss defined, we can descend the gradient,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "e8190a6e-acc7-4b15-99eb-ed9ad2731df7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "loss=270608.5625: 100%|█████████████████████████████████████████████████████████████████████████████| 10/10 [00:00<00:00, 81.80it/s]\n"
     ]
    }
   ],
   "source": [
    "jit_grad_loss = jax.jit(jax.grad(loss))\n",
    "\n",
    "num_iter = 10\n",
    "lr = 1e-3\n",
    "j_predictor = jp.ones((2, num_features, num_params))\n",
    "for i in (pbar := tqdm.trange(num_iter)):\n",
    "    j_predictor = j_predictor - lr * jit_grad_loss(j_predictor)\n",
    "    pbar.set_description(f'loss={loss(j_predictor)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9c38aba-af7d-4cb3-b5bf-0590032e570c",
   "metadata": {},
   "source": [
    "the loss stablizes quickly at ~270000 because there's no nonlinearity and the loss landscape is convex (i.e. easy).\n",
    "\n",
    "a simple way to evaluate the result is to z-score the true params based on the estimated mu & sd,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "dbdaa2d8-72a1-4044-810c-5824d6047362",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAADJCAYAAADmZGDwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAf8UlEQVR4nO3dfVSUdf7/8RcgDJICocLoaobZmqaFYOLofr1JBDfajfS0a5pRmaZhG7LrXadMrc2yG82yqK20PSfL3NU0rQxBaC0k8+aEWqwWm6aBJxXwFlA+vz/6MSdCBoeY4QKfj3Pm1FzX+5p5X28vr15dM1z4GGOMAAAAAIvybeoGAAAAAFcIrAAAALA0AisAAAAsjcAKAAAASyOwAgAAwNIIrAAAALA0AisAAAAsjcAKAAAASyOwAgAAwNIIrAAAALA0twProUOHdMcdd6hdu3Zq3bq1+vTpoy+++MK53hijOXPmqGPHjmrdurXi4uK0b9++Gq9x7NgxjRs3TsHBwQoNDdWECRN08uTJGjVffvml/u///k+BgYHq0qWLFi5c2MBdBAAAQHPWyp3i48ePa9CgQRo2bJg+/PBDdejQQfv27dPll1/urFm4cKGWLFmiN998U5GRkXrkkUeUkJCgvXv3KjAwUJI0btw4/fDDD8rIyFBlZaXuvvtuTZo0SStWrJAklZWVKT4+XnFxcUpPT1d+fr7uuecehYaGatKkSRfVa1VVlQ4fPqy2bdvKx8fHnd0EAACAFxhjdOLECXXq1Em+vi6uoxo3zJw50/zud7+rc31VVZWx2+3m6aefdi4rKSkxNpvNvP3228YYY/bu3WskmW3btjlrPvzwQ+Pj42MOHTpkjDHmpZdeMpdffrkpLy+v8d49evS46F4PHjxoJPHgwYMHDx48ePCw+OPgwYMuc51bV1jXrVunhIQE3XbbbcrJydFvfvMb3X///Zo4caIkqbCwUEVFRYqLi3NuExISotjYWOXm5mrMmDHKzc1VaGio+vXr56yJi4uTr6+v8vLydOuttyo3N1eDBw9WQECAsyYhIUFPPfWUjh8/XuOKbrXy8nKVl5c7nxtjnD21bdvWnd1skMrKSm3evFnDhg2Tv7+/x9+vuWE+9WNGrjEf15iPa8ynfszINebjWkPnc+LECUVGRtab1dwKrN9++61efvllpaWl6aGHHtK2bdv0l7/8RQEBAUpOTlZRUZEkKSIiosZ2ERERznVFRUUKDw+v2USrVgoLC6tRExkZWes1qtddKLAuWLBA8+bNq7U8NzdXQUFB7uxmgwUFBSkvL88r79UcMZ/6MSPXmI9rzMc15lM/ZuQa83GtIfM5ffq0JNX79U23AmtVVZX69eunJ554QpLUt29f7d69W+np6UpOTnarwcY2e/ZspaWlOZ+XlZWpS5cuio+PV3BwsMffv7KyUhkZGRoxYgT/53UBzKd+zMg15uMa83GN+dSPGbnGfFxr6HzKysouqs6twNqxY0f16tWrxrKePXvq3//+tyTJbrdLkoqLi9WxY0dnTXFxsaKiopw1R44cqfEa586d07Fjx5zb2+12FRcX16ipfl5d80s2m002m63Wcn9/f68eWN5+v+aG+dSPGbnGfFxjPq4xn/oxI9eYj2vuzudia926rdWgQYNUUFBQY9l///tfde3aVZIUGRkpu92uzMxM5/qysjLl5eXJ4XBIkhwOh0pKSrR9+3ZnTVZWlqqqqhQbG+us+eSTT1RZWemsycjIUI8ePS74dQAAAAC0XG4F1mnTpmnr1q164okntH//fq1YsUKvvvqqUlJSJP30/YPU1FQ9/vjjWrdunfLz83XnnXeqU6dOSkpKkvTTFdmRI0dq4sSJ+vzzz/Xpp59q6tSpGjNmjDp16iRJGjt2rAICAjRhwgTt2bNHK1eu1PPPP1/jI38AAABcGtz6SsANN9ygNWvWaPbs2Zo/f74iIyO1ePFijRs3zlkzY8YMnTp1SpMmTVJJSYl+97vf6aOPPnLeg1WS3nrrLU2dOlXDhw+Xr6+vRo8erSVLljjXh4SE6OOPP1ZKSopiYmLUvn17zZkz56LvwQoAAICWw63AKkk333yzbr755jrX+/j4aP78+Zo/f36dNWFhYc5fElCX6667Tv/5z3/cbQ8Amr0rZ21wexubn9HC/h5oBgAswO1fzQoAAAB4E4EVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpvyqwPvnkk/Lx8VFqaqpz2dmzZ5WSkqJ27dqpTZs2Gj16tIqLi2tsd+DAASUmJiooKEjh4eGaPn26zp07V6MmOztb0dHRstls6t69u5YvX/5rWgUAAEAz1eDAum3bNr3yyiu67rrraiyfNm2a3n//fa1atUo5OTk6fPiwRo0a5Vx//vx5JSYmqqKiQp999pnefPNNLV++XHPmzHHWFBYWKjExUcOGDdOuXbuUmpqqe++9Vxs3bmxouwAAAGimWjVko5MnT2rcuHH6xz/+occff9y5vLS0VK+//rpWrFihG2+8UZK0bNky9ezZU1u3btWAAQP08ccfa+/evdq0aZMiIiIUFRWlxx57TDNnztTcuXMVEBCg9PR0RUZG6tlnn5Uk9ezZU1u2bNGiRYuUkJBwwZ7Ky8tVXl7ufF5WViZJqqysVGVlZUN20y3V7+GN92qOmE/9mJFrl9J8bH7G/W18f9omZv5HKq/ycXv73XMvfG5tKS6l46ehmJFrzMe1hs7nYut9jDFunxmTk5MVFhamRYsWaejQoYqKitLixYuVlZWl4cOH6/jx4woNDXXWd+3aVampqZo2bZrmzJmjdevWadeuXc71hYWF6tatm3bs2KG+fftq8ODBio6O1uLFi501y5YtU2pqqkpLSy/Y09y5czVv3rxay1esWKGgoCB3dxEAAAAedvr0aY0dO1alpaUKDg6us87tK6zvvPOOduzYoW3bttVaV1RUpICAgBphVZIiIiJUVFTkrImIiKi1vnqdq5qysjKdOXNGrVu3rvXes2fPVlpamvN5WVmZunTpovj4eJcDaCyVlZXKyMjQiBEj5O/v7/H3a26YT/2YkWuX0nx6z3X/6082X6PH+lXpkS98ucJ6AZfS8dNQzMg15uNaQ+dT/Yl4fdwKrAcPHtSDDz6ojIwMBQYGurOpx9lsNtlstlrL/f39vXpgefv9mhvmUz9m5NqlMJ/y8+4HTue2VT4N2r6lz7TapXD8/FrMyDXm45q787nYWrd+6Gr79u06cuSIoqOj1apVK7Vq1Uo5OTlasmSJWrVqpYiICFVUVKikpKTGdsXFxbLb7ZIku91e664B1c/rqwkODr7g1VUAAAC0XG4F1uHDhys/P1+7du1yPvr166dx48Y5/93f31+ZmZnObQoKCnTgwAE5HA5JksPhUH5+vo4cOeKsycjIUHBwsHr16uWs+flrVNdUvwYAAAAuHW59JaBt27bq3bt3jWWXXXaZ2rVr51w+YcIEpaWlKSwsTMHBwXrggQfkcDg0YMAASVJ8fLx69eql8ePHa+HChSoqKtLDDz+slJQU50f6kydP1osvvqgZM2bonnvuUVZWlt59911t2LChMfYZAAAAzUiDbmvlyqJFi+Tr66vRo0ervLxcCQkJeumll5zr/fz8tH79ek2ZMkUOh0OXXXaZkpOTNX/+fGdNZGSkNmzYoGnTpun5559X586d9dprr9V5SysAAAC0XL86sGZnZ9d4HhgYqKVLl2rp0qV1btO1a1d98MEHLl936NCh2rlz569tDwAAAM3cr/rVrAAAAICnEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJZGYAUAAIClEVgBAABgaQRWAAAAWBqBFQAAAJbmVmBdsGCBbrjhBrVt21bh4eFKSkpSQUFBjZqzZ88qJSVF7dq1U5s2bTR69GgVFxfXqDlw4IASExMVFBSk8PBwTZ8+XefOnatRk52drejoaNlsNnXv3l3Lly9v2B4CAACgWXMrsObk5CglJUVbt25VRkaGKisrFR8fr1OnTjlrpk2bpvfff1+rVq1STk6ODh8+rFGjRjnXnz9/XomJiaqoqNBnn32mN998U8uXL9ecOXOcNYWFhUpMTNSwYcO0a9cupaam6t5779XGjRsbYZcBAADQnLRyp/ijjz6q8Xz58uUKDw/X9u3bNXjwYJWWlur111/XihUrdOONN0qSli1bpp49e2rr1q0aMGCAPv74Y+3du1ebNm1SRESEoqKi9Nhjj2nmzJmaO3euAgIClJ6ersjISD377LOSpJ49e2rLli1atGiREhISLthbeXm5ysvLnc/LysokSZWVlaqsrHRnNxuk+j288V7NEfOpHzNy7VKaj83PuL+Nr6nxT3e19LleSsdPQzEj15iPaw2dz8XW+xhjGnZ2k7R//35dffXVys/PV+/evZWVlaXhw4fr+PHjCg0NddZ17dpVqampmjZtmubMmaN169Zp165dzvWFhYXq1q2bduzYob59+2rw4MGKjo7W4sWLnTXLli1TamqqSktLL9jL3LlzNW/evFrLV6xYoaCgoIbuIgAAADzk9OnTGjt2rEpLSxUcHFxnnVtXWH+uqqpKqampGjRokHr37i1JKioqUkBAQI2wKkkREREqKipy1kRERNRaX73OVU1ZWZnOnDmj1q1b1+pn9uzZSktLcz4vKytTly5dFB8f73IAjaWyslIZGRkaMWKE/P39Pf5+zQ3zqR8zcu1Smk/vue5//cnma/RYvyo98oWvyqt83N5+99wLf3rVUlxKx09DMSPXmI9rDZ1P9Sfi9WlwYE1JSdHu3bu1ZcuWhr5Eo7LZbLLZbLWW+/v7e/XA8vb7NTfMp37MyLVLYT7l590PnM5tq3watH1Ln2m1S+H4+bWYkWvMxzV353OxtQ26rdXUqVO1fv16bd68WZ07d3Yut9vtqqioUElJSY364uJi2e12Z80v7xpQ/by+muDg4AteXQUAAEDL5VZgNcZo6tSpWrNmjbKyshQZGVljfUxMjPz9/ZWZmelcVlBQoAMHDsjhcEiSHA6H8vPzdeTIEWdNRkaGgoOD1atXL2fNz1+juqb6NQAAAHDpcOsrASkpKVqxYoXWrl2rtm3bOr9zGhISotatWyskJEQTJkxQWlqawsLCFBwcrAceeEAOh0MDBgyQJMXHx6tXr14aP368Fi5cqKKiIj388MNKSUlxfqQ/efJkvfjii5oxY4buueceZWVl6d1339WGDRsaefcBAABgdW5dYX355ZdVWlqqoUOHqmPHjs7HypUrnTWLFi3SzTffrNGjR2vw4MGy2+1avXq1c72fn5/Wr18vPz8/ORwO3XHHHbrzzjs1f/58Z01kZKQ2bNigjIwMXX/99Xr22Wf12muv1XlLKwAAALRcbl1hvZg7YAUGBmrp0qVaunRpnTVdu3bVBx984PJ1hg4dqp07d7rTHgAAAFqgBv3QFQAAAOAtBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKURWAEAAGBpBFYAAABYGoEVAAAAlkZgBQAAgKVZOrAuXbpUV155pQIDAxUbG6vPP/+8qVsCAACAl7Vq6gbqsnLlSqWlpSk9PV2xsbFavHixEhISVFBQoPDw8KZuDwBcunLWhqZuAQBaDMsG1ueee04TJ07U3XffLUlKT0/Xhg0b9MYbb2jWrFm16svLy1VeXu58XlpaKkk6duyYKisrPd5vZWWlTp8+raNHj8rf39/j79fcMJ/6MSPXmtt8Wp075d33qzI6fbpKrSp9db7Kx+3tjx496oGurKO5HT9NgRm5xnxca+h8Tpw4IUkyxris8zH1VTSBiooKBQUF6V//+peSkpKcy5OTk1VSUqK1a9fW2mbu3LmaN2+eF7sEAABAYzh48KA6d+5c53pLXmH98ccfdf78eUVERNRYHhERoa+//vqC28yePVtpaWnO51VVVTp27JjatWsnHx/3rza4q6ysTF26dNHBgwcVHBzs8fdrbphP/ZiRa8zHNebjGvOpHzNyjfm41tD5GGN04sQJderUyWWdJQNrQ9hsNtlsthrLQkNDvd5HcHAwB7ILzKd+zMg15uMa83GN+dSPGbnGfFxryHxCQkLqrbHkXQLat28vPz8/FRcX11heXFwsu93eRF0BAACgKVgysAYEBCgmJkaZmZnOZVVVVcrMzJTD4WjCzgAAAOBtlv1KQFpampKTk9WvXz/1799fixcv1qlTp5x3DbAam82mRx99tNbXEvAT5lM/ZuQa83GN+bjGfOrHjFxjPq55ej6WvEtAtRdffFFPP/20ioqKFBUVpSVLlig2Nrap2wIAAIAXWTqwAgAAAJb8DisAAABQjcAKAAAASyOwAgAAwNIIrAAAALA0AmsD/e9//9OECRMUGRmp1q1b66qrrtKjjz6qiooKl9udPXtWKSkpateundq0aaPRo0fX+gUJLcXf//53DRw4UEFBQRf9W8fuuusu+fj41HiMHDnSs402kYbMxxijOXPmqGPHjmrdurXi4uK0b98+zzbahI4dO6Zx48YpODhYoaGhmjBhgk6ePOlym6FDh9Y6hiZPnuyljj1r6dKluvLKKxUYGKjY2Fh9/vnnLutXrVqla665RoGBgerTp48++OADL3XaNNyZz/Lly2sdJ4GBgV7s1rs++eQT/eEPf1CnTp3k4+Oj9957r95tsrOzFR0dLZvNpu7du2v58uUe77OpuDuf7OzsWsePj4+PioqKvNOwly1YsEA33HCD2rZtq/DwcCUlJamgoKDe7RrzHERgbaCvv/5aVVVVeuWVV7Rnzx4tWrRI6enpeuihh1xuN23aNL3//vtatWqVcnJydPjwYY0aNcpLXXtXRUWFbrvtNk2ZMsWt7UaOHKkffvjB+Xj77bc91GHTash8Fi5cqCVLlig9PV15eXm67LLLlJCQoLNnz3qw06Yzbtw47dmzRxkZGVq/fr0++eQTTZo0qd7tJk6cWOMYWrhwoRe69ayVK1cqLS1Njz76qHbs2KHrr79eCQkJOnLkyAXrP/vsM91+++2aMGGCdu7cqaSkJCUlJWn37t1e7tw73J2P9NOvkPz5cfLdd995sWPvOnXqlK6//notXbr0ouoLCwuVmJioYcOGadeuXUpNTdW9996rjRs3erjTpuHufKoVFBTUOIbCw8M91GHTysnJUUpKirZu3aqMjAxVVlYqPj5ep06dqnObRj8HGTSahQsXmsjIyDrXl5SUGH9/f7Nq1Srnsq+++spIMrm5ud5osUksW7bMhISEXFRtcnKyueWWWzzaj9Vc7HyqqqqM3W43Tz/9tHNZSUmJsdls5u233/Zgh01j7969RpLZtm2bc9mHH35ofHx8zKFDh+rcbsiQIebBBx/0Qofe1b9/f5OSkuJ8fv78edOpUyezYMGCC9b/6U9/MomJiTWWxcbGmvvuu8+jfTYVd+fjznmppZFk1qxZ47JmxowZ5tprr62x7M9//rNJSEjwYGfWcDHz2bx5s5Fkjh8/7pWerObIkSNGksnJyamzprHPQVxhbUSlpaUKCwurc/327dtVWVmpuLg457JrrrlGV1xxhXJzc73RYrOQnZ2t8PBw9ejRQ1OmTNHRo0ebuiVLKCwsVFFRUY3jJyQkRLGxsS3y+MnNzVVoaKj69evnXBYXFydfX1/l5eW53Patt95S+/bt1bt3b82ePVunT5/2dLseVVFRoe3bt9f4s/f19VVcXFydf/a5ubk16iUpISGhRR4rDZmPJJ08eVJdu3ZVly5ddMstt2jPnj3eaLdZuJSOn18jKipKHTt21IgRI/Tpp582dTteU1paKkkuM09jH0OW/dWszc3+/fv1wgsv6JlnnqmzpqioSAEBAbW+rxgREdFiv/firpEjR2rUqFGKjIzUN998o4ceeki///3vlZubKz8/v6Zur0lVHyMRERE1lrfU46eoqKjWx2utWrVSWFiYy/0dO3asunbtqk6dOunLL7/UzJkzVVBQoNWrV3u6ZY/58ccfdf78+Qv+2X/99dcX3KaoqOiSOVYaMp8ePXrojTfe0HXXXafS0lI988wzGjhwoPbs2aPOnTt7o21Lq+v4KSsr05kzZ9S6desm6swaOnbsqPT0dPXr10/l5eV67bXXNHToUOXl5Sk6Orqp2/OoqqoqpaamatCgQerdu3eddY19DuIK6y/MmjXrgl+k/vnjlyfAQ4cOaeTIkbrttts0ceLEJurcOxoyH3eMGTNGf/zjH9WnTx8lJSVp/fr12rZtm7KzsxtvJzzI0/NpCTw9o0mTJikhIUF9+vTRuHHj9M9//lNr1qzRN99804h7gebO4XDozjvvVFRUlIYMGaLVq1erQ4cOeuWVV5q6NTQDPXr00H333aeYmBgNHDhQb7zxhgYOHKhFixY1dWsel5KSot27d+udd97x6vtyhfUX/vrXv+quu+5yWdOtWzfnvx8+fFjDhg3TwIED9eqrr7rczm63q6KiQiUlJTWushYXF8tut/+atr3G3fn8Wt26dVP79u21f/9+DR8+vNFe11M8OZ/qY6S4uFgdO3Z0Li8uLlZUVFSDXrMpXOyM7HZ7rR+YOXfunI4dO+bW35fY2FhJP30KctVVV7ndrxW0b99efn5+te4o4urcYbfb3apvzhoyn1/y9/dX3759tX//fk+02OzUdfwEBwdf8ldX69K/f39t2bKlqdvwqKlTpzp/ALa+TyIa+xxEYP2FDh06qEOHDhdVe+jQIQ0bNkwxMTFatmyZfH1dX7COiYmRv7+/MjMzNXr0aEk//YThgQMH5HA4fnXv3uDOfBrD999/r6NHj9YIaFbmyflERkbKbrcrMzPTGVDLysqUl5fn9p0YmtLFzsjhcKikpETbt29XTEyMJCkrK0tVVVXOEHoxdu3aJUnN5hi6kICAAMXExCgzM1NJSUmSfvpYLjMzU1OnTr3gNg6HQ5mZmUpNTXUuy8jIaDbnGnc0ZD6/dP78eeXn5+umm27yYKfNh8PhqHULopZ6/DSWXbt2NevzjCvGGD3wwANas2aNsrOzFRkZWe82jX4OatCPasF8//33pnv37mb48OHm+++/Nz/88IPz8fOaHj16mLy8POeyyZMnmyuuuMJkZWWZL774wjgcDuNwOJpiFzzuu+++Mzt37jTz5s0zbdq0MTt37jQ7d+40J06ccNb06NHDrF692hhjzIkTJ8zf/vY3k5ubawoLC82mTZtMdHS0ufrqq83Zs2ebajc8xt35GGPMk08+aUJDQ83atWvNl19+aW655RYTGRlpzpw50xS74HEjR440ffv2NXl5eWbLli3m6quvNrfffrtz/S//ju3fv9/Mnz/ffPHFF6awsNCsXbvWdOvWzQwePLipdqHRvPPOO8Zms5nly5ebvXv3mkmTJpnQ0FBTVFRkjDFm/PjxZtasWc76Tz/91LRq1co888wz5quvvjKPPvqo8ff3N/n5+U21Cx7l7nzmzZtnNm7caL755huzfft2M2bMGBMYGGj27NnTVLvgUSdOnHCeYySZ5557zuzcudN89913xhhjZs2aZcaPH++s//bbb01QUJCZPn26+eqrr8zSpUuNn5+f+eijj5pqFzzK3fksWrTIvPfee2bfvn0mPz/fPPjgg8bX19ds2rSpqXbBo6ZMmWJCQkJMdnZ2jbxz+vRpZ42nz0EE1gZatmyZkXTBR7XCwkIjyWzevNm57MyZM+b+++83l19+uQkKCjK33nprjZDbkiQnJ19wPj+fhySzbNkyY4wxp0+fNvHx8aZDhw7G39/fdO3a1UycONH5H5yWxt35GPPTra0eeeQRExERYWw2mxk+fLgpKCjwfvNecvToUXP77bebNm3amODgYHP33XfXCPS//Dt24MABM3jwYBMWFmZsNpvp3r27mT59uiktLW2iPWhcL7zwgrniiitMQECA6d+/v9m6datz3ZAhQ0xycnKN+nfffdf89re/NQEBAebaa681GzZs8HLH3uXOfFJTU521ERER5qabbjI7duxogq69o/o2TL98VM8kOTnZDBkypNY2UVFRJiAgwHTr1q3GuailcXc+Tz31lLnqqqtMYGCgCQsLM0OHDjVZWVlN07wX1JV3fn5MePoc5PP/GwEAAAAsibsEAAAAwNIIrAAAALA0AisAAAAsjcAKAAAASyOwAgAAwNIIrAAAALA0AisAAAAsjcAKAAAASyOwAgAAwNIIrAAAALA0AisAAAAs7f8BZS+vhdk7NRwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu = j_features @ j_predictor[0]\n",
    "sd = j_features @ j_predictor[1]\n",
    "z = (j_params - mu) / sd\n",
    "figure(figsize=(8,2)); hist(z.ravel(), np.r_[-2:2:0.1]), grid(1);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76a2b662-9166-416a-92fd-2491ce970bb6",
   "metadata": {},
   "source": [
    "with z scores close to zero, then we can be confident that this has worked."
   ]
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}