{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"navigation": true
},
"source": [
"\n",
"< []() | [Chapter 1. The Golem of Prague](01-the-golem-of-prague.html) >"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"!pip install -q numpyro arviz"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"\n",
"import arviz as az\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"\n",
"import jax.numpy as jnp\n",
"from jax import random\n",
"\n",
"import numpyro\n",
"import numpyro.distributions as dist\n",
"import numpyro.optim as optim\n",
"from numpyro.infer import SVI, Trace_ELBO\n",
"\n",
"if \"SVG\" in os.environ:\n",
" %config InlineBackend.figure_formats = [\"svg\"]\n",
"az.style.use(\"arviz-darkgrid\")\n",
"numpyro.set_platform(\"cpu\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 0.1"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2017-12-07T01:31:33.358577Z",
"start_time": "2017-12-07T01:31:33.356029Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"All models are wrong, but some are useful.\n"
]
}
],
"source": [
"print(\"All models are wrong, but some are useful.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 0.2"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2017-12-07T01:31:33.433222Z",
"start_time": "2017-12-07T01:31:33.359596Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"DeviceArray(200.00002, dtype=float32)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = jnp.arange(1, 3)\n",
"x = x * 10\n",
"x = jnp.log(x)\n",
"x = jnp.sum(x)\n",
"x = jnp.exp(x)\n",
"x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 0.3"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-inf\n",
"-921.03406\n"
]
}
],
"source": [
"print(jnp.log(0.01**200))\n",
"print(200 * jnp.log(0.01))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 0.4"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2017-12-07T01:31:33.798552Z",
"start_time": "2017-12-07T01:31:33.508571Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 1000/1000 [00:00<00:00, 1165.66it/s, init loss: 30629.4453, avg. loss [951-1000]: 5608.2529]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"{'a': DeviceArray(-17.579102, dtype=float32), 'b': DeviceArray(3.9324093, dtype=float32)}\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAO0lEQVR4nO3dfXyT1cH/8W8KFIkttUIjUpQnLVYUJrYUgcnQDrxt9QZHgdfgBsQqPwUcDEGdOh83HN5sCnQTZYIgOCjCdKIvkZWJKFCgKKgtT9oqBRvEUlqqUMn1+6N3IrUppCRXk1z5vP/RnpwkJ7kOV78czoPNMAxDAAAAAAIuKtgNAAAAAKyKsA0AAACYhLANAAAAmISwDQAAAJiEsA0AAACYhLANAAAAmISwDQAAAJiEsA0AAACYpHmwG2AF5eXlwW6C6eLi4lRRURHsZiAM0FfgC/oJfEVfga+C0Vfi4+PPWoeRbfgkKoquAt/QV+AL+gl8RV+Br0K1r4RmqwAAAAALIGwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAACCgnE5DBTsMOZ1GsJsSdM2D3QAAAABYx5trDM2abcjlkqKipBnTpMwMW7CbFTSMbAMAACAgnM4fg7YkuVzSrNmRPcJN2AYAAEBAHCiVJ2i7uVy15ZGKsA0AAICA6JBYO3XkdFFRteWRirANAGAxE4CAcDhsmjHN5gnctXO2bXI4InfONgskASDCsZgJQCBlZtjUO7V26kiHREV00JYY2QaAiMZiJgBmcDhs6nVNZI9ouxG2ASCCsZgJAMxF2AaACMZiJgAwF2EbACIYi5kAwFwskASACMdiJgAwD2EbACCHwyaHI9itAADriZhpJC+++KK6deumbt266aOPPvJap6qqSjNnztTAgQN11VVXaeDAgZo5c6aqqqqatrEAAACwhIgI2/v379ecOXNkt9sbrFNdXa3Ro0dr0aJF6ty5s8aNG6euXbtq0aJFGj16tKqrq5uwxQAAALACy4ftU6dO6f7779cVV1yh9PT0BustWLBAhYWFys7O1ksvvaT77rtPCxYs0MSJE1VYWKgFCxY0YasBAABgBZYP2y+++KKKior0xz/+Uc2aNfNaxzAM5ebmym63a+LEiXUemzBhguLi4rRy5UoZBoc8AAAAwHeWDtt79uzRvHnzdPfdd+vyyy9vsF5xcbGcTqd69epVb6pJy5YtlZKSorKyMpWUlJjdZAAAAFiIZcP2Dz/8oAceeEBdu3bVXXfddca67hDdqVMnr4937NixTj0AAADAF5bd+u/555/X7t27tWLFCrVo0eKMdSsrKyVJMTExXh93l7vr/VRcXJyifnoEmwXFx8cHuwkIE/QV+IJ+Al/RV+CrUOwrlgzbRUVFev755zV+/Hh1797d9PerqKgw/T2CLT4+XuXl5cFuBsIAfQW+oJ/AV/QV+CoYfcWXcG/J4dj7779fl1xyiSZPnuxT/djYWElqcD9td7m7HgAAAOALy45sS9LVV1/t9fERI0ZIknJycpSenu6Zk11cXOy1vnuutrseAAAA4AtLhu1hw4Z5Ld+2bZuKi4t1ww036MILL1RiYqKk2oWRDodDBQUFqq6urrMjyYkTJ7Rt2zY5HA7CNgAAABrFkmH7D3/4g9fyBx54QMXFxZowYYJ+9rOfecptNpuysrKUk5OjnJwcTZ8+3fPY/PnzVVFRoYkTJ8pms5nddAAAAFiIJcP2ucjOzlZeXp7nJMnu3burqKhIGzZsUHJysrKzs4PdRAAAAIQZSy6QPBd2u11LlizRuHHj9Pnnn2vhwoXau3evxo0bpyVLltQ77AYAAAA4G5vBGeR+i4Qtidh6Cb6ir8AX9BP4ir4CX7H1HwAAABBhCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASQjbAAAAgEkI2wAAAIBJCNsAAACASZoHuwFmOHbsmObMmaNdu3bpwIEDqqioUHx8vDp37qxRo0Zp0KBBstlsdZ5TVVWluXPnau3atTp8+LASEhI0aNAgTZ48WTExMUH6JAAAAAhnlhzZLi8v12uvvaZWrVrpxhtv1Pjx43X99ddr3759uvfee/X73/++Tv3q6mqNHj1aixYtUufOnTVu3Dh17dpVixYt0ujRo1VdXR2kTwIAAIBwZsmR7Q4dOmjr1q1q3rzux6uqqtKIESO0YsUKjRkzRpdffrkkacGCBSosLFR2dramT5/uqT9nzhzl5ORowYIFuvfee5v0MwAAACD8WXJku1mzZvWCtiTFxMSof//+kqSSkhJJkmEYys3Nld1u18SJE+vUnzBhguLi4rRy5UoZhmF+wwEAaCSn01DBDkNOJ7+ngFBkybDdkBMnTmjz5s2y2Wy67LLLJEnFxcVyOp3q1auX7HZ7nfotW7ZUSkqKysrKPOEcAIBQ8eYaQ8NGGrp3au1/31xD4AZCjSWnkbgdO3ZML7/8slwul44cOaINGzbo0KFDmjRpkjp16iTpxxFu988/1bFjR0+9huoAANDUnE5Ds2Ybcrlqf3a5pFmzDfVOlRwO25mfDKDJWD5sz5s3z/NzixYtNGPGDI0fP95TVllZKUkN7jjiLnfX8yYuLk5RUdb/R4L4+PhgNwFhgr4CX9BP/LNnb41crmN1ylwuqeJYrLp1axGkVpmDvgJfhWJfsXTY7tChg3bv3q1Tp07p0KFDeuutt/SXv/xFO3bs0LPPPut1Xve5qKioCMjrhLL4+HiVl5cHuxkIA/QV+IJ+4r8L4gxFRckzsi1JUVFSXOtKlZdbZ2SbvgJfBaOv+BLurT8cq9oFkx06dNBdd92lKVOm6N1339WKFSskSbGxsZJqdyrxxl3urgcAQChwOGyaMc0m9z+sRkVJM6bZmEIChJiICNunc+9Gkp+fL+nHOdnFxcVe67vndLvrAQAQKjIzbFr5D5vm/KX2v5kZBG0g1Fh6Gok3ZWVlkmpHu6XahZEOh0MFBQWqrq6usyPJiRMntG3bNjkcDsI2ACAkORw2ORzBbgWAhlhyZLuwsNDrgsajR4/qL3/5iyTp+uuvlyTZbDZlZWWpurpaOTk5derPnz9fFRUVysrKqne8OwAAAHA2lhzZXrVqlVauXKm0tDS1b99erVq10sGDB/Wf//xH1dXVGjx4sG655RZP/ezsbOXl5XlOkuzevbuKioq0YcMGJScnKzs7O4ifBgAAAOHKkmF78ODBqqqq0kcffaStW7fq+++/V1xcnK699loNGTJEGRkZdUaq7Xa7lixZonnz5umdd95Rfn6+2rZtq3HjxmnSpEn1DrsBAAAAfGEzOIfcb5GwJRFbL8FX9BX4gn4CX9FX4Cu2/gMAAAAiDGEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAIKfTUMEOQ06nEeymoJG4dkBoax7sBgAAguvNNYZmzTbkcklRUdKMaVJmhi3YzYIPuHZA6GNkGwAimNP5Y1iTJJdLmjWbUdJwwLUDwgNhGwAi2IFSecKam8tVW47QxrUDwgNhGwAiWIfE2ukHp4uKqi1HaOPaAeGBsA0AEczhsGnGNJsntNXO+7XJ4WDeb6jj2gHhgQWSABDhMjNs6p1aO/2gQ6IIa2GEaweEPsI2AEAOh00OR7BbgXPBtQNCmyXDdllZmd5++21t2LBBn3/+ub755hvFxcWpV69eys7OVs+ePes9p6qqSnPnztXatWt1+PBhJSQkaNCgQZo8ebJiYmKC8CkAAAAQ7iwZtpcsWaIXX3xRl156qfr27as2bdqopKRE69at07p16zR79mzdfPPNnvrV1dUaPXq0CgsL1a9fP2VkZKioqEiLFi3Sli1btGzZMtnt9iB+IgAAAIQjS4btHj16aOnSpUpJSalTvm3bNo0bN06PP/640tPTFR0dLUlasGCBCgsLlZ2drenTp3vqz5kzRzk5OVqwYIHuvffeJv0MAAAACH+W3I1k0KBB9YK2JKWkpCgtLU1Hjx7V7t27JUmGYSg3N1d2u10TJ06sU3/ChAmKi4vTypUrZRgcEgAAAIDGsWTYPpPmzZvX+W9xcbGcTqd69epVb6pIy5YtlZKSorKyMpWUlDR5WwEAABDeIipsHzx4UB9++KESEhKUlJQkSZ4Q3alTJ6/P6dixY516AAAAgK8sOWfbm5qaGs2YMUMnT57Ufffdp2bNmkmSKisrJanBHUfc5e563sTFxSnqp8d4WVB8fHywm4AwQV+BL+gn8BV9Bb4Kxb4SEWHb5XLpd7/7nbZu3arhw4dryJAhAX39ioqKgL5eKIqPj1d5eXmwm4EwQF+BL+gn8BV9Bb4KRl/xJdxbfjjWMAw9/PDDeuONN3Trrbfq8ccfr/N4bGyspNp9tr1xl7vrAQAAAL6ydNh2j2i/9tpryszM1NNPP11vuod7TnZxcbHX13DP1XbXAwAAAHxl2bDtcrn00EMPadWqVbr55ps1a9Yszzzt03Xq1EkOh0MFBQWqrq6u89iJEye0bds2ORwOwjYAAAAazZJh+/SgfdNNN+mZZ57xGrQlyWazKSsrS9XV1crJyanz2Pz581VRUaGsrCzZbLamaDoAAAAsxJILJHNycrRq1SrZ7XZ16tRJf/vb3+rVSU9PV3JysiQpOztbeXl5npMku3fvrqKiIm3YsEHJycnKzs5u6o8AAAAAC7Bk2C4tLZUkVVdX6/nnn/daJzEx0RO27Xa7lixZonnz5umdd95Rfn6+2rZtq3HjxmnSpEn1DrsBAAAAfGEzOIfcb5GwJRFbL8FX9BX4gn4CX9FX4Cu2/gMAAAAiDGEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAB5XQaKthhyOnkoPLmwW4AAAAArOPNNYZmzTbkcklRUdKMaVJmhi3YzQoaRrYBAAAQEE7nj0FbklwuadbsyB7hJmwDAAAgIA6UyhO03Vyu2vJIRdgGAABAQHRIrJ06crqoqNrySEXYBgAAQEA4HDbNmGbzBO7aOds2ORyRO2ebBZIAAAAImMwMm7p0MbTrE+nqq6QrkyM3aEuEbQAAAARQ7W4kOm03EoPdSAAAAAB/sRtJfYRtAAAABAS7kdRH2AYAAEBAsBtJfYRtAAAABAS7kdTHAkkAAAAETGaGTb1Ta6eOdEhURAdtibANAACAAHM4bHI4gt2K0MA0EgAAAMAkhG0AAADAJIRtAAAAwCSEbQAAAMAkhG0AAADAJIRtAAAAwCSEbQAAAMAkhG0AAADAJIRtRDyn01DBDkNOpxHspoQ9p9PQlvwavssAsHK/pJ8AiCScIImI9uYaQ7NmG3K5pKgoaca02mNm0Xg/fpfH+C79ZOV+ST8BEGlshmEwtOCn8vLyYDfBdPHx8Zb7nE6noWEjawONW1SUtPIfNjkc/PJvDL7LwLHyd2nlzwbzWPH3D8wRjL4SHx9/1jpMI0HEOlCqOr/0pdqfD5QGpz3hjO8ycKz8XVr5swFAQwjbiFgdEmtH1U4XFVVbjsbhuwwcK3+XVv5sANAQwnaYsvLiqabicNg0Y5rN88u/dv4o/5x9LvguA8fK36WVPxsANIQ52wHQ1POD6i+espm+wMjKc+acTkMHSmtH1/il7x+n01DFsVjFta7ku/STlfsl/QSNYeXfPwisUJ2zTdgOgKa8sMFaYMTNDr6ir8AX9BP4ir4CX4Vq2Lbs1n+vv/66tm/frk8++UR79uxRTU2NZs6cqdtuu81r/aqqKs2dO1dr167V4cOHlZCQoEGDBmny5MmKiYlp4tY37EwLjByO4LQJAAAA3jUqbN94443n/EY2m03r1q075+c31nPPPafS0lLFx8fL4XCotLTh5e7V1dUaPXq0CgsL1a9fP2VkZKioqEiLFi3Sli1btGzZMtnt9iZr+5m4Fxj9dGSbBUYAAAChp1Fh+0yBtSE2m03BmKny1FNPqWPHjkpMTNQLL7yg2bNnN1h3wYIFKiwsVHZ2tqZPn+4pnzNnjnJycrRgwQLde++9TdHss6pdYKR6c7aZ9wgAMJuV1xIAZmlU2C4qKqpX9sQTT+iNN97QqFGjlJGRofbt20uSDh48qDVr1mjZsmW65ZZb9MgjjwSmxT7q27evT/UMw1Bubq7sdrsmTpxY57EJEybolVde0cqVKzV58mTZbKFxY8nMsKl3qrjhAQCajJVPNgXM5Nec7UWLFik3N1crVqxQcnJynceSkpKUlJSkm266ScOHD9cll1yi22+/3a/GmqG4uFhOp1P9+/evN1WkZcuWSklJ0b///W+VlJSoU6dOwWmkFw6HjTnaAIAm4XT+GLSl2qmMs2Yb6p3KgA9wNn6F7eXLl6tPnz71gvbpkpOT1adPH61YsSIkw3ZJSYkkNRikO3bs6KnXUJ24uDhF/fSkBgvyZcUtINFX4Bv6SfjYs7dGLtexOmUul1RxLFbdurUw/f3pK/BVKPYVv8J2aWmpunXrdtZ6559//jnN924KlZWVktTgjiPucnc9byoqKgLfsBDD1kvwFX0FvqCfhJcL4gyvi/PjWleqvJxzHhAaQnXrP7+GY9u0aaP8/HwdP368wTpVVVXasmWLLrzwQn/eCgAABAmnfwLnzq+wffPNN+vbb7/VHXfcoY8//rje4x9//LHuvPNOHT16VBkZGf68lWliY2Ml1f6lwBt3ubseAACRKDPDppX/sGnOX2r/y+JIwDd+TSOZNGmStm/fro8++kgjR47URRddpIsvvlg2m00HDx5UWVmZDMNQz549NWnSpEC1OaDcc7KLi4u9Pu6e0+2uBwBApGJxPtB4foXtVq1aafHixVq4cKFeffVVff311/r66689j7dr104jRozQHXfcoejoaL8ba4ZOnTrJ4XCooKBA1dXVdXYkOXHihLZt2yaHw0HYBgAAQKP5fVx7dHS0JkyYoAkTJujQoUNyOp0yDEMOh8Oz53Yos9lsysrKUk5OjnJycuocajN//nxVVFRo4sSJIbPHNgAAAMKH32H7dBdffLEuvvjiQL7kOcvNzdX27dslSXv27PGU5efnS5LS09OVnp4uScrOzlZeXp7nJMnu3burqKhIGzZsUHJysrKzs4PzIQAAABDWAhq2Q8n27du1evXqOmUFBQUqKCiQJCUmJnrCtt1u15IlSzRv3jy98847ys/PV9u2bTVu3DhNmjSp3mE3AAAAgC9shmEYvlZ+8MEHZbPZ9Nvf/lZt27bVgw8+6Psb2Wz64x//eE6NDHWRsP8n+5zCV/QV+IJ+Al/RV+CrUN1nu1Ej26tXr5bNZtOdd96ptm3b1hs5PhMrh20AAADAm0aF7cWLF0uSZ+Gj+2cAAAAA9TUqbPfu3fuMPwMAAAD4kV8nSAIAAABomF+7kXz33XcqLy/XBRdcUGfHjsrKSr3wwgvas2eP2rdvr/Hjx+uSSy7xu7EAAABAOPErbP/tb3/Tiy++qBUrVujqq6+WJJ08eVIjRozQF198IfdGJ2vXrtXrr7+utm3b+t9iAAAAIEz4NY1k06ZN6tChgydoS9Kbb76pzz//XGlpafr73/+usWPH6siRI1q0aJG/bQUAAADCil9h+9ChQ+rUqVOdsnXr1ikqKkpPP/20+vXrpwcffFCdO3fWe++9589bAQAAAGHHr7BdUVGh1q1b1ynbsWOHkpKS1K5dO09Zt27ddOjQIX/eCgAAAAg7foXthIQEOZ1Oz8979+5VeXm5UlNT69Sz2Wz+vA0AAAAQlvwK28nJydqxY4cKCwslSYsWLZLNZtPAgQPr1CspKZHD4fDnrQAAAICw49duJHfddZfWr1+vX/3qV4qNjVVFRYWSk5PVp08fT50jR46oqKhIGRkZfjcWAAAACCd+jWz37NlTf/3rX3Xttdeqbdu2uvXWW/W3v/1NUVE/vuy//vUvnX/++fr5z3/ud2MBAACAcGIz3Jth45yVl5cHuwmmi4+Pj4jPCf/RV+AL+gl8RV+Br4LRV+Lj489ah+PaAQAAAJP4NWfb7dtvv9Ubb7yhXbt26ejRo+rTp4/uvPNOSdKePXv01VdfqW/fvmrVqlUg3g4AAAAIC36H7TVr1uiRRx7Rd999J8MwZLPZ6uw8UlJSonvvvVczZ87UkCFD/H07AAAAIGz4NY1k27Ztmj59uqKjo/Xggw9q5cqV+ukU8F/84heKjY3Vu+++61dDAQAAgHDj18j2/Pnz1bx5cy1atEhXXHGF1zotWrRQly5dtG/fPn/eCgAAAAg7fo1sf/zxx+rZs2eDQdutXbt2dU6aBGBNTqehLfk1cjrZ5AgA0HRC+fePXyPb33//vU9bnlRVVXFkO2Bxb64xNGu2IZfrmKKipBnTpMwM/twDAMwV6r9//BrZbt++vXbv3n3GOj/88IN2796tjh07+vNWAEKY0+m+0dX+7HJJs2YbITnCAACwjnD4/eNX2B44cKC+/PJLLV26tME6Cxcu1DfffKP09HR/3gpACDtQKs+Nzs3lqi0HAMAs4fD7x69pJHfeeafWrFmjp556Sh999JFuvPFGSbX7bq9fv17r1q3T6tWrdfHFF2vMmDEBaTCA0NMhUYqKqnvDi4qqLQcAwCzh8PvH7+PaP//8c/3mN7/R3r17ZbPZPHttS5JhGOrSpYvmzZunLl26BKTBoSgSjpHluFyczY9z5vR/c+ZsITVnDqGFewp8RV/B2QTz948vaxf9GtmeOXOmWrdurddff115eXn68MMPVVpaqlOnTqldu3bq27evBg8erGbNmvnzNvDC6TR0oLT2b24OB4EGwZeZYVPvVKniWKziWlfSLwEATSLUf//4FbaXLl2qG264QVFRUUpPT2dedhOp/ze40Fp1i8jlcNjUrVsLlZfTHwEATSeUf//4tUDyoosuqndiJMwVDqtuAQAAUMuvsP3LX/5S+fn5qqqqClR7cBbhsOoWAAAAtfwK25MmTVL79u1111136bPPPgtUm3AG7lW3pwu1VbcAAACo5dec7XvuuUfR0dEqKCjQr371KyUkJKh9+/aKjo6uV9dms+nll1/25+2g2jlJM6ap3qrbUFsMAO9Y2BpYfJ+Atf8cOJ2G9uyt0QVxhuU+GyKHX2E7Pz/f8/+GYcjpdMrpdHqty3HtgeNedWvVm6tVsbA1sPg+AWv/OQj1I7gBX/m1z3ZpaeMmCicmWnOuQyTs/8k+p/5xOg0NG2nU23R/5T+s968STdFXIun7tCruKf6z8p8DK382mCcY9xXT99m2angGAu1MC1sdjuC0KZzxfQLW/nNg5c+GyOPXAkkAvmFha2DxfQLW/nNg5c+GyEPYBppA7cJWm+eXBwtb/cP3CVj7z4GVPxsij19ztlErEuYdMr8yMKy8a4BbU/aVSPg+rYp7SuBY+c+B02mE7BHcCD2WnLMNoHEcDhvzDQOI7xOw9p+DUD6CG/AV00hOs3PnTt15551KTU3Vz372Mw0bNkz/+te/gt0sAAAAhClGtv/Pli1bdMcdd6hFixbKyMhQbGys1q5dq/vuu0+lpaX6f//v/wW7iQAAAAgzhG1JP/zwgx5++GHZbDYtXbpUV155pSRp4sSJGjlypObOnaubbrpJnTp1Cm5DAQAAEFaYRiJp8+bN+vLLL5WZmekJ2pIUExOje+65Rz/88INWrVoVxBYCiDSfFbr0jxUufVboOntlhBSn01DBDkNOJ/sPAGBkW9KPx87379+/3mP9+vWrUwcAzPaHmS69/c6PP//XYJceepCxkXBg5ePTAZwb7t6SiouLJUkdO3as91hcXJzi4+NVUlLSxK0CEIk+K6wbtCXp7XfECHcYcDp/DNpS7YmHs2Yzwg1EOka2JVVVVUmSYmNjvT4eExOjr7/+usHnx8XFKeqnR11ZkC97SQISfcUfe/d9J6m6Xvn+/Xb169uq6RtkIqv1kz17a+RyHatT5nJJFcdi1a1biyC1yhqs1ldgnlDsK4TtAKioqAh2E0zHARTwFX3FP5df5n0Eu2vXapWXf9/ErTGPFfvJBXGGoqLkGdmWaqeSxLWuZJ9oP1ixr8AcoXqojfWHY30QExMjSaqsrPT6eFVVVYOj3gAiQ1MtersyOUr/Nbhu2X8Nri03E4v6/McR4wC8YWRb8mzpV1JSoquuuqrOYxUVFSovL9c111wThJYBCAVNvejtoQejNHSIS7s+ka6+yvygzaK+wMnMsKl3qix7fDqAxmNkW1JqaqokaePGjfUe++CDDyRJvXv3btI2AQgNwVr0dmVylEZkRTXJiDaL+gLL4bCp1zWMaAOoRdiWdN111+mSSy7Rm2++qcLCQk95VVWV/vrXv6p58+YaOnRoEFsIIFgOlNadgyvV/nygNDjtCTSrfz4ACDamkUhq3ry5nnrqKWVnZ+vXv/61MjMzFRMTo7Vr1+rAgQOaMmWKOnfuHOxmAgiCDonyuuitQ2Lw2hRIVv98ABBsjGz/nz59+mjZsmW69tpr9fbbb2vZsmW64IIL9Mwzz+juu+8OdvMABInVF71Z/fMBoYpFyZHDZhgGV9lPkbAlEVsvwVdW7StOp2HpRW9N/fms2k8QeFbsK/UXJdtYlBwAobr1H9NIAMAHDodNDkewW2Eeq38+IFQ0tCi5d6o1/yIPppEAAAA0GRYlRx7CNgAAQBNxL0o+HYuSrY2wDQAA0ERYlBx5mLMNIGCcTkN79tbogjiDXxwA0ABOGo0shG0AAfHj6vpjHPkNAGfBouTIwTQSAH7jyG8AALwjbAPwG6vrAQDwjrANwG+srgcAwDvCNgC/sboeAADvWCAJICDcq+srjsUqrnUlQRsAABG2AQSQw2FTt24tVF5O0AYAQGIaCQAAAGAawjYAAABgEsI2Ip7TaahgB3tC48zoJ+GLaxe+nE5DW/JruHYIa8zZRkT78dRDceohGkQ/CV9cu/DFqbSwCka2EbE49RC+oJ+EL65d+OLawUoI24hYnHoIX9BPwhfXLnxx7WAlhG1ELE49hC/oJ+GLaxe+uHawEsI2IhanHsIX9JPwxbULX1w7WInNMAwmQPmpvLw82E0wXXx8vGU/p9Np6EBp7YgJN3L/WbWv0E8Cqyn7CdcufDmdBqfShqlg/LkLxu+f+Pj4s9ZhNxJEPIfDJocj2K1AqKOfhC+uXfjiVNrwxC5AdTGNBAAAAAHBTjL1EbYBAAAQEOwkUx9hGwAAAAHBTjL1EbYBC+OY6vDFtQMQjthJpj4WSAIWxQKV8MW1AxDOMjNs6p0qdgH6P4xsAxbEApXwxbUDYAUOh029ronsEW03wjZgQSxQCV9cOwCwFsI2YEEsUAlfXDsAsBbCNmBBkbJAxYqLCCPl2gFApGCBJGBRVl+gYuVFhFa/dgAQSQjbgIVZ9ZjqhhYR9k61TjC16rUDgEjDNBIAYYdFhACAcEHYBhB2WEQIAAgXhG0AYYdFhACAcMGcbQBhqakXETqdBgsW4RP6SuA4nYb27K3RBXEG3yXCFmEbQNhqqkWEVt75BIFFXwmcH7/LY3yXCGs2wzCss0Ht/9m6davy8vL0ySef6LPPPlNVVZWGDh2qp59+usHnuFwuLVu2TMuXL1dJSYnsdrvS0tI0depUderU6YzvV15eHuBPEHri4+Mj4nPCf1brK06noWEjjToLMqOipJX/YNqKP6zWTyT6SiDxXeJcBOO+Eh8ff9Y6lpyz/dprr+mll17Srl275PBx2OvRRx/Vk08+KZfLpdGjR2vAgAHKy8vTsGHDtG/fPpNbDCBUsfMJfEVfCRy+S1iJJaeRjBo1SnfccYe6dOmiXbt2acSIEWesv3nzZq1YsUIpKSlauHChoqOjJUlDhgzR7bffrscee0yvvPJKUzQdQIhx73zy0xE2dj7BT9FXAofvElZiyZHtq6++WpdffrmaNWvmU/3c3FxJ0pQpUzxBW5Kuu+469e/fX1u3btUXX3xhSlsBhDZ2PoGv6CuBw3cJK7HkyHZjbdmyRXa7Xb169ar3WP/+/fX+++9r69at6ty5cxBaByDYOD4dvqKvBI77u6w4Fqu41pV8lwhbER+2q6urdfjwYSUlJXkdCXcvjiwuLm7ahgEIKRyfDl/RVwLH4bCpW7cWKi8naCN8RXzYrqyslCTFxMR4fdxdXlVV1eBrxMXFKeqnx9lZkC8rbgGJvgLf0E/gK/oKfBWKfSVkw3ZaWpqOHj3qc/3FixcrLS3NvAadQUVFRVDetylZcZsumIO+Al/QT+Ar+gp8Fapb/4Vs2M7MzNTx48d9rt+2bdtzep/Y2FhJDY9cu8sbGvkGACvg1MPw1ZTXjn4CNF7Ihu1HHnmkSd7HbrcrISFBBw4c0KlTp+rN23bP1T7bwTYAEK449TB8NeW1o58A58b6E4190Lt3b1VXV6ugoKDeYxs3bpQkpaamNnWzAMB0TuePAUqq3dd41mxDTqflDhe2nKa8dvQT4NwRtiUNHz5ckvTss8/q5MmTnvJNmzZp48aNSk1NZds/AJbESX3hqymvHf0EOHchO43EH9u2bdPKlSslSd9++60kafv27XrggQckSV26dNFdd93lqd+nTx9lZWUpNzdXQ4cO1YABA3TkyBG99dZbiomJ0WOPPdbknwEAmgIn9YWvprx29BPg3FlyZPvLL7/U6tWrtXr1ar333nv1yt5///16z3niiSf08MMPy2azacmSJfrPf/6jgQMHKjc3V5dddllTfwQAaBKc1Be+mvLa0U+Ac2czDIMJV36KhC2J2HoJvqKvhKem3mWCfhI4Vt+NhL4CX7H1HwAgZHHqYfhqymtHPwEaz5LTSAAAAIBQQNgGAAAATELYBgAAAExC2AYAWJrTaahgBwewhCOn09CW/BquHcIaCyQBAJbFEePh68drd4xrh7DGyDYAwJI4Yjx8ce1gJYRtAIAlccR4+OLawUoI2wAAS3IfMX46jhgPD1w7WAlhGwBgSRwxHr64drASFkgCACwrM8Om3qlq8iPG4T/3tas4Fqu41pVcO4QtwjYAwNI4Yjx8ORw2devWQuXlBG2EL6aRAAAAACYhbAMAAAAmIWwDAAAAJiFsA4APOPI7cDiCG0AkYYEkAJwFR34HDkdwA4g0jGwDwBlwbHTg8F0CiESEbQA4A46NDhy+SwCRiLANAGfAsdGBw3cJIBIRtgHgDDg2OnD4LgFEIhZIAsBZcOR34HAEN4BIQ9gGAB9w5HfgcAQ3gEjCNBIAAADAJIRtAAAAwCSEbQAAAojTRgGcjjnbAAAECKeNAvgpRrYBAAgATsgE4A1hGwCAAOCETADeELYBAAgATsgE4A1hGwCAAOCETADesEASAIAA4bRRhCqn06BfBglhGwCAAOK0UYQadskJLqaRAAAAWBS75AQfYRsAAMCi2CUn+AjbAAAAFsUuOcFH2AYAALAodskJPhZIAgAAWBi75AQXYRsAAMDi2CUneCwXtqurq/Xuu+8qLy9PRUVFOnTokKKjo3XFFVdo5MiRyszM9Po8l8ulZcuWafny5SopKZHdbldaWpqmTp2qTp06Ne2HAAAAgCVYbs72tm3bNGPGDG3evFnJyckaO3asBg0apN27d2vatGl68sknvT7v0Ucf1ZNPPimXy6XRo0drwIABysvL07Bhw7Rv374m/hQAAACwApthGJbaaLGoqEh79+7VTTfdpBYtWnjKv/nmGw0fPlylpaXKzc1Vjx49PI9t3rxZY8eOVUpKihYuXKjo6GhJ0qZNm3T77bcrJSVFr7zySoPvWV5ebt4HChHx8fER8TnhP/oKfEE/ga/oK/BVMPpKfHz8WetYbmT7iiuu0C233FInaEtS27ZtNWLECEnS1q1b6zyWm5srSZoyZYonaEvSddddp/79+2vr1q364osvTG45AAAArMZyYftMmjevnaLerFmzOuVbtmyR3W5Xr1696j2nf//+kuoHdAAAAOBsLLdAsiGnTp3SP//5T9lsNvXt29dTXl1drcOHDyspKaleCJfkWRxZXFzc4GvHxcUp6qc7xluQL/9UAkj0FfiGfgJf0Vfgq1DsKxETtp977jnt2bNHv/rVr5SUlOQpr6yslCTFxMR4fZ67vKqqqsHXrqioCGBLQxNz5uAr+gp8QT+Br+gr8FWoztkO2bCdlpamo0eP+lx/8eLFSktL8/rY8uXLNX/+fF155ZV66KGHAtRCAAAA4MxCNmxnZmbq+PHjPtdv27at1/LXXntNjz76qJKSkvTSSy/p/PPPr/N4bGyspIZHrt3lDY18AwAax+k0tGdvjS6IMzjJLsw4nQanEAKNFLJh+5FHHvH7NVauXKlHHnlEl112mV5++WWvQ/12u10JCQk6cOCATp06VW/etnuuNgfbAID/3lxjaNZsQy7XMUVFSTOm1R4ljdD347UT1w5oBMuu6lu5cqUefvhhdenSRS+//LIuvPDCBuv27t1b1dXVKigoqPfYxo0bJUmpqammtRUAIoHT+WNYkySXS5o125DTaanjHiyJawecO0uG7dzc3DpBu02bNmesP3z4cEnSs88+q5MnT3rKN23apI0bNyo1NVWdO3c2tc0AYHUHSuUJa24uV205QhvXDjh3ITuN5Fxt2rRJjzzyiAzDUEpKil599dV6dZKTk5Wenu75uU+fPsrKylJubq6GDh2qAQMG6MiRI3rrrbcUExOjxx57rAk/AQBYU4fE2ukHp4e2qKjacoQ2rh1w7iwXtg8dOiT3CfTLly/3Wmfo0KF1wrYkPfHEE+rWrZuWL1+uJUuWyG63a+DAgZo6dSqj2gAQAA6HTTOm6Sfzfm0stAsDXDvg3NkMdzLFOYuE/T/Z5xS+oq/gbJxOQxXHYhXXupKwFmaCsRsJ9xT4in22AQBQbUjr1q2FyssJ2uHG4bDJ4Qh2K4DwYskFkgAAAEAoIGwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmaR7sBgAAAHjjdBras7dGF8QZcjhswW4OcE4I2wAAIOS8ucbQrNmGXK5jioqSZkyTMjMI3Ag/TCMBAAAhxel0B+3an10uadZsQ06nEdyGAefAkiPbL7zwgjZv3qz9+/ervLxcrVq1UmJiom655RaNHDlSrVq1qvccl8ulZcuWafny5SopKZHdbldaWpqmTp2qTp06Nf2HAAAgQh0olSdou7lcteUOR3DaBJwrS45s/+Mf/1BFRYX69u2rMWPGKCMjQydPntTTTz+tkSNH6rvvvqv3nEcffVRPPvmkXC6XRo8erQEDBigvL0/Dhg3Tvn37gvApAACITB0SpaifJJSoqNpyINxYcmT77bffVsuWLeuVz5gxQ6+//rpWrVqlUaNGeco3b96sFStWKCUlRQsXLlR0dLQkaciQIbr99tv12GOP6ZVXXmmy9gMAEMkcDptmTJNnKkntnG0biyQRliwZtr0FbUkaPHiwXn/9dZWUlNQpz83NlSRNmTLFE7Ql6brrrlP//v31/vvv64svvlDnzp3NazQAAPDIzLCpd6pUcSxWca0rCdoIW5acRtKQ9957T5J0+eWX1ynfsmWL7Ha7evXqVe85/fv3lyRt3brV/AYCAAAPh8Om3qktCNoIa5Yc2XZbtGiRKisrdezYMRUUFOiTTz5R//79NWTIEE+d6upqHT58WElJSWrWrFm913AvjiwuLm6aRgMAAMAyLB22Fy9erNLSUs/Pt956qx577DG1aNHCU1ZZWSlJiomJ8foa7vKqqqoG3ycuLk5RP13JYUHx8fHBbgLCBH0FvqCfwFf0FfgqFPtKyIbttLQ0HT161Of6ixcvVlpaWp2yvLw8SdLhw4e1ZcsWPfPMMxo+fLj+/ve/q127dgFra0VFRcBeK1TFx8ervLw82M1AGKCvwBf0E/iKvgJfBaOv+BLuQzZsZ2Zm6vjx4z7Xb9u2bYOPJSQkKDMzU5deeqmysrL09NNP69lnn5UkxcbGSmp45Npd3tDINwAAANCQkA3bjzzySMBfs0ePHoqLi1N+fr6nzG63KyEhQQcOHNCpU6fqzdt2z9XmYBsAAAA0lvUnGp/m+PHjqqysrBeoe/furerqahUUFNR7zsaNGyVJqampTdJGAAAAWIflwnZpaakOHDhQr7ympkZ//OMf5XK5dP3119d5bPjw4ZKkZ599VidPnvSUb9q0SRs3blRqaip7bAMAAKDRQnYaybkqLCzU5MmTlZKSoo4dOyo+Pl7ffPONNm3apEOHDqlz586aOnVqnef06dNHWVlZys3N1dChQzVgwAAdOXJEb731lmJiYvTYY48F58MAAAAgrNkMwzCC3YhAOnjwoF5++WVt3bpVpaWlqqyslN1uV9euXZWenq5Ro0bJbrfXe57L5dLSpUu1fPlylZSUyG63Ky0tTVOnTj3rqHYkrJJmNTh8RV+BL+gn8BV9Bb4K1d1ILBe2gyESbgLc7OAr+gp8QT+Br+gr8FWohm3LzdkGAAAAQgVhGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJu5EAAAAAJmFkGwAAADAJYRsAAAAwCWEbAAAAMAlhGwAAADAJYRsAAAAwSfNgNwCh64YbblBpaanXx0aMGKEnnniiiVuEYHr99de1fft2ffLJJ9qzZ49qamo0c+ZM3XbbbV7rV1VVae7cuVq7dq0OHz6shIQEDRo0SJMnT1ZMTEwTtx5NpTH9ZO7cuZo3b57X14mOjtauXbvMbi6CpKysTG+//bY2bNigzz//XN98843i4uLUq1cvZWdnq2fPnvWewz0lMjW2r4TifYWwjTOKjY3V2LFj65VfddVVQWgNgum5555TaWmp4uPj5XA4GvyLmCRVV1dr9OjRKiwsVL9+/ZSRkaGioiItWrRIW7Zs0bJly2S325uw9WgqjeknbkOHDlViYmKdsmbNmpnVRISAJUuW6MUXX9Sll16qvn37qk2bNiopKdG6deu0bt06zZ49WzfffLOnPveUyNXYvuIWUvcVA2jAwIEDjYEDBwa7GQgRH3zwgXHgwAHDMAxj/vz5RlJSkvHaa695rfvcc88ZSUlJxqxZs7yWP/fcc6a3F8HRmH4yZ84cIykpydi8eXNTNhEh4J133jG2bt1ar3zr1q1G9+7djd69exsnTpzwlHNPiVyN7SuheF9hzjYAn/Tt27feKIE3hmEoNzdXdrtdEydOrPPYhAkTFBcXp5UrV8rgPC1L8rWfILINGjRIKSkp9cpTUlKUlpamo0ePavfu3ZK4p0S6xvSVUMU0EpzRyZMntXr1apWVlal169bq1auXrrjiimA3CyGsuLhYTqdT/fv3r/fPui1btlRKSor+/e9/q6SkRJ06dQpOIxFStm3bpp07d6pZs2bq0qWL+vbtq+jo6GA3C0HSvHnzOv/lnoKG/LSvnC6U7iuEbZzR4cOH9cADD9Qp+/nPf65Zs2bpwgsvDFKrEMpKSkokqcFfeh07dvTU4xcjJGnOnDl1fk5ISNCf/vQn9evXL0gtQrAcPHhQH374oRISEpSUlCSJewq889ZXThdK9xWmkaBBt912m5YsWaJNmzZp+/btWrFiha6//nq9//77uueee/gnO3hVWVkpSQ3uDuAud9dD5EpOTtaf/vQn5eXlaefOnVq7dq1+85vfqLKyUnfffbeKioqC3UQ0oZqaGs2YMUMnT57Ufffd51nMxj0FP9VQX5FC875C2EaDJk2apN69e+vCCy9UTEyMevbsqfnz5+vaa6/Vjh079N577wW7iQDCWHp6uoYMGaLExES1bNlSHTt21D333KOHHnpIJ06c0F//+tdgNxFNxOVy6Xe/+522bt2q4cOHa8iQIcFuEkLU2fpKKN5XCNtolKioKM9+uQUFBUFuDUJRbGyspNo9cb1xl7vrAT81ZMgQNW/enHtMhDAMQw8//LDeeOMN3XrrrXr88cfrPM49BW5n6ytnEsz7CmEbjRYfHy9J+u6774LcEoQi9/zJ4uJir4+751+66wE/FR0drfPPP1/ff/99sJsCk7lHKV977TVlZmbq6aefVlRU3WjCPQWSb33lTIJ5XyFso9F27twpSWzvBa86deokh8OhgoICVVdX13nsxIkT2rZtmxwOB78Y0aDi4mJVVFRwj7E4l8ulhx56SKtWrdLNN9+sWbNmeT10hHsKfO0rZxLM+wphG17t27dPx44dq1e+bds2LVy4UNHR0Ro0aFAQWoZQZ7PZlJWVperqauXk5NR5bP78+aqoqFBWVpZsNluQWohQUFVV5XWhUkVFhR566CFJUkZGRlM3C03k9PB000036ZlnnmkwPHFPiWyN6Suhel+xGWwpAS/mzp2rBQsW6LrrrlNiYqKio6O1Z88effDBB4qKitLjjz+urKysYDcTTSg3N1fbt2+XJO3Zs0effvqpevXq5RlNSk9PV3p6uqTao5V//etfe45W7t69u4qKirRhwwYlJydztLKF+dpPDhw4oBtvvFFXXXWVkpKS1KZNG5WVlWnDhg06evSo+vXrp+eff579ti1q7ty5mjdvnux2u8aMGeN1n+T09HQlJydL4p4SyRrTV0L1vsI+2/AqLS1N+/fv12effab8/HydPHlSbdq00c0336xx48apR48ewW4imtj27du1evXqOmUFBQWexSaJiYmesG2327VkyRLNmzdP77zzjvLz89W2bVuNGzdOkyZN4peihfnaTy644AKNGjVKH330kdavX6/Kykq1atVKSUlJuvXWW5WVldXofyZG+CgtLZVUG6Kff/55r3USExM9YZt7SuRqTF8J1fsKI9sAAACASZizDQAAAJiEsA0AAACYhLANAAAAmISwDQAAAJiEsA0AAACYhLANAAAAmISwDQAAAJiEsA0AAACYhLANAAAAmISwDQAIC926ddMNN9wQ7GYAQKMQtgEAAACTELYBAAAAkxC2AQAAAJM0D3YDAACBsX//fj3//PPasWOHysrKdP755+uiiy5S7969deedd8rhcGjLli0aM2aMhg4dqt/+9rf685//rA0bNqiqqkpdu3bV2LFjNWTIEK+v/+233+rFF19UXl6eDh48qPPOO089e/bUhAkTlJqa6vU5u3fv1gsvvKD8/HyVl5frggsu0M9//nNNnDhRHTp0qFe/urpaOTk5WrNmjY4cOaLExESNGDFC48aNC+A3BQBNh7ANABbw6aef6te//rVOnDihHj16qEePHjp+/Li++uorLV68WOnp6XI4HJ76R48e1YgRI3Ty5En17t1bx44d05YtW3T//ffrwIEDmjRpUp3X379/v26//XaVlZXp0ksv1YABA3T06FFt3rxZH3zwgWbNmqVbbrmlznPeeecdTZs2TTU1NerevbuuueYaffXVV1q1apXy8vL0yiuv6PLLL/fUP3nypMaPH68dO3YoPj5eAwcO1PHjxzV79mx9+eWX5n6BAGASwjYAWMCSJUv0/fffa+7cuRo0aFCdx/bv36/Y2Ng6ZevXr1e/fv00b9482e12SdLOnTs1duxY5eTk6MYbb1RycrIk6dSpU5oyZYrKysr00EMP6X/+539ks9kkSZ999pluv/12/f73v1ffvn3Vpk0bSdJXX32l+++/X+edd54WLlxYZ+T7n//8p+6//349+OCDWrlypad84cKF2rFjh3r06KGXXnrJ0+ZPP/1UY8aMCfA3BgBNgznbAGABR44ckST16dOn3mNdu3atM6otSTabTQ8//LAnaEtSjx49NGrUKLlcLr366que8vXr12vPnj3KzMzUmDFjPEFbkq688krdc889qq6u1htvvOEpX7x4sb777jtNnz693hSTIUOGKD09Xbt27dKnn37qKXe/54MPPljnLwfdu3fXqFGjGvV9AECoIGwDgAV0795dkjRjxgzt3LlTLpfrjPWvvPJKdenSpV55ZmamJGn79u2esg8++ECSdOONN3p9rWuvvVaStGvXLk/Zhx9+2KjnHDx4UIcOHdJFF12kXr161aufkZFxxs8DAKGKaSQAYAHZ2dnavn271q9fr/Xr1ys2NlY9e/bUL37xCw0dOlQxMTF16rdv397r6yQmJkqSnE6np6y0tFSSNHXqVE2dOrXBNpSXl9d7Tr9+/c7Ybvdz3O/XULsaKgeAUEfYBgALiImJ0eLFiz2BOz8/X5s2bdLGjRs1f/58LVu2TJdeeuk5vfapU6ckSddff71nTrY3p4+Unzp1SjabrcGdTdzcCyQNwzhjvdOnrgBAOCFsA4BF2Gw2paSkKCUlRVLtVn1/+MMf9Oabb+rPf/6znn32WU/dgwcPen0N94j06XO827VrJ0kaOXJkg9NCfqpdu3b68ssv9fDDD9cbVffG/X5naxcAhBvmbAOARV144YWeLfz27NlT57HCwkJ98cUX9Z6zZs0aSaozb7pv376SpHXr1vn83tddd12jnpOYmKh27dqprKxMO3bsqPf4W2+95fN7A0AoIWwDgAW8+uqr+uqrr+qVb9iwQZJ08cUX1yl3uVx66qmn9N1333nKPvnkEy1dulRRUVEaMWKEp3zw4MHq0qWLVq9erRdeeEE1NTV1XuvkyZNau3atdu/e7SkbP368zjvvPM2cOVN5eXn12nX06FEtXbpU33//vafM/Z5/+tOfVFVV5SkvLCzU0qVLffoeACDU2IyzTZQDAIS8//7v/1ZRUZEuu+wyde3aVc2aNdMXX3yhwsJCnXfeeVq0aJGuueYazwmSAwcO1O7du1VTU6OUlBRVVlZqy5Ytqqmp0d13360pU6bUef39+/crOztbBw8eVEJCgrp166aYmBh9/fXX+vzzz3Xs2DHl5OQoPT3d85y1a9dq+vTp+v7779W5c2d17dpVhmHo4MGD2rdvn2pqarR161a1bt1aUm1oHz16tD7++GPFx8crLS1Nx48f1+bNmzVs2DC9+uqrSkxM9BreASBUEbYBwALy8vK0bt067dy5U2VlZaqpqdFFF12ktLQ03XHHHerYsaMk1TmuferUqfrf//1fbdy4UVVVVerSpYvGjh2r2267zet7VFRUaMmSJXr33Xf15ZdfyjAMJSQk6LLLLtMvf/lLDR48WOeff36d5xQXF+ull17Shx9+qLKyMrVs2VIOh0M9e/bU4MGDNWDAgDqLH48fP6558+ZpzZo1+vbbb5WYmKisrCyNHz9eycnJhG0AYYewDQAR5PSw/fTTTwe7OQBgeczZBgAAAExC2AYAAABMQtgGAAAATMKcbQAAAMAkjGwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACYhbAMAAAAmIWwDAAAAJiFsAwAAACb5//V0Wwko+6l/AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Load the data:\n",
"# car braking distances in feet paired with speeds in km/h\n",
"# see cars.info() for details\n",
"cars = pd.read_csv(\"../data/cars.csv\", index_col=0)\n",
"\n",
"\n",
"# fit a linear regression of distance on speed\n",
"def model(speed, dist_):\n",
" mu = numpyro.param(\"a\", 0.0) + numpyro.param(\"b\", 1.0) * speed\n",
" numpyro.sample(\"dist\", dist.Normal(mu, 1), obs=dist_)\n",
"\n",
"\n",
"svi = SVI(\n",
" model,\n",
" lambda speed, dist_: None,\n",
" optim=optim.Adam(1),\n",
" loss=Trace_ELBO(),\n",
" speed=cars.speed.values,\n",
" dist_=cars.dist.values,\n",
")\n",
"svi_result = svi.run(random.PRNGKey(0), 1000)\n",
"params = svi_result.params\n",
"\n",
"# estimated coefficients from the model\n",
"print(params)\n",
"\n",
"# plot residuals against speed\n",
"resid = cars.dist - (params[\"a\"] + params[\"b\"] * cars.speed.values)\n",
"az.plot_pair({\"speed\": cars.speed, \"resid\": resid})\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"ExecuteTime": {
"end_time": "2017-12-07T01:31:33.798552Z",
"start_time": "2017-12-07T01:31:33.508571Z"
}
},
"source": [
"### Code 0.5"
]
},
{
"cell_type": "markdown",
"metadata": {
"ExecuteTime": {
"end_time": "2017-12-07T01:31:33.798552Z",
"start_time": "2017-12-07T01:31:33.508571Z"
}
},
"source": [
"```sh\n",
"pip install numpyro arviz daft networkx\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"navigation": true
},
"source": [
"\n",
"< []() | [Chapter 1. The Golem of Prague](01-the-golem-of-prague.html) >"
]
}
],
"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.6"
},
"nikola": {
"title": "Preface"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}
![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\"\n",){kind=icon}
\n",
" \n",
" \n",
"