{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "(pymc_pytensor)=\n", "\n", "# PyMC and PyTensor\n", "\n", "**Authors:** [Ricardo Vieira](https://github.com/ricardoV94) and [Juan Orduz](https://juanitorduz.github.io/)\n", "\n", "In this notebook we want to give an introduction of how PyMC models translate to PyTensor graphs. The purpose is not to give a detailed description of all [`pytensor`](https://github.com/pymc-devs/pytensor)'s capabilities but rather focus on the main concepts to understand its connection with PyMC. For a more detailed description of the project please refer to the official documentation." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "heading_collapsed": true }, "source": [ "## Prepare Notebook\n", "\n", "First import the required libraries." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "hidden": true }, "outputs": [], "source": [ "import pytensor\n", "import pytensor.tensor as pt\n", "import pymc as pm\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy.stats" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction to PyTensor\n", "\n", "We start by looking into `pytensor`. According to their documentation\n", "\n", "> PyTensor is a Python library that allows one to define, optimize, and efficiently evaluate mathematical expressions involving multi-dimensional arrays." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "![pytensor logo](https://raw.githubusercontent.com/pymc-devs/pytensor/main/doc/images/PyTensor_RGB.svg)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### A simple example\n", "\n", "To begin, we define some pytensor tensors and show how to perform some basic operations." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "x type: TensorType(float64, ())\n", "x name = x\n", "---\n", "y type: TensorType(float64, (?,))\n", "y name = y\n", "\n" ] } ], "source": [ "x = pt.scalar(name=\"x\")\n", "y = pt.vector(name=\"y\")\n", "\n", "print(\n", " f\"\"\"\n", "x type: {x.type}\n", "x name = {x.name}\n", "---\n", "y type: {y.type}\n", "y name = {y.name}\n", "\"\"\"\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have defined the `x` and `y` tensors, we can create a new one by adding them together." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "z = x + y\n", "z.name = \"x + y\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "To make the computation a bit more complex let us take the logarithm of the resulting tensor." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "w = pt.log(z)\n", "w.name = \"log(x + y)\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can use the {func}`~pytensor.dprint` function to print the computational graph of any given tensor." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elemwise{log,no_inplace} [id A] 'log(x + y)'\n", " |Elemwise{add,no_inplace} [id B] 'x + y'\n", " |InplaceDimShuffle{x} [id C]\n", " | |x [id D]\n", " |y [id E]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pytensor.dprint(w)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Note that this graph does not do any computation (yet!). It is simply defining the sequence of steps to be done. We can use {func}`~pytensor.function` to define a callable object so that we can push values trough the graph." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "f = pytensor.function(inputs=[x, y], outputs=w)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now that the graph is compiled, we can push some concrete values:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0., 1.])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x=0, y=[1, np.e])" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ ":::{tip}\n", "Sometimes we just want to debug, we can use {meth}`~pytensor.graph.basic.Variable.eval` for that:\n", ":::" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0., 1.])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w.eval({x: 0, y: [1, np.e]})" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "You can set intermediate values as well" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0., 1.])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w.eval({z: [1, np.e]})" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### PyTensor is clever!\n", "\n", "One of the most important features of `pytensor` is that it can automatically optimize the mathematical operations inside a graph. Let's consider a simple example:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elemwise{true_div,no_inplace} [id A] 'a / b'\n", " |a [id B]\n", " |b [id C]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = pt.scalar(name=\"a\")\n", "b = pt.scalar(name=\"b\")\n", "\n", "c = a / b\n", "c.name = \"a / b\"\n", "\n", "pytensor.dprint(c)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now let us multiply `b` times `c`. This should result in simply `a`." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elemwise{mul,no_inplace} [id A] 'b * c'\n", " |b [id B]\n", " |Elemwise{true_div,no_inplace} [id C] 'a / b'\n", " |a [id D]\n", " |b [id B]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d = b * c\n", "d.name = \"b * c\"\n", "\n", "pytensor.dprint(d)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The graph shows the full computation, but once we compile it the operation becomes the identity on `a` as expected." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DeepCopyOp [id A] 'a' 0\n", " |a [id B]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "g = pytensor.function(inputs=[a, b], outputs=d)\n", "\n", "pytensor.dprint(g)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### What is in an PyTensor graph?\n", "\n", "The following diagram shows the basic structure of an `pytensor` graph.\n", "\n", "![pytensor graph](https://raw.githubusercontent.com/pymc-devs/pytensor/main/doc/tutorial/apply.png)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can can make these concepts more tangible by explicitly indicating them in the first example from the section above. Let us compute the graph components for the tensor `z`. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "z type: TensorType(float64, (?,))\n", "z name = x + y\n", "z owner = Elemwise{add,no_inplace}(InplaceDimShuffle{x}.0, y)\n", "z owner inputs = [InplaceDimShuffle{x}.0, y]\n", "z owner op = Elemwise{add,no_inplace}\n", "z owner output = [x + y]\n", "\n" ] } ], "source": [ "print(\n", " f\"\"\"\n", "z type: {z.type}\n", "z name = {z.name}\n", "z owner = {z.owner}\n", "z owner inputs = {z.owner.inputs}\n", "z owner op = {z.owner.op}\n", "z owner output = {z.owner.outputs}\n", "\"\"\"\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The following code snippet helps us understand these concepts by going through the computational graph of `w`. The actual code is not as important here, the focus is on the outputs." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---\n", "Checking variable log(x + y) of type TensorType(float64, (?,))\n", " > Op is Elemwise{log,no_inplace}\n", " > Input 0 is x + y\n", "---\n", "Checking variable x + y of type TensorType(float64, (?,))\n", " > Op is Elemwise{add,no_inplace}\n", " > Input 0 is InplaceDimShuffle{x}.0\n", " > Input 1 is y\n", "---\n", "Checking variable InplaceDimShuffle{x}.0 of type TensorType(float64, (1,))\n", " > Op is InplaceDimShuffle{x}\n", " > Input 0 is x\n", "---\n", "Checking variable y of type TensorType(float64, (?,))\n", " > y is a root variable\n", "---\n", "Checking variable x of type TensorType(float64, ())\n", " > x is a root variable\n" ] } ], "source": [ "# start from the top\n", "stack = [w]\n", "\n", "while stack:\n", " print(\"---\")\n", " var = stack.pop(0)\n", " print(f\"Checking variable {var} of type {var.type}\")\n", " # check variable is not a root variable\n", " if var.owner is not None:\n", " print(f\" > Op is {var.owner.op}\")\n", " # loop over the inputs\n", " for i, input in enumerate(var.owner.inputs):\n", " print(f\" > Input {i} is {input}\")\n", " stack.append(input)\n", " else:\n", " print(f\" > {var} is a root variable\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Note that this is very similar to the output of {func}`~pytensor.dprint` function introduced above." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elemwise{log,no_inplace} [id A] 'log(x + y)'\n", " |Elemwise{add,no_inplace} [id B] 'x + y'\n", " |InplaceDimShuffle{x} [id C]\n", " | |x [id D]\n", " |y [id E]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pytensor.dprint(w)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Graph manipulation 101\n", "\n", "Another interesting feature of PyTensor is the ability to manipulate the computational graph, something that is not possible with TensorFlow or PyTorch. Here we continue with the example above in order to illustrate the main idea around this technique." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[x, y]" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# get input tensors\n", "list(pytensor.graph.graph_inputs(graphs=[w]))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "As a simple example, let's add an {func}`~pytensor.tensor.exp` before the {func}`~pytensor.tensor.log` (to get the identity function)." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "parent_of_w = w.owner.inputs[0] # get z tensor\n", "new_parent_of_w = pt.exp(parent_of_w) # modify the parent of w\n", "new_parent_of_w.name = \"exp(x + y)\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Note that the graph of `w` has actually not changed:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elemwise{log,no_inplace} [id A] 'log(x + y)'\n", " |Elemwise{add,no_inplace} [id B] 'x + y'\n", " |InplaceDimShuffle{x} [id C]\n", " | |x [id D]\n", " |y [id E]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pytensor.dprint(w)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "To modify the graph we need to use the {func}`~pytensor.clone_replace` function, which *returns a copy of the initial subgraph with the corresponding substitutions.*" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elemwise{log,no_inplace} [id A] 'log(exp(x + y))'\n", " |Elemwise{exp,no_inplace} [id B] 'exp(x + y)'\n", " |Elemwise{add,no_inplace} [id C] 'x + y'\n", " |InplaceDimShuffle{x} [id D]\n", " | |x [id E]\n", " |y [id F]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "new_w = pytensor.clone_replace(output=[w], replace={parent_of_w: new_parent_of_w})[0]\n", "new_w.name = \"log(exp(x + y))\"\n", "pytensor.dprint(new_w)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can test the modified graph by passing some input to the new graph." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1. , 2.71828183])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "new_w.eval({x: 0, y: [1, np.e]})" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the new graph is just the identity function." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ ":::{note}\n", "Again, note that `pytensor` is clever enough to omit the `exp` and `log` once we compile the function.\n", ":::" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Elemwise{add,no_inplace} [id A] 'x + y' 1\n", " |InplaceDimShuffle{x} [id B] 0\n", " | |x [id C]\n", " |y [id D]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f = pytensor.function(inputs=[x, y], outputs=new_w)\n", "\n", "pytensor.dprint(f)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1. , 2.71828183])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x=0, y=[1, np.e])" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### PyTensor RandomVariables\n", "\n", "Now that we have seen pytensor's basics we want to move in the direction of random variables.\n", "\n", "How do we generate random numbers in [`numpy`](https://github.com/numpy/numpy)? To illustrate it we can sample from a normal distribution:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAIQCAYAAACbhEYhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCrklEQVR4nO3deVRV9f7/8dcBZHAAhBAkFdDKqbJyQNJyotDMKa28maKZ5pxSqXRT05tiZWmaQ/YttZs2aVpmYeZ4yzG9ZVlZJg7pBTUFFBMU9u+PFufnEVRAZMPH52OtvZb7s4fz3vts4OXnfM7eDsuyLAEAAAAGcLO7AAAAAKC4EG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQboES5nA49Pzzz5foa6akpKhbt24KDAyUw+HQtGnTSvT1cWX27dsnh8Oh+fPnF2n7+fPny+FwaN++fc62li1bqmXLlsVS3+VceM0///zzcjgcOnbsWIm8fnh4uHr37l0ir3W1XOk1AFxLCLcok3744Qd169ZNYWFh8vb21vXXX6977rlHM2bMsLu0UmnEiBFauXKl4uPj9e9//1tt27a1uySUQRs3btTzzz+v1NRUu0vJozTXBqBkedhdAFBYGzduVKtWrVSjRg3169dPISEhOnjwoDZv3qzXXntNQ4cOtbvEUmfNmjXq1KmTnn76abtLQSnx5ZdfFnqbjRs3avz48erdu7f8/f0LvN1ff/0lD4+r++fmUrXt3r1bbm5luy8nLCxMf/31l8qVK2d3KUCpR7hFmTNx4kT5+flp27Ztef6IHTlyxJ6iSrkjR44UKIxkZGSoQoUKV7+gMu7MmTPy9PQs04HJ09Pzqu4/JydHWVlZ8vb2lre391V9rcvx8vKy9fWLg8PhsP08AmVF2f3NjGvW77//rvr16+cb1qpUqeIyP2/ePLVu3VpVqlSRl5eX6tWrp9mzZ+fZLjw8XPfff7/WrVunRo0aycfHR7fccovWrVsnSfr44491yy23yNvbWw0bNtR///tfl+179+6tihUrau/evYqJiVGFChUUGhqqCRMmyLKsyx7ToUOH9Nhjjyk4OFheXl6qX7++3n777TzrzZgxQ/Xr11f58uVVuXJlNWrUSIsWLbrofnPHWlqWpZkzZ8rhcMjhcLgsW79+vQYNGqQqVaqoWrVqzm1nzZql+vXry8vLS6GhoRo8eHCej3xbtmypm2++WTt37lSLFi1Uvnx53XDDDVq8eLEkaf369YqMjJSPj49q166tr7766rLnIisrS2PHjlXDhg3l5+enChUq6K677tLatWsvu630/9/Lr7/+Wk2aNJG3t7dq1qypd955J8+6e/fu1YMPPqiAgACVL19eTZs21YoVK1zWWbdunRwOh95//30999xzuv7661W+fHmlp6c73/cDBw7o/vvvV8WKFXX99ddr5syZkv4ePtO6dWtVqFBBYWFhed6r48eP6+mnn9Ytt9yiihUrytfXV+3atdP3339foGPNz65du9S6dWv5+PioWrVqeuGFF5STk5NnvfzG3F7q+nr++ef1zDPPSJIiIiKc11LuOF6Hw6EhQ4Zo4cKFzusmMTHRuSy/cebHjh3TQw89JF9fXwUGBurJJ5/UmTNnnMsvNc70/H1errb8xtwW5r3/8MMPNXHiRFWrVk3e3t5q06aN9uzZk6emC/Xu3Vvh4eF52nPHHJ9v1apVat68ufz9/VWxYkXVrl1bzz777CXPRe71d+jQIXXu3FkVK1ZUUFCQnn76aWVnZ7vs/88//1TPnj3l6+srf39/xcbG6vvvvy/QON7c3xXffPON4uLiFBQUpAoVKqhLly46evSoy7oXe68vfA9y9/n1119r2LBhCgoKkr+/v5544gllZWUpNTVVvXr1UuXKlVW5cmWNHDnS5Xdp7vmYMmWKpk6dqrCwMPn4+KhFixb68ccfnevNmzdPDocjz+9sSZo0aZLc3d116NChSx4/yh56blHmhIWFadOmTfrxxx918803X3Ld2bNnq379+urYsaM8PDy0fPlyDRo0SDk5ORo8eLDLunv27NEjjzyiJ554Qo8++qimTJmiDh06aM6cOXr22Wc1aNAgSVJCQoIeeuihPB91Zmdnq23btmratKleeuklJSYmaty4cTp37pwmTJhw0RpTUlLUtGlTZzgICgrSF198ob59+yo9PV3Dhw+XJL355psaNmyYunXr5gwBO3fu1JYtW/TII4/ku++7775b//73v9WzZ0/dc8896tWrV551Bg0apKCgII0dO1YZGRmS/v7jO378eEVHR2vgwIHavXu3Zs+erW3btumbb75x+Wj0xIkTuv/++9W9e3c9+OCDmj17trp3766FCxdq+PDhGjBggB555BG9/PLL6tatmw4ePKhKlSpd9Hykp6fr//7v//SPf/xD/fr108mTJ/XWW28pJiZGW7du1W233XbRbXPt2bNH3bp1U9++fRUbG6u3335bvXv3VsOGDVW/fn3neb/zzjt1+vRpDRs2TIGBgVqwYIE6duyoxYsXq0uXLi77/Ne//iVPT089/fTTyszMdPZ8Zmdnq127drr77rv10ksvaeHChRoyZIgqVKigf/7zn+rRo4ceeOABzZkzR7169VJUVJQiIiIk/R2wli1bpgcffFARERFKSUnRG2+8oRYtWuinn35SaGjoZY/1fMnJyWrVqpXOnTun0aNHq0KFCpo7d658fHwuu+3lrq8HHnhAv/76q9577z1NnTpV1113nSQpKCjIuY81a9boww8/1JAhQ3TdddflG+zO99BDDyk8PFwJCQnavHmzpk+frhMnTuT7H5FLKUht5yvsez958mS5ubnp6aefVlpaml566SX16NFDW7ZsKVSdF7Nr1y7df//9uvXWWzVhwgR5eXlpz549+uabby67bXZ2tmJiYhQZGakpU6boq6++0iuvvKJatWpp4MCBkv7uRe/QoYO2bt2qgQMHqk6dOvrkk08UGxtbqDqHDh2qypUra9y4cdq3b5+mTZumIUOG6IMPPijScefuMyQkROPHj9fmzZs1d+5c+fv7a+PGjapRo4YmTZqkzz//XC+//LJuvvnmPL/D3nnnHZ08eVKDBw/WmTNn9Nprr6l169b64YcfFBwcrG7dumnw4MFauHChbr/9dpdtFy5cqJYtW+r6668vcv0opSygjPnyyy8td3d3y93d3YqKirJGjhxprVy50srKysqz7unTp/O0xcTEWDVr1nRpCwsLsyRZGzdudLatXLnSkmT5+PhY+/fvd7a/8cYbliRr7dq1zrbY2FhLkjV06FBnW05OjtW+fXvL09PTOnr0qLNdkjVu3DjnfN++fa2qVatax44dc6mpe/fulp+fn/MYOnXqZNWvX/8yZyd/kqzBgwe7tM2bN8+SZDVv3tw6d+6cs/3IkSOWp6ende+991rZ2dnO9tdff92SZL399tvOthYtWliSrEWLFjnbfvnlF0uS5ebmZm3evNnZnns+582bd8laz507Z2VmZrq0nThxwgoODrYee+yxyx5r7nu5YcMGl2Py8vKynnrqKWfb8OHDLUnWf/7zH2fbyZMnrYiICCs8PNx57GvXrrUkWTVr1sxzPeW+75MmTXKp1cfHx3I4HNb777+f57yc/96fOXPG5RxblmUlJSVZXl5e1oQJE1zaCnLuco9py5YtLsfu5+dnSbKSkpKc7S1atLBatGjhnC/I9fXyyy/n2U+u3Pd8165d+S47/7jHjRtnSbI6duzost6gQYMsSdb3339vWdalj/vCfV6qtrCwMCs2NtY5X9j3vm7dui7X5GuvvWZJsn744Yc8r3W+2NhYKywsLE977vHnmjp1qiXJ5ffEhfI7F7nX3/nXimVZ1u233241bNjQOb9kyRJLkjVt2jRnW3Z2ttW6desCXVe5vyuio6OtnJwcZ/uIESMsd3d3KzU11dl24fuS68L3IHefMTExLvuMioqyHA6HNWDAAGfbuXPnrGrVqrlcr7nnw8fHx/rjjz+c7Vu2bLEkWSNGjHC2/eMf/7BCQ0NdftZ27NhRoGNH2cSwBJQ599xzjzZt2qSOHTvq+++/10svvaSYmBhdf/31+vTTT13WPb/HKi0tTceOHVOLFi20d+9epaWluaxbr149RUVFOecjIyMlSa1bt1aNGjXytO/duzdPbUOGDHH+O7cnNisr66Ifx1uWpSVLlqhDhw6yLEvHjh1zTjExMUpLS9OOHTskSf7+/vrjjz+0bdu2Ap2ngurXr5/c3d2d81999ZWysrI0fPhwl57pfv36ydfXN89HtxUrVlT37t2d87Vr15a/v7/q1q3rPFfSpc/b+dzd3Z29ojk5OTp+/LjOnTunRo0aOc/F5dSrV0933XWXcz4oKEi1a9d2ee3PP/9cTZo0UfPmzV2OpX///tq3b59++uknl33GxsZetAf08ccfd/7b399ftWvXVoUKFfTQQw8523PPy/k1eHl5Oc9xdna2/vzzT+dH0gU91vN9/vnnatq0qZo0aeJy7D169LjstsVxfbVo0UL16tUr8PoXfnqS+2XQzz//vMg1FERh3/s+ffq4jFHOvbYudy0XVO4Qq08++STfISSXM2DAAJf5u+66y6W2xMRElStXTv369XO2ubm55Tn/l9O/f3+X4RR33XWXsrOztX///kLXnKtv374u+4yMjJRlWerbt6+zzd3dXY0aNcr3fHfu3Nml57VJkyaKjIx0uYZ69eqlw4cPuwxtWrhwoXx8fNS1a9ci147Si3CLMqlx48b6+OOPdeLECW3dulXx8fE6efKkunXr5vKH6ZtvvlF0dLQqVKggf39/BQUFOcexXRhuzw+wkuTn5ydJql69er7tJ06ccGl3c3NTzZo1XdpuuukmSXK5v+j5jh49qtTUVM2dO1dBQUEuU58+fST9/y/JjRo1ShUrVlSTJk104403avDgwQX62PJycj8iz5X7h6p27dou7Z6enqpZs2aeP2TVqlXLM37Qz8+vwOctPwsWLNCtt94qb29vBQYGKigoSCtWrMjznl3Mhe+lJFWuXNnltffv35/nGCWpbt26zuXnu/A85fL29s7z8befn99Fz8v5NeTk5Gjq1Km68cYb5eXlpeuuu05BQUHauXNngY/1fPv379eNN96Ypz2/47xQcVxfFztHF3NhrbVq1ZKbm9tFf16KS2Hf+wuvp8qVK0sq2LVcEA8//LCaNWumxx9/XMHBwerevbs+/PDDAgXd/K6//K71qlWrqnz58i7r3XDDDYWq82qch8L83s3vdfK73m+66SaXa+iee+5R1apVtXDhQkl//9y999576tSp0yWHSKHsItyiTPP09FTjxo01adIkzZ49W2fPntVHH30k6e8vnrVp00bHjh3Tq6++qhUrVmjVqlUaMWKEJOX5w3F+72VB2q0CfFHscnJrePTRR7Vq1ap8p2bNmkn6+w/v7t279f7776t58+ZasmSJmjdvrnHjxl1RDQUZj3kpxX3e3n33XfXu3Vu1atXSW2+9pcTERK1atUqtW7cucK/W1XjPLnaeruT4J02apLi4ON1999169913tXLlSq1atUr169cvUg/elSiO6+tKr6UL/zNw4XyuC78sdbUV9XoqaP0+Pj7asGGDvvrqK/Xs2VM7d+7Uww8/rHvuueeyx3qx2q6GK/m5uthxFObnp6g/v+7u7nrkkUe0ZMkSnTlzRmvXrtXhw4f16KOPFml/KP34QhmM0ahRI0nS//73P0nS8uXLlZmZqU8//dSld6Cg37ovrJycHO3du9fZWytJv/76qyRd9Is1QUFBqlSpkrKzsxUdHX3Z16hQoYIefvhhPfzww8rKytIDDzygiRMnKj4+vthuExQWFibp73uDnt8TnZWVpaSkpALVeSUWL16smjVr6uOPP3YJB1ca4i8UFham3bt352n/5ZdfnMuvtsWLF6tVq1Z66623XNpTU1OdX4oqjLCwMP3222952vM7zvxc7vq6WFgrqt9++82lt3fPnj3Kyclx/rzk9gxeeJeO/D4GL0xtJfXeV65cOd+HSuRXv5ubm9q0aaM2bdro1Vdf1aRJk/TPf/5Ta9euveKfubCwMK1du1anT5926b0tyB0fCiu/Y87KynL+Xi5u+V3vv/76a57fub169dIrr7yi5cuX64svvlBQUJBiYmKuSk2wHz23KHPWrl2b7//gc8dY5X7cmPs///PXTUtL07x5865aba+//rrz35Zl6fXXX1e5cuXUpk2bfNd3d3dX165dtWTJEpfb1+Q6/zY7f/75p8syT09P1atXT5Zl6ezZs8V0BFJ0dLQ8PT01ffp0l3P31ltvKS0tTe3bty+218pPfu/bli1btGnTpmJ9nfvuu09bt2512W9GRobmzp2r8PDwQo0dLSp3d/c81/JHH31U5FsT3Xfffdq8ebO2bt3qbDt69Kjz49hLKcj1lXsP5OJ6CljuLdNy5T5hsF27dpIkX19fXXfdddqwYYPLerNmzcqzr8LUVlLvfa1atZSWlqadO3c62/73v/9p6dKlLusdP348z7a5dwXJzMy84jpiYmJ09uxZvfnmm862nJycPOe/ONSqVSvP+zV37tyr1tu+bNkyl5+XrVu3asuWLc5rKNett96qW2+9Vf/3f/+nJUuWqHv37lf9wSKwD+8sypyhQ4fq9OnT6tKli+rUqaOsrCxt3LhRH3zwgcLDw51jVe+99155enqqQ4cOeuKJJ3Tq1Cm9+eabqlKlylXpRfD29lZiYqJiY2MVGRmpL774QitWrNCzzz570VsSSX/fZmjt2rWKjIxUv379VK9ePR0/flw7duzQV1995fzDd++99yokJETNmjVTcHCwfv75Z73++utq3759sY4bCwoKUnx8vMaPH6+2bduqY8eO2r17t2bNmqXGjRtf9Y/y7r//fn388cfq0qWL2rdvr6SkJM2ZM0f16tXTqVOniu11Ro8erffee0/t2rXTsGHDFBAQoAULFigpKUlLliwpkQc03H///ZowYYL69OmjO++8Uz/88IMWLlyYZ+x2QY0cOdL5eOUnn3zSeSuwsLAwl4CVn4JcXw0bNpQk/fOf/1T37t1Vrlw5dejQocgP/khKSlLHjh3Vtm1bbdq0Se+++64eeeQRNWjQwLnO448/rsmTJ+vxxx9Xo0aNtGHDBucnIucrTG0l9d53795do0aNUpcuXTRs2DCdPn1as2fP1k033eTyhcEJEyZow4YNat++vcLCwnTkyBHNmjVL1apVc/nSW1F17txZTZo00VNPPaU9e/aoTp06+vTTT52/W4qzR/7xxx/XgAED1LVrV91zzz36/vvvtXLlyiJ9ElEQN9xwg5o3b66BAwcqMzNT06ZNU2BgoEaOHJln3V69ejmf0siQBLMRblHmTJkyRR999JE+//xzzZ07V1lZWapRo4YGDRqk5557zvnN49q1a2vx4sV67rnn9PTTTyskJEQDBw5UUFCQHnvssWKvy93dXYmJiRo4cKCeeeYZVapUSePGjdPYsWMvuV1wcLC2bt2qCRMm6OOPP9asWbMUGBio+vXr68UXX3Su98QTT2jhwoV69dVXderUKVWrVk3Dhg3Tc889V+zH8vzzzysoKEivv/66RowYoYCAAPXv31+TJk266o//7N27t5KTk/XGG29o5cqVqlevnt5991199NFHzodqFIfg4GBt3LhRo0aN0owZM3TmzBndeuutWr58+VXvnc717LPPKiMjQ4sWLdIHH3ygO+64QytWrNDo0aOLtL+qVatq7dq1Gjp0qCZPnqzAwEANGDBAoaGhLt8+z09Brq/GjRvrX//6l+bMmaPExETl5OQoKSmpyOH2gw8+0NixYzV69Gh5eHhoyJAhevnll13WGTt2rI4eParFixfrww8/VLt27fTFF1/keWBLYWorqfc+MDBQS5cuVVxcnEaOHKmIiAglJCTot99+cwm3HTt21L59+/T222/r2LFjuu6669SiRQuNHz/e+QWrK+Hu7q4VK1boySef1IIFC+Tm5qYuXbpo3LhxatasWbE++axfv35KSkpyjpe/6667tGrVqot+enWlevXqJTc3N02bNk1HjhxRkyZN9Prrr6tq1ap51u3Ro4dGjRqlWrVqudxRBOZxWMXxrRjgGte7d28tXry4WHsWAeBqWrZsmbp06aKvv/7a+cXVsmLfvn2KiIjQyy+/7OyNvZxjx46patWqGjt2rMaMGXOVK4SdGHMLAIDh/vrrL5f57OxszZgxQ76+vrrjjjtsqqpkzZ8/X9nZ2erZs6fdpeAqY1gCAACGGzp0qP766y9FRUUpMzNTH3/8sTZu3KhJkyZd8S3cSrs1a9bop59+0sSJE9W5c+fLPhYaZR/hFgAAw7Vu3VqvvPKKPvvsM505c0Y33HCDZsyY4fJURVNNmDBBGzduVLNmzZx35IDZGHMLAAAAYzDmFgAAAMYg3AIAAMAYjLnV309qOXz4sCpVqlTsj5cEAADAlbMsSydPnlRoaOglH7ZCuJV0+PBhVa9e3e4yAAAAcBkHDx5UtWrVLrqccCs5Hy158OBB+fr62lwNAAAALpSenq7q1atf9pHzhFv9/+dq+/r6Em4BAABKscsNIeULZQAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMW8Pthg0b1KFDB4WGhsrhcGjZsmV51vn555/VsWNH+fn5qUKFCmrcuLEOHDjgXH7mzBkNHjxYgYGBqlixorp27aqUlJQSPAoAAACUFraG24yMDDVo0EAzZ87Md/nvv/+u5s2bq06dOlq3bp127typMWPGyNvb27nOiBEjtHz5cn300Udav369Dh8+rAceeKCkDgEAAACliMOyLMvuIqS/b8i7dOlSde7c2dnWvXt3lStXTv/+97/z3SYtLU1BQUFatGiRunXrJkn65ZdfVLduXW3atElNmzYt0Gunp6fLz89PaWlpPMQBAACgFCpoXiu1Y25zcnK0YsUK3XTTTYqJiVGVKlUUGRnpMnRh+/btOnv2rKKjo51tderUUY0aNbRp0yYbqgYAAICdSm24PXLkiE6dOqXJkyerbdu2+vLLL9WlSxc98MADWr9+vSQpOTlZnp6e8vf3d9k2ODhYycnJF913Zmam0tPTXSYAAACUfR52F3AxOTk5kqROnTppxIgRkqTbbrtNGzdu1Jw5c9SiRYsi7zshIUHjx48vljoBAABQepTantvrrrtOHh4eqlevnkt73bp1nXdLCAkJUVZWllJTU13WSUlJUUhIyEX3HR8fr7S0NOd08ODBYq8fAAAAJa/UhltPT081btxYu3fvdmn/9ddfFRYWJklq2LChypUrp9WrVzuX7969WwcOHFBUVNRF9+3l5SVfX1+XCQAAAGWfrcMSTp06pT179jjnk5KS9N133ykgIEA1atTQM888o4cfflh33323WrVqpcTERC1fvlzr1q2TJPn5+alv376Ki4tTQECAfH19NXToUEVFRRX4TgkAAAAwh623Alu3bp1atWqVpz02Nlbz58+XJL399ttKSEjQH3/8odq1a2v8+PHq1KmTc90zZ87oqaee0nvvvafMzEzFxMRo1qxZlxyWcCFuBQYAAFC6FTSvlZr73NqJcAsAAFC6lfn73AIAAACFRbgFAACAMUrtfW4B4FoUPnqF3SXka9/k9naXAAAFQs8tAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGMPD7gIAwA7ho1fYXQIA4Cqg5xYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAM7pYAALis0nh3iX2T29tdAoBSiJ5bAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMWwNtxs2bFCHDh0UGhoqh8OhZcuWXXTdAQMGyOFwaNq0aS7tx48fV48ePeTr6yt/f3/17dtXp06durqFAwAAoFSyNdxmZGSoQYMGmjlz5iXXW7p0qTZv3qzQ0NA8y3r06KFdu3Zp1apV+uyzz7Rhwwb179//apUMAACAUszW+9y2a9dO7dq1u+Q6hw4d0tChQ7Vy5Uq1b+96T8Off/5ZiYmJ2rZtmxo1aiRJmjFjhu677z5NmTIl3zAMAAAAc5XqMbc5OTnq2bOnnnnmGdWvXz/P8k2bNsnf398ZbCUpOjpabm5u2rJlS0mWCgAAgFKgVD+h7MUXX5SHh4eGDRuW7/Lk5GRVqVLFpc3Dw0MBAQFKTk6+6H4zMzOVmZnpnE9PTy+eggEAAGCrUttzu337dr322muaP3++HA5Hse47ISFBfn5+zql69erFun8AAADYo9SG2//85z86cuSIatSoIQ8PD3l4eGj//v166qmnFB4eLkkKCQnRkSNHXLY7d+6cjh8/rpCQkIvuOz4+Xmlpac7p4MGDV/NQAAAAUEJK7bCEnj17Kjo62qUtJiZGPXv2VJ8+fSRJUVFRSk1N1fbt29WwYUNJ0po1a5STk6PIyMiL7tvLy0teXl5Xr3gAAADYwtZwe+rUKe3Zs8c5n5SUpO+++04BAQGqUaOGAgMDXdYvV66cQkJCVLt2bUlS3bp11bZtW/Xr109z5szR2bNnNWTIEHXv3p07JQAAAFyDbB2W8O233+r222/X7bffLkmKi4vT7bffrrFjxxZ4HwsXLlSdOnXUpk0b3XfffWrevLnmzp17tUoGAABAKWZrz23Lli1lWVaB19+3b1+etoCAAC1atKgYqwIAAEBZVWq/UAYAAAAUFuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYHnYXAMB84aNX2F0CAOAaQc8tAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDFsDbcbNmxQhw4dFBoaKofDoWXLljmXnT17VqNGjdItt9yiChUqKDQ0VL169dLhw4dd9nH8+HH16NFDvr6+8vf3V9++fXXq1KkSPhIAAACUBraG24yMDDVo0EAzZ87Ms+z06dPasWOHxowZox07dujjjz/W7t271bFjR5f1evTooV27dmnVqlX67LPPtGHDBvXv37+kDgEAAACliMOyLMvuIiTJ4XBo6dKl6ty580XX2bZtm5o0aaL9+/erRo0a+vnnn1WvXj1t27ZNjRo1kiQlJibqvvvu0x9//KHQ0NACvXZ6err8/PyUlpYmX1/f4jgcAOcJH73C7hJgoH2T29tdAoASVNC8VqbG3KalpcnhcMjf31+StGnTJvn7+zuDrSRFR0fLzc1NW7Zsueh+MjMzlZ6e7jIBAACg7Csz4fbMmTMaNWqU/vGPfzjTenJysqpUqeKynoeHhwICApScnHzRfSUkJMjPz885Va9e/arWDgAAgJJRJsLt2bNn9dBDD8myLM2ePfuK9xcfH6+0tDTndPDgwWKoEgAAAHbzsLuAy8kNtvv379eaNWtcxliEhIToyJEjLuufO3dOx48fV0hIyEX36eXlJS8vr6tWMwAAAOxRqntuc4Ptb7/9pq+++kqBgYEuy6OiopSamqrt27c729asWaOcnBxFRkaWdLkAAACwma09t6dOndKePXuc80lJSfruu+8UEBCgqlWrqlu3btqxY4c+++wzZWdnO8fRBgQEyNPTU3Xr1lXbtm3Vr18/zZkzR2fPntWQIUPUvXv3At8pAQAAAOawNdx+++23atWqlXM+Li5OkhQbG6vnn39en376qSTptttuc9lu7dq1atmypSRp4cKFGjJkiNq0aSM3Nzd17dpV06dPL5H6AQAAULrYGm5btmypS91mtyC34A0ICNCiRYuKsywAAACUUaV6zC0AAABQGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADG8LC7AAAAiiJ89Aq7S8jXvsnt7S4BuKbRcwsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxbA23GzZsUIcOHRQaGiqHw6Fly5a5LLcsS2PHjlXVqlXl4+Oj6Oho/fbbby7rHD9+XD169JCvr6/8/f3Vt29fnTp1qgSPAgAAAKWFreE2IyNDDRo00MyZM/Nd/tJLL2n69OmaM2eOtmzZogoVKigmJkZnzpxxrtOjRw/t2rVLq1at0meffaYNGzaof//+JXUIAAAAKEU87Hzxdu3aqV27dvkusyxL06ZN03PPPadOnTpJkt555x0FBwdr2bJl6t69u37++WclJiZq27ZtatSokSRpxowZuu+++zRlyhSFhoaW2LEAAADAfqV2zG1SUpKSk5MVHR3tbPPz81NkZKQ2bdokSdq0aZP8/f2dwVaSoqOj5ebmpi1btlx035mZmUpPT3eZAAAAUPaV2nCbnJwsSQoODnZpDw4Odi5LTk5WlSpVXJZ7eHgoICDAuU5+EhIS5Ofn55yqV69ezNUDAADADqU23F5N8fHxSktLc04HDx60uyQAAAAUg1IbbkNCQiRJKSkpLu0pKSnOZSEhITpy5IjL8nPnzun48ePOdfLj5eUlX19flwkAAABlX6kNtxEREQoJCdHq1audbenp6dqyZYuioqIkSVFRUUpNTdX27dud66xZs0Y5OTmKjIws8ZoBAABgL1vvlnDq1Cnt2bPHOZ+UlKTvvvtOAQEBqlGjhoYPH64XXnhBN954oyIiIjRmzBiFhoaqc+fOkqS6deuqbdu26tevn+bMmaOzZ89qyJAh6t69O3dKAAAAuAbZGm6//fZbtWrVyjkfFxcnSYqNjdX8+fM1cuRIZWRkqH///kpNTVXz5s2VmJgob29v5zYLFy7UkCFD1KZNG7m5ualr166aPn16iR8LAAAA7OewLMuyuwi7paeny8/PT2lpaYy/Ba6C8NEr7C4BKDH7Jre3uwTASAXNa6V2zC0AAABQWIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjOFhdwEAik/46BV2lwAAgK3ouQUAAIAxCLcAAAAwBuEWAAAAxihSuG3durVSU1PztKenp6t169ZXWhMAAABQJEUKt+vWrVNWVlae9jNnzug///nPFRcFAAAAFEWh7pawc+dO579/+uknJScnO+ezs7OVmJio66+/vviqAwAAAAqhUOH2tttuk8PhkMPhyHf4gY+Pj2bMmFFsxQEAAACFUahwm5SUJMuyVLNmTW3dulVBQUHOZZ6enqpSpYrc3d2LvUgAAACgIAoVbsPCwiRJOTk5V6UYAAAA4EoU+Qllv/32m9auXasjR47kCbtjx4694sIAAACAwipSuH3zzTc1cOBAXXfddQoJCZHD4XAuczgchFsAAADYokjh9oUXXtDEiRM1atSo4q4HAAAAKLIi3ef2xIkTevDBB4u7FgAAAOCKFCncPvjgg/ryyy+LuxYAAADgihRpWMINN9ygMWPGaPPmzbrllltUrlw5l+XDhg0rluIAAACAwnBYlmUVdqOIiIiL79Dh0N69e6+oqJKWnp4uPz8/paWlydfX1+5ygCILH73C7hKAa96+ye3tLgEwUkHzWpF6bpOSkopcGAAAAHC1FGnMLQAAAFAaFann9rHHHrvk8rfffrtIxQAAAABXokjh9sSJEy7zZ8+e1Y8//qjU1FS1bt26WAoDAAAACqtI4Xbp0qV52nJycjRw4EDVqlXriosCAAAAiqLYxty6ubkpLi5OU6dOLa5dAgAAAIVSrF8o+/3333Xu3Lni3CUAAABQYEUalhAXF+cyb1mW/ve//2nFihWKjY0tlsIkKTs7W88//7zeffddJScnKzQ0VL1799Zzzz0nh8PhfO1x48bpzTffVGpqqpo1a6bZs2frxhtvLLY6AAAAUDYUKdz+97//dZl3c3NTUFCQXnnllcveSaEwXnzxRc2ePVsLFixQ/fr19e2336pPnz7y8/NzPgXtpZde0vTp07VgwQJFRERozJgxiomJ0U8//SRvb+9iqwUAAAClX5HC7dq1a4u7jnxt3LhRnTp1Uvv2fz/tJTw8XO+99562bt0q6e9e22nTpum5555Tp06dJEnvvPOOgoODtWzZMnXv3r1E6gQAAEDpcEVjbo8ePaqvv/5aX3/9tY4ePVpcNTndeeedWr16tX799VdJ0vfff6+vv/5a7dq1k/T3k9KSk5MVHR3t3MbPz0+RkZHatGnTRfebmZmp9PR0lwkAAABlX5F6bjMyMjR06FC98847ysnJkSS5u7urV69emjFjhsqXL18sxY0ePVrp6emqU6eO3N3dlZ2drYkTJ6pHjx6SpOTkZElScHCwy3bBwcHOZflJSEjQ+PHji6VGAAAAlB5F6rmNi4vT+vXrtXz5cqWmpio1NVWffPKJ1q9fr6eeeqrYivvwww+1cOFCLVq0SDt27NCCBQs0ZcoULViw4Ir2Gx8fr7S0NOd08ODBYqoYAAAAdipSz+2SJUu0ePFitWzZ0tl23333ycfHRw899JBmz55dLMU988wzGj16tHPs7C233KL9+/crISFBsbGxCgkJkSSlpKSoatWqzu1SUlJ02223XXS/Xl5e8vLyKpYaAQAAUHoUqef29OnTeYYCSFKVKlV0+vTpKy7q/Ndxc3Mt0d3d3TkUIiIiQiEhIVq9erVzeXp6urZs2aKoqKhiqwMAAABlQ5HCbVRUlMaNG6czZ8442/766y+NHz++WENlhw4dNHHiRK1YsUL79u3T0qVL9eqrr6pLly6SJIfDoeHDh+uFF17Qp59+qh9++EG9evVSaGioOnfuXGx1AAAAoGwo0rCEadOmqW3btqpWrZoaNGgg6e87GXh5eenLL78stuJmzJihMWPGaNCgQTpy5IhCQ0P1xBNPaOzYsc51Ro4cqYyMDPXv31+pqalq3ry5EhMTucctAADANchhWZZVlA1Pnz6thQsX6pdffpEk1a1bVz169JCPj0+xFlgS0tPT5efnp7S0NPn6+tpdDlBk4aNX2F0CcM3bN7m93SUARipoXitSz21CQoKCg4PVr18/l/a3335bR48e1ahRo4qyWwAAAOCKFGnM7RtvvKE6derkaa9fv77mzJlzxUUBAAAARVGkcJucnOxy661cQUFB+t///nfFRQEAAABFUaRwW716dX3zzTd52r/55huFhoZecVEAAABAURRpzG2/fv00fPhwnT17Vq1bt5YkrV69WiNHjizWJ5QBAAAAhVGkcPvMM8/ozz//1KBBg5SVlSVJ8vb21qhRoxQfH1+sBQIAAAAFVaRw63A49OKLL2rMmDH6+eef5ePjoxtvvJFH2gIAAMBWRQq3uSpWrKjGjRsXVy0AAADAFSnSF8oAAACA0ohwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGB52FwAAgEnCR6+wu4Q89k1ub3cJQImh5xYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAY5T6cHvo0CE9+uijCgwMlI+Pj2655RZ9++23zuWWZWns2LGqWrWqfHx8FB0drd9++83GigEAAGCXUh1uT5w4oWbNmqlcuXL64osv9NNPP+mVV15R5cqVneu89NJLmj59uubMmaMtW7aoQoUKiomJ0ZkzZ2ysHAAAAHbwsLuAS3nxxRdVvXp1zZs3z9kWERHh/LdlWZo2bZqee+45derUSZL0zjvvKDg4WMuWLVP37t1LvGYAAADYp1T33H766adq1KiRHnzwQVWpUkW333673nzzTefypKQkJScnKzo62tnm5+enyMhIbdq06aL7zczMVHp6ussEAACAsq9Uh9u9e/dq9uzZuvHGG7Vy5UoNHDhQw4YN04IFCyRJycnJkqTg4GCX7YKDg53L8pOQkCA/Pz/nVL169at3EAAAACgxpTrc5uTk6I477tCkSZN0++23q3///urXr5/mzJlzRfuNj49XWlqaczp48GAxVQwAAAA7lepwW7VqVdWrV8+lrW7dujpw4IAkKSQkRJKUkpLisk5KSopzWX68vLzk6+vrMgEAAKDsK9XhtlmzZtq9e7dL26+//qqwsDBJf3+5LCQkRKtXr3YuT09P15YtWxQVFVWitQIAAMB+pfpuCSNGjNCdd96pSZMm6aGHHtLWrVs1d+5czZ07V5LkcDg0fPhwvfDCC7rxxhsVERGhMWPGKDQ0VJ07d7a3eAAAAJS4Uh1uGzdurKVLlyo+Pl4TJkxQRESEpk2bph49ejjXGTlypDIyMtS/f3+lpqaqefPmSkxMlLe3t42VAwAAwA4Oy7Isu4uwW3p6uvz8/JSWlsb4W5Rp4aNX2F0CgFJo3+T2dpcAXLGC5rVSPeYWAAAAKAzCLQAAAIxBuAUAAIAxCLcAAAAwRqm+WwJQmvHlLQAASh96bgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjlKlwO3nyZDkcDg0fPtzZdubMGQ0ePFiBgYGqWLGiunbtqpSUFPuKBAAAgG3KTLjdtm2b3njjDd16660u7SNGjNDy5cv10Ucfaf369Tp8+LAeeOABm6oEAACAncpEuD116pR69OihN998U5UrV3a2p6Wl6a233tKrr76q1q1bq2HDhpo3b542btyozZs321gxAAAA7FAmwu3gwYPVvn17RUdHu7Rv375dZ8+edWmvU6eOatSooU2bNl10f5mZmUpPT3eZAAAAUPZ52F3A5bz//vvasWOHtm3blmdZcnKyPD095e/v79IeHBys5OTki+4zISFB48ePL+5SAQAAYLNS3XN78OBBPfnkk1q4cKG8vb2Lbb/x8fFKS0tzTgcPHiy2fQMAAMA+pTrcbt++XUeOHNEdd9whDw8PeXh4aP369Zo+fbo8PDwUHBysrKwspaamumyXkpKikJCQi+7Xy8tLvr6+LhMAAADKvlI9LKFNmzb64YcfXNr69OmjOnXqaNSoUapevbrKlSun1atXq2vXrpKk3bt368CBA4qKirKjZAAAANioVIfbSpUq6eabb3Zpq1ChggIDA53tffv2VVxcnAICAuTr66uhQ4cqKipKTZs2taNkAAAA2KhUh9uCmDp1qtzc3NS1a1dlZmYqJiZGs2bNsrssAAAA2MBhWZZldxF2S09Pl5+fn9LS0hh/iwILH73C7hIAoED2TW5vdwnAFStoXivVXygDAAAACoNwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjeNhdAAAAuLrCR6+wu4R87Zvc3u4SYCB6bgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDE87C4AAABcm8JHr7C7hDz2TW5vdwm4QvTcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDFKfbhNSEhQ48aNValSJVWpUkWdO3fW7t27XdY5c+aMBg8erMDAQFWsWFFdu3ZVSkqKTRUDAADALqU+3K5fv16DBw/W5s2btWrVKp09e1b33nuvMjIynOuMGDFCy5cv10cffaT169fr8OHDeuCBB2ysGgAAAHYo9Q9xSExMdJmfP3++qlSpou3bt+vuu+9WWlqa3nrrLS1atEitW7eWJM2bN09169bV5s2b1bRpUzvKBgAAgA1Kfc/thdLS0iRJAQEBkqTt27fr7Nmzio6Odq5Tp04d1ahRQ5s2bcp3H5mZmUpPT3eZAAAAUPaVqXCbk5Oj4cOHq1mzZrr55pslScnJyfL09JS/v7/LusHBwUpOTs53PwkJCfLz83NO1atXv9qlAwAAoASUqXA7ePBg/fjjj3r//fevaD/x8fFKS0tzTgcPHiymCgEAAGCnUj/mNteQIUP02WefacOGDapWrZqzPSQkRFlZWUpNTXXpvU1JSVFISEi++/Ly8pKXl9fVLhkAAAAlrNT33FqWpSFDhmjp0qVas2aNIiIiXJY3bNhQ5cqV0+rVq51tu3fv1oEDBxQVFVXS5QIAAMBGpb7ndvDgwVq0aJE++eQTVapUyTmO1s/PTz4+PvLz81Pfvn0VFxengIAA+fr6aujQoYqKiuJOCQAAANeYUh9uZ8+eLUlq2bKlS/u8efPUu3dvSdLUqVPl5uamrl27KjMzUzExMZo1a1YJVwoAAAC7OSzLsuwuwm7p6eny8/NTWlqafH197S4HFwgfvcLuEgAA14h9k9vbXQIuoqB5rdSPuQUAAAAKinALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGB52F4DSJXz0CrtLAAAAKDJ6bgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABjDw+4CrlXho1fYXQIAALhAaf37vG9ye7tLKDPouQUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDuyUAAACUctzFoeDouQUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjGBNuZ86cqfDwcHl7eysyMlJbt261uyQAAACUMCPC7QcffKC4uDiNGzdOO3bsUIMGDRQTE6MjR47YXRoAAABKkBHh9tVXX1W/fv3Up08f1atXT3PmzFH58uX19ttv210aAAAASlCZf0JZVlaWtm/frvj4eGebm5uboqOjtWnTpny3yczMVGZmpnM+LS1NkpSenn51iz1PTubpEnstAACAq6Eks1Pua1mWdcn1yny4PXbsmLKzsxUcHOzSHhwcrF9++SXfbRISEjR+/Pg87dWrV78qNQIAAJjIb1rJv+bJkyfl5+d30eVlPtwWRXx8vOLi4pzzOTk5On78uAIDA+VwOIr1tdLT01W9enUdPHhQvr6+xbpv5I9zXvI45yWPc17yOOf24LyXvNJ6zi3L0smTJxUaGnrJ9cp8uL3uuuvk7u6ulJQUl/aUlBSFhITku42Xl5e8vLxc2vz9/a9WiZIkX1/fUnWBXAs45yWPc17yOOclj3NuD857ySuN5/xSPba5yvwXyjw9PdWwYUOtXr3a2ZaTk6PVq1crKirKxsoAAABQ0sp8z60kxcXFKTY2Vo0aNVKTJk00bdo0ZWRkqE+fPnaXBgAAgBJkRLh9+OGHdfToUY0dO1bJycm67bbblJiYmOdLZnbw8vLSuHHj8gyDwNXDOS95nPOSxzkveZxze3DeS15ZP+cO63L3UwAAAADKiDI/5hYAAADIRbgFAACAMQi3AAAAMAbhFgAAAMYg3Jagjh07qkaNGvL29lbVqlXVs2dPHT582O6yjLZv3z717dtXERER8vHxUa1atTRu3DhlZWXZXZrRJk6cqDvvvFPly5e/6g9IuVbNnDlT4eHh8vb2VmRkpLZu3Wp3SUbbsGGDOnTooNDQUDkcDi1btszukoyWkJCgxo0bq1KlSqpSpYo6d+6s3bt3212W0WbPnq1bb73V+eCGqKgoffHFF3aXVSSE2xLUqlUrffjhh9q9e7eWLFmi33//Xd26dbO7LKP98ssvysnJ0RtvvKFdu3Zp6tSpmjNnjp599lm7SzNaVlaWHnzwQQ0cONDuUoz0wQcfKC4uTuPGjdOOHTvUoEEDxcTE6MiRI3aXZqyMjAw1aNBAM2fOtLuUa8L69es1ePBgbd68WatWrdLZs2d17733KiMjw+7SjFWtWjVNnjxZ27dv17fffqvWrVurU6dO2rVrl92lFRq3ArPRp59+qs6dOyszM1PlypWzu5xrxssvv6zZs2dr7969dpdivPnz52v48OFKTU21uxSjREZGqnHjxnr99dcl/f1UxurVq2vo0KEaPXq0zdWZz+FwaOnSpercubPdpVwzjh49qipVqmj9+vW6++677S7nmhEQEKCXX35Zffv2tbuUQqHn1ibHjx/XwoULdeeddxJsS1haWpoCAgLsLgMokqysLG3fvl3R0dHONjc3N0VHR2vTpk02VgZcPWlpaZLE7+4Skp2drffff18ZGRmKioqyu5xCI9yWsFGjRqlChQoKDAzUgQMH9Mknn9hd0jVlz549mjFjhp544gm7SwGK5NixY8rOzs7zBMbg4GAlJyfbVBVw9eTk5Gj48OFq1qyZbr75ZrvLMdoPP/ygihUrysvLSwMGDNDSpUtVr149u8sqNMLtFRo9erQcDsclp19++cW5/jPPPKP//ve/+vLLL+Xu7q5evXqJkSGFV9jzLkmHDh1S27Zt9eCDD6pfv342VV52FeWcA8CVGjx4sH788Ue9//77dpdivNq1a+u7777Tli1bNHDgQMXGxuqnn36yu6xCY8ztFTp69Kj+/PPPS65Ts2ZNeXp65mn/448/VL16dW3cuLFMdvvbqbDn/fDhw2rZsqWaNm2q+fPny82N/9cVVlGudcbcFr+srCyVL19eixcvdhnzGRsbq9TUVD4NKgGMuS05Q4YM0SeffKINGzYoIiLC7nKuOdHR0apVq5beeOMNu0spFA+7CyjrgoKCFBQUVKRtc3JyJEmZmZnFWdI1oTDn/dChQ2rVqpUaNmyoefPmEWyL6EqudRQfT09PNWzYUKtXr3aGq5ycHK1evVpDhgyxtzigmFiWpaFDh2rp0qVat24dwdYmOTk5ZTKjEG5LyJYtW7Rt2zY1b95clStX1u+//64xY8aoVq1a9NpeRYcOHVLLli0VFhamKVOm6OjRo85lISEhNlZmtgMHDuj48eM6cOCAsrOz9d1330mSbrjhBlWsWNHe4gwQFxen2NhYNWrUSE2aNNG0adOUkZGhPn362F2asU6dOqU9e/Y455OSkvTdd98pICBANWrUsLEyMw0ePFiLFi3SJ598okqVKjnHk/v5+cnHx8fm6swUHx+vdu3aqUaNGjp58qQWLVqkdevWaeXKlXaXVngWSsTOnTutVq1aWQEBAZaXl5cVHh5uDRgwwPrjjz/sLs1o8+bNsyTlO+HqiY2Nzfecr1271u7SjDFjxgyrRo0alqenp9WkSRNr8+bNdpdktLVr1+Z7TcfGxtpdmpEu9nt73rx5dpdmrMcee8wKCwuzPD09raCgIKtNmzbWl19+aXdZRcKYWwAAABiDwYcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGOP/AUPbzNadgBFpAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "a = np.random.normal(loc=0, scale=1, size=1_000)\n", "\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "ax.hist(a, color=\"C0\", bins=15)\n", "ax.set(title=\"Samples from a normal distribution using numpy\", ylabel=\"count\");" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now let's try to do it in PyTensor." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "TensorType(float64, ())" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = pt.random.normal(loc=0, scale=1, name=\"y\")\n", "y.type" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Next, we show the graph using {func}`~pytensor.dprint`." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normal_rv{0, (0, 0), floatX, False}.1 [id A] 'y'\n", " |RandomGeneratorSharedVariable() [id B]\n", " |TensorConstant{[]} [id C]\n", " |TensorConstant{11} [id D]\n", " |TensorConstant{0} [id E]\n", " |TensorConstant{1} [id F]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pytensor.dprint(y)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The inputs are always in the following order:\n", "1. `rng` shared variable\n", "2. `size`\n", "3. `dtype` (number code)\n", "4. `arg1`\n", "5. `arg2`\n", "6. `argn`" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We *could* sample by calling {meth}`~pytensor.graph.basic.Variable.eval`. on the random variable." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(-1.4186441)" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.eval()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Note however that these samples are always the same!" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sample 0: -1.4186441029543038\n", "Sample 1: -1.4186441029543038\n", "Sample 2: -1.4186441029543038\n", "Sample 3: -1.4186441029543038\n", "Sample 4: -1.4186441029543038\n", "Sample 5: -1.4186441029543038\n", "Sample 6: -1.4186441029543038\n", "Sample 7: -1.4186441029543038\n", "Sample 8: -1.4186441029543038\n", "Sample 9: -1.4186441029543038\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {y.eval()}\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We always get the same samples! This has to do with the random seed step in the graph, i.e. `RandomGeneratorSharedVariable` (we will not go deeper into this subject here). We will show how to generate different samples with `pymc` below." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## PyMC\n", "\n", "![pymc logo](https://raw.githubusercontent.com/pymc-devs/pymc/main/docs/logos/PyMC.png)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "To do so, we start by defining a `pymc` normal distribution." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normal_rv{0, (0, 0), floatX, False}.1 [id A]\n", " |RandomGeneratorSharedVariable() [id B]\n", " |TensorConstant{[]} [id C]\n", " |TensorConstant{11} [id D]\n", " |TensorConstant{0} [id E]\n", " |TensorConstant{1.0} [id F]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = pm.Normal.dist(mu=0, sigma=1)\n", "pytensor.dprint(x)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Observe that `x` is just a normal `RandomVariable` and which is the same as `y` except for the `rng`." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can try to generate samples by calling {meth}`~pytensor.graph.basic.Variable.eval` as above." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sample 0: 1.3064743941879295\n", "Sample 1: 1.3064743941879295\n", "Sample 2: 1.3064743941879295\n", "Sample 3: 1.3064743941879295\n", "Sample 4: 1.3064743941879295\n", "Sample 5: 1.3064743941879295\n", "Sample 6: 1.3064743941879295\n", "Sample 7: 1.3064743941879295\n", "Sample 8: 1.3064743941879295\n", "Sample 9: 1.3064743941879295\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {x.eval()}\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "As before we get the same value for all iterations. The correct way to generate random samples is using {func}`~pymc.draw`." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAIQCAYAAACbhEYhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCc0lEQVR4nO3deVyU9f7//+egMiCrKGsiopamqZUmUuaeu2appZnhkpaipdTJ6Jzczik82mKZS/V16XzSY1kuZSfLBbUSzTQzLT0ulJbiLigmKFy/P/wxxxGURWDg3eN+u83txvW+3nPN67rmYnjynvdcY7MsyxIAAABgADdXFwAAAAAUF8ItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi1QRthsNk2YMKFUH/Po0aPq3bu3qlatKpvNpmnTppXq4+PG/PLLL7LZbJo/f36R7j9//nzZbDb98ssvjrbWrVurdevWxVJffq4+5ydMmCCbzaYTJ06UyuPXrFlTAwcOLJXHKik3eg4AJiLcwig//vijevfurYiICHl4eOimm27Sfffdp+nTp7u6tDJpzJgx+uKLLxQfH6//+7//U6dOnVxdEsqhjRs3asKECTpz5oyrS8mlLNcGoGRUdHUBQHHZuHGj2rRpoxo1amjo0KEKCQnRoUOHtGnTJr3xxhsaNWqUq0ssc9auXav7779fzz77rKtLQRnx5ZdfFvo+Gzdu1MSJEzVw4ED5+/sX+H5//PGHKlYs2T9D16ttz549cnMr32M8ERER+uOPP1SpUiVXlwKUGYRbGOOll16Sn5+ftmzZkuuP2LFjx1xTVBl37NixAoWR9PR0eXl5lXxB5dyFCxfk7u5ergOTu7t7iW4/OztbmZmZ8vDwkIeHR4k+Vn7sdrtLH7842Gw2lx9HoKwpv6/AwFX279+vBg0a5BnWgoKCnJbnzZuntm3bKigoSHa7XfXr19esWbNy3a9mzZrq1q2b1q1bp6ZNm8rT01MNGzbUunXrJElLlixRw4YN5eHhoSZNmuj77793uv/AgQPl7e2tAwcOqGPHjvLy8lJYWJgmTZoky7Ly3afff/9dgwcPVnBwsOx2uxo0aKC5c+fm6jd9+nQ1aNBAlStXVpUqVdS0aVMtXLjwmtvNmWtpWZZmzJghm80mm83mtG79+vUaMWKEgoKCVL16dcd9Z86cqQYNGshutyssLEyxsbG53vJt3bq1brvtNu3YsUOtWrVS5cqVVadOHX300UeSpPXr1ysqKkqenp6qW7euVq9ene+xyMzM1Lhx49SkSRP5+fnJy8tL9957rxITE/O9r/S/5/Lrr79Ws2bN5OHhoVq1aulf//pXrr4HDhxQnz59FBAQoMqVK6t58+b67LPPnPqsW7dONptNixYt0t/+9jfddNNNqly5stLS0hzP+8GDB9WtWzd5e3vrpptu0owZMyRdnj7Ttm1beXl5KSIiItdzderUKT377LNq2LChvL295evrq86dO+uHH34o0L7mZdeuXWrbtq08PT1VvXp1/eMf/1B2dnaufnnNub3e+TVhwgT95S9/kSRFRkY6zqWcebw2m00jR47UggULHOfNypUrHevymmd+4sQJPfTQQ/L19VXVqlX19NNP68KFC47115tneuU286strzm3hXnuP/zwQ7300kuqXr26PDw81K5dO+3bty9XTVcbOHCgatasmas9Z87xlVatWqUWLVrI399f3t7eqlu3rl544YXrHouc8+/3339Xz5495e3trcDAQD377LPKyspy2v7Jkyc1YMAA+fr6yt/fXzExMfrhhx8KNI8357Viw4YNeuKJJ1S1alX5+vrqscce0+nTpx39YmJiVK1aNV28eDHXNjp06KC6des6lnPOl8WLF6t+/fry9PRUdHS0fvzxR0nS22+/rTp16sjDw0OtW7d2mi+eY/PmzerSpYuqVKkiLy8vNWrUSG+88cZ19wVmYeQWxoiIiFBSUpJ27typ22677bp9Z82apQYNGqhHjx6qWLGiPv30U40YMULZ2dmKjY116rtv3z498sgjeuKJJ/Too4/qlVdeUffu3TV79my98MILGjFihCQpISFBDz30UK63OrOystSpUyc1b95cU6ZM0cqVKzV+/HhdunRJkyZNumaNR48eVfPmzR0v9oGBgfr88881ZMgQpaWlafTo0ZKkd999V0899ZR69+7tCAE7duzQ5s2b9cgjj+S57ZYtW+r//u//NGDAAN1333167LHHcvUZMWKEAgMDNW7cOKWnp0u6/Md34sSJat++vYYPH649e/Zo1qxZ2rJli7755hunt0ZPnz6tbt26qW/fvurTp49mzZqlvn37asGCBRo9erSefPJJPfLII5o6dap69+6tQ4cOycfH55rHIy0tTf/v//0/9evXT0OHDtXZs2c1Z84cdezYUd9++61uv/32a943x759+9S7d28NGTJEMTExmjt3rgYOHKgmTZqoQYMGjuN+99136/z583rqqadUtWpVvffee+rRo4c++ugjPfDAA07b/Pvf/y53d3c9++yzysjIcIx8ZmVlqXPnzmrZsqWmTJmiBQsWaOTIkfLy8tJf//pX9e/fXw8++KBmz56txx57TNHR0YqMjJR0OWAtW7ZMffr0UWRkpI4ePaq3335brVq10k8//aSwsLB89/VKKSkpatOmjS5duqTnn39eXl5eeuedd+Tp6ZnvffM7vx588EH997//1b///W+9/vrrqlatmiQpMDDQsY21a9fqww8/1MiRI1WtWrU8g92VHnroIdWsWVMJCQnatGmT3nzzTZ0+fTrPf0SupyC1Xamwz/3kyZPl5uamZ599VqmpqZoyZYr69++vzZs3F6rOa9m1a5e6deumRo0aadKkSbLb7dq3b5+++eabfO+blZWljh07KioqSq+88opWr16tV199VbVr19bw4cMlXR5F7969u7799lsNHz5c9erV0/LlyxUTE1OoOkeOHCl/f39NmDDB8Zrw66+/Ov4JGDBggP71r3/piy++ULdu3Rz3S0lJ0dq1azV+/Hin7X311Vf65JNPHK/FCQkJ6tatm5577jnNnDlTI0aM0OnTpzVlyhQNHjxYa9euddx31apV6tatm0JDQ/X0008rJCREP//8s1asWKGnn366UPuFcswCDPHll19aFSpUsCpUqGBFR0dbzz33nPXFF19YmZmZufqeP38+V1vHjh2tWrVqObVFRERYkqyNGzc62r744gtLkuXp6Wn9+uuvjva3337bkmQlJiY62mJiYixJ1qhRoxxt2dnZVteuXS13d3fr+PHjjnZJ1vjx4x3LQ4YMsUJDQ60TJ0441dS3b1/Lz8/PsQ/333+/1aBBg3yOTt4kWbGxsU5t8+bNsyRZLVq0sC5duuRoP3bsmOXu7m516NDBysrKcrS/9dZbliRr7ty5jrZWrVpZkqyFCxc62nbv3m1Jstzc3KxNmzY52nOO57x5865b66VLl6yMjAynttOnT1vBwcHW4MGD893XnOdyw4YNTvtkt9utZ555xtE2evRoS5L11VdfOdrOnj1rRUZGWjVr1nTse2JioiXJqlWrVq7zKed5f/nll51q9fT0tGw2m7Vo0aJcx+XK5/7ChQtOx9iyLCs5Odmy2+3WpEmTnNoKcuxy9mnz5s1O++7n52dJspKTkx3trVq1slq1auVYLsj5NXXq1FzbyZHznO/atSvPdVfu9/jx4y1JVo8ePZz6jRgxwpJk/fDDD5ZlXX+/r97m9WqLiIiwYmJiHMuFfe5vvfVWp3PyjTfesCRZP/74Y67HulJMTIwVERGRqz1n/3O8/vrrliSn14mr5XUscs6/K88Vy7KsO+64w2rSpIlj+eOPP7YkWdOmTXO0ZWVlWW3bti3QeZXzWtGkSROn19kpU6ZYkqzly5c7tlm9enXr4Ycfdrr/a6+9ZtlsNuvAgQOONkmW3W53er5yXltDQkKstLQ0R3t8fLzTc3vp0iUrMjLSioiIsE6fPu30WNnZ2dfdF5iFaQkwxn333aekpCT16NFDP/zwg6ZMmaKOHTvqpptu0ieffOLU98oRq9TUVJ04cUKtWrXSgQMHlJqa6tS3fv36io6OdixHRUVJktq2basaNWrkaj9w4ECu2kaOHOn4OWckNjMz85pvx1uWpY8//ljdu3eXZVk6ceKE49axY0elpqZq27ZtkiR/f3/99ttv2rJlS4GOU0ENHTpUFSpUcCyvXr1amZmZGj16tNPI9NChQ+Xr65vrrVtvb2/17dvXsVy3bl35+/vr1ltvdRwr6frH7UoVKlRwjIpmZ2fr1KlTunTpkpo2beo4FvmpX7++7r33XsdyYGCg6tat6/TY//nPf9SsWTO1aNHCaV+GDRumX375RT/99JPTNmNiYq45Avr44487fvb391fdunXl5eWlhx56yNGec1yurMFutzuOcVZWlk6ePOl4S7qg+3ql//znP2revLmaNWvmtO/9+/fP977FcX61atVK9evXL3D/q989yfkw6H/+858i11AQhX3uBw0a5DRHOefcyu9cLqicKVbLly/PcwpJfp588kmn5XvvvdeptpUrV6pSpUoaOnSoo83NzS3X8c/PsGHDnN61GT58uCpWrOh4vtzc3NS/f3998sknOnv2rKPfggULdPfddzvescjRrl07p9H9nNeIXr16Ob27c/Vrx/fff6/k5GSNHj061/S0q6d7wGyEWxjlrrvu0pIlS3T69Gl9++23io+P19mzZ9W7d2+nP0zffPON2rdvLy8vL/n7+yswMNAxj+3qcHtlgJUkPz8/SVJ4eHie7VfONZMuv7DXqlXLqe2WW26RpDzni0nS8ePHdebMGb3zzjsKDAx0ug0aNEjS/z4kN3bsWHl7e6tZs2a6+eabFRsbW6C3LfNz9R+cX3/9VZKc5sdJlz+AVKtWLcf6HNWrV8/1B8XPz6/Axy0v7733nho1aiQPDw9VrVpVgYGB+uyzz3I9Z9dy9XMpSVWqVHF67F9//TXXPkrSrbfe6lh/pauPUw4PD49cb3/7+fld87hcWUN2drZef/113XzzzbLb7apWrZoCAwO1Y8eOAu/rlX799VfdfPPNudrz2s+rFcf5da1jdC1X11q7dm25ubld8/eluBT2ub/6fKpSpYqkgp3LBfHwww/rnnvu0eOPP67g4GD17dtXH374YYGCbl7nX17nemhoqCpXruzUr06dOoWq8+rny9vbW6GhoU7P12OPPaY//vhDS5culXT5ShVbt27VgAEDcm2vqK+5+/fvl6R8p6XBfIRbGMnd3V133XWXXn75Zc2aNUsXL17U4sWLJV1+AWzXrp1OnDih1157TZ999plWrVqlMWPGSFKuPxxXjl4WpN0qwAfF8pNTw6OPPqpVq1blebvnnnskXf7Du2fPHi1atEgtWrTQxx9/rBYtWuSax1ZYBZmPeT3Ffdzef/99DRw4ULVr19acOXO0cuVKrVq1Sm3bti3wqFZJPGfXOk43sv8vv/yy4uLi1LJlS73//vv64osvtGrVKjVo0KBII3g3ojjOrxs9l67+Z+Bao3BXf1iqpBX1fCpo/Z6entqwYYNWr16tAQMGaMeOHXr44Yd133335buv16rNVerXr68mTZro/fffl3T599nd3d3pXYwcrnjNhVn4QBmM17RpU0nSkSNHJEmffvqpMjIy9MknnziNEBT0U/eFlZ2drQMHDjhGayXpv//9ryRd84M1gYGB8vHxUVZWltq3b5/vY3h5eenhhx/Www8/rMzMTD344IN66aWXFB8fX2yXCYqIiJB0ecTlypHozMxMJScnF6jOG/HRRx+pVq1aWrJkiVM4uNEQf7WIiAjt2bMnV/vu3bsd60vaRx99pDZt2mjOnDlO7WfOnHF8KKowIiIitHfv3lztee1nXvI7v4r7Ld+9e/c6jfbu27dP2dnZjt+XnBHSq6/ScfXIqlS4t6NL67mvUqVKnl8qkVf9bm5uateundq1a6fXXntNL7/8sv76178qMTHxhn/nIiIilJiYqPPnzzuN3hbkig9X2rt3r9q0aeNYPnfunI4cOaIuXbo49XvssccUFxenI0eOaOHCheratavjuSwOtWvXliTt3LmzxF+PULYxcgtjJCYm5vkffM68r5y3G3P++7+yb2pqqubNm1ditb311luOny3L0ltvvaVKlSqpXbt2efavUKGCevXqpY8//lg7d+7Mtf748eOOn0+ePOm0zt3dXfXr15dlWXleeqeo2rdvL3d3d7355ptOx27OnDlKTU1V165di+2x8pLX87Z582YlJSUV6+N06dJF3377rdN209PT9c4776hmzZqFmjtaVBUqVMh1Li9evFi///57kbbXpUsXbdq0Sd9++62j7fjx41qwYEG+9y3I+ZVzDeTi+hawnEum5cj5hsHOnTtLknx9fVWtWjVt2LDBqd/MmTNzbaswtZXWc1+7dm2lpqZqx44djrYjR4443rLPcerUqVz3zbkqSEZGxg3X0bFjR128eFHvvvuuoy07OzvX8c/PO++84/RaM2vWLF26dMnxfOXo16+fbDabnn76aR04cECPPvroje3AVe68805FRkZq2rRpuZ5vRnf/XBi5hTFGjRql8+fP64EHHlC9evWUmZmpjRs36oMPPlDNmjUdc1U7dOggd3d3de/eXU888YTOnTund999V0FBQY7R3eLk4eGhlStXKiYmRlFRUfr888/12Wef6YUXXrjmJYmky5cZSkxMVFRUlIYOHar69evr1KlT2rZtm1avXu34w9ehQweFhITonnvuUXBwsH7++We99dZb6tq163UvrVVYgYGBio+P18SJE9WpUyf16NFDe/bs0cyZM3XXXXcV+x+qq3Xr1k1LlizRAw88oK5duyo5OVmzZ89W/fr1de7cuWJ7nOeff17//ve/1blzZz311FMKCAjQe++9p+TkZH388cel8gUN3bp106RJkzRo0CDdfffd+vHHH7VgwYJcc7cL6rnnnnN8vfLTTz/tuBRYRESEU8DKS0HOryZNmkiS/vrXv6pv376qVKmSunfvXuQv/khOTlaPHj3UqVMnJSUl6f3339cjjzyixo0bO/o8/vjjmjx5sh5//HE1bdpUGzZscLwjcqXC1FZaz33fvn01duxYPfDAA3rqqad0/vx5zZo1S7fccovTBwYnTZqkDRs2qGvXroqIiNCxY8c0c+ZMVa9e3elDb0XVs2dPNWvWTM8884z27dunevXq6ZNPPnG8thR01DszM1Pt2rVzXApx5syZatGihXr06OHULzAwUJ06ddLixYvl7+9f7P8Qu7m5adasWerevbtuv/12DRo0SKGhodq9e7d27dqlL774olgfD2VYqV+fASghn3/+uTV48GCrXr16lre3t+Xu7m7VqVPHGjVqlHX06FGnvp988onVqFEjy8PDw6pZs6b1z3/+05o7d26uSwZFRERYXbt2zfVYyuMSWjmX5Jk6daqjLSYmxvLy8rL2799vdejQwapcubIVHBxsjR8/PtelnnTVJYwsy7KOHj1qxcbGWuHh4ValSpWskJAQq127dtY777zj6PP2229bLVu2tKpWrWrZ7Xardu3a1l/+8hcrNTU132OW137kXN5ny5Yted7nrbfesurVq2dVqlTJCg4OtoYPH57rsjutWrXK8/JRhTmeV8vOzrZefvllKyIiwrLb7dYdd9xhrVix4pqXVSroY1996SvLsqz9+/dbvXv3tvz9/S0PDw+rWbNm1ooVK5z65FwOavHixbm2mfO85/VYBTkuFy5csJ555hkrNDTU8vT0tO655x4rKSkpV60FvRSYZVnWjh07rFatWlkeHh7WTTfdZP3973+35syZk++lwAp6fv3973+3brrpJsvNzc1pm9d7bq8+53MuhfXTTz9ZvXv3tnx8fKwqVapYI0eOtP744w+n+54/f94aMmSI5efnZ/n4+FgPPfSQdezYsTx/j65V29WXArOsG3vuC/N8fPnll9Ztt91mubu7W3Xr1rXef//9XJcCW7NmjXX//fdbYWFhlru7uxUWFmb169fP+u9//3vdx7zW+Xf19i3Lso4fP2498sgjlo+Pj+Xn52cNHDjQ+uabbyxJTpesy0vOa8X69eutYcOGWVWqVLG8vb2t/v37WydPnszzPh9++KElyRo2bFie6wv62mpZ134evv76a+u+++6zfHx8LC8vL6tRo0bW9OnTr7svMIvNshirB0rKwIED9dFHHxXryCIAlKRly5bpgQce0Ndff+344Gpe5s+fr0GDBmnLli2OzzbkZ/ny5erZs6c2bNjgdFk+oDgx5xYAgD+pP/74w2k5KytL06dPl6+vr+68885if7x3331XtWrVKpZpFcC1MOcWAIA/qVGjRumPP/5QdHS0MjIytGTJEm3cuFEvv/zyDV/C7UqLFi3Sjh079Nlnn+mNN97gSxVQogi3AAD8SbVt21avvvqqVqxYoQsXLqhOnTqaPn2607cqFod+/frJ29tbQ4YM0YgRI4p128DVmHMLAAAAYzDnFgAAAMYg3AIAAMAYzLnV5W9kOXz4sHx8fJjkDgAAUAZZlqWzZ88qLCzsul+qQriVdPjwYYWHh7u6DAAAAOTj0KFDql69+jXXE24lx1dIHjp0SL6+vi6uBgAAAFdLS0tTeHh4vl8tT7jV/74/29fXl3ALAABQhuU3hZQPlAEAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGNUdOWDJyQkaMmSJdq9e7c8PT11991365///Kfq1q3r6HPhwgU988wzWrRokTIyMtSxY0fNnDlTwcHBjj4HDx7U8OHDlZiYKG9vb8XExCghIUEVK7p09wCg8Cb4ubqC3CakuroCACgwl47crl+/XrGxsdq0aZNWrVqlixcvqkOHDkpPT3f0GTNmjD799FMtXrxY69ev1+HDh/Xggw861mdlZalr167KzMzUxo0b9d5772n+/PkaN26cK3YJAAAALmSzLMtydRE5jh8/rqCgIK1fv14tW7ZUamqqAgMDtXDhQvXu3VuStHv3bt16661KSkpS8+bN9fnnn6tbt246fPiwYzR39uzZGjt2rI4fPy53d/d8HzctLU1+fn5KTU2Vr69vie4jAFwXI7cAkKeC5rUyNec2NfXyC2hAQIAkaevWrbp48aLat2/v6FOvXj3VqFFDSUlJkqSkpCQ1bNjQaZpCx44dlZaWpl27dpVi9QAAAHC1MjMpNTs7W6NHj9Y999yj2267TZKUkpIid3d3+fv7O/UNDg5WSkqKo8+VwTZnfc66vGRkZCgjI8OxnJaWVly7AQAAABcqMyO3sbGx2rlzpxYtWlTij5WQkCA/Pz/HLTw8vMQfEwAAACWvTITbkSNHasWKFUpMTFT16tUd7SEhIcrMzNSZM2ec+h89elQhISGOPkePHs21PmddXuLj45Wamuq4HTp0qBj3BgAAAK7i0nBrWZZGjhyppUuXau3atYqMjHRa36RJE1WqVElr1qxxtO3Zs0cHDx5UdHS0JCk6Olo//vijjh075uizatUq+fr6qn79+nk+rt1ul6+vr9MNAAAA5Z9L59zGxsZq4cKFWr58uXx8fBxzZP38/OTp6Sk/Pz8NGTJEcXFxCggIkK+vr0aNGqXo6Gg1b95cktShQwfVr19fAwYM0JQpU5SSkqK//e1vio2Nld1ud+XuAQAAoJS5NNzOmjVLktS6dWun9nnz5mngwIGSpNdff11ubm7q1auX05c45KhQoYJWrFih4cOHKzo6Wl5eXoqJidGkSZNKazcAAABQRpSp69y6Cte5BVBmcJ1bAMhTubzOLQAAAHAjCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMEZFVxcAAC4xwc/VFQAASgAjtwAAADCGS8Pthg0b1L17d4WFhclms2nZsmVO6202W563qVOnOvrUrFkz1/rJkyeX8p4AAACgLHBpuE1PT1fjxo01Y8aMPNcfOXLE6TZ37lzZbDb16tXLqd+kSZOc+o0aNao0ygcAAEAZ49I5t507d1bnzp2vuT4kJMRpefny5WrTpo1q1arl1O7j45OrLwAAAP58ys2c26NHj+qzzz7TkCFDcq2bPHmyqlatqjvuuENTp07VpUuXrrutjIwMpaWlOd0AAABQ/pWbqyW899578vHx0YMPPujU/tRTT+nOO+9UQECANm7cqPj4eB05ckSvvfbaNbeVkJCgiRMnlnTJAAAAKGU2y7IsVxchXf7w2NKlS9WzZ88819erV0/33Xefpk+fft3tzJ07V0888YTOnTsnu92eZ5+MjAxlZGQ4ltPS0hQeHq7U1FT5+voWeR8AlCNcCqzgJqS6ugIAUFpamvz8/PLNa+Vi5Parr77Snj179MEHH+TbNyoqSpcuXdIvv/yiunXr5tnHbrdfM/gCAACg/CoXc27nzJmjJk2aqHHjxvn23b59u9zc3BQUFFQKlQEAAKAscenI7blz57Rv3z7HcnJysrZv366AgADVqFFD0uUh6MWLF+vVV1/Ndf+kpCRt3rxZbdq0kY+Pj5KSkjRmzBg9+uijqlKlSqntBwAAAMoGl4bb7777Tm3atHEsx8XFSZJiYmI0f/58SdKiRYtkWZb69euX6/52u12LFi3ShAkTlJGRocjISI0ZM8axHQAAAPy5lJkPlLlSQScoAzAIHygrOD5QBqAMKGheKxdzbgEAAICCINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABjDpeF2w4YN6t69u8LCwmSz2bRs2TKn9QMHDpTNZnO6derUyanPqVOn1L9/f/n6+srf319DhgzRuXPnSnEvAAAAUFa4NNymp6ercePGmjFjxjX7dOrUSUeOHHHc/v3vfzut79+/v3bt2qVVq1ZpxYoV2rBhg4YNG1bSpQMAAKAMqujKB+/cubM6d+583T52u10hISF5rvv555+1cuVKbdmyRU2bNpUkTZ8+XV26dNErr7yisLCwYq8ZAP50Jvi5uoK8TUh1dQUAyqAyP+d23bp1CgoKUt26dTV8+HCdPHnSsS4pKUn+/v6OYCtJ7du3l5ubmzZv3uyKcgEAAOBCLh25zU+nTp304IMPKjIyUvv379cLL7ygzp07KykpSRUqVFBKSoqCgoKc7lOxYkUFBAQoJSXlmtvNyMhQRkaGYzktLa3E9gEAAAClp0yH2759+zp+btiwoRo1aqTatWtr3bp1ateuXZG3m5CQoIkTJxZHiQAAAChDyvy0hCvVqlVL1apV0759+yRJISEhOnbsmFOfS5cu6dSpU9ecpytJ8fHxSk1NddwOHTpUonUDAACgdJSrcPvbb7/p5MmTCg0NlSRFR0frzJkz2rp1q6PP2rVrlZ2draioqGtux263y9fX1+kGAACA8s+l0xLOnTvnGIWVpOTkZG3fvl0BAQEKCAjQxIkT1atXL4WEhGj//v167rnnVKdOHXXs2FGSdOutt6pTp04aOnSoZs+erYsXL2rkyJHq27cvV0oAAAD4E3LpyO13332nO+64Q3fccYckKS4uTnfccYfGjRunChUqaMeOHerRo4duueUWDRkyRE2aNNFXX30lu93u2MaCBQtUr149tWvXTl26dFGLFi30zjvvuGqXAAAA4EI2y7IsVxfhamlpafLz81NqaipTFIA/i7J67VYUHNe5Bf5UCprXytWcWwAAAOB6CLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGCMiq4uAMCfwAQ/V1cAAPiTYOQWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjuDTcbtiwQd27d1dYWJhsNpuWLVvmWHfx4kWNHTtWDRs2lJeXl8LCwvTYY4/p8OHDTtuoWbOmbDab023y5MmlvCcAAAAoC1wabtPT09W4cWPNmDEj17rz589r27ZtevHFF7Vt2zYtWbJEe/bsUY8ePXL1nTRpko4cOeK4jRo1qjTKBwAAQBnj0kuBde7cWZ07d85znZ+fn1atWuXU9tZbb6lZs2Y6ePCgatSo4Wj38fFRSEhIidYKAACAsq9czblNTU2VzWaTv7+/U/vkyZNVtWpV3XHHHZo6daouXbp03e1kZGQoLS3N6QYAAIDyr9x8icOFCxc0duxY9evXT76+vo72p556SnfeeacCAgK0ceNGxcfH68iRI3rttdeuua2EhARNnDixNMoGAABAKbJZlmW5ughJstlsWrp0qXr27Jlr3cWLF9WrVy/99ttvWrdunVO4vdrcuXP1xBNP6Ny5c7Lb7Xn2ycjIUEZGhmM5LS1N4eHhSk1Nve62ARQR31CGkjAh1dUVAChFaWlp8vPzyzevlfmR24sXL+qhhx7Sr7/+qrVr1+YbPqOionTp0iX98ssvqlu3bp597Hb7NYMvAAAAyq8yHW5zgu3evXuVmJioqlWr5nuf7du3y83NTUFBQaVQIQAAAMoSl4bbc+fOad++fY7l5ORkbd++XQEBAQoNDVXv3r21bds2rVixQllZWUpJSZEkBQQEyN3dXUlJSdq8ebPatGkjHx8fJSUlacyYMXr00UdVpUoVV+0WAAAAXMSlc27XrVunNm3a5GqPiYnRhAkTFBkZmef9EhMT1bp1a23btk0jRozQ7t27lZGRocjISA0YMEBxcXGFmnZQ0DkcAIqIObcoCcy5Bf5UysWc29atW+t62Tq/3H3nnXdq06ZNxV0WAAAAyqlydZ1bAAAA4HoItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxihRu27ZtqzNnzuRqT0tLU9u2bW+0JgAAAKBIihRu161bp8zMzFztFy5c0FdffXXDRQEAAABFUbEwnXfs2OH4+aefflJKSopjOSsrSytXrtRNN91UfNUBAAAAhVCocHv77bfLZrPJZrPlOf3A09NT06dPL7biAAAAgMIoVLhNTk6WZVmqVauWvv32WwUGBjrWubu7KygoSBUqVCj2IgEAAICCKFS4jYiIkCRlZ2eXSDEAAADAjShUuL3S3r17lZiYqGPHjuUKu+PGjbvhwgAAuK4Jfq6uIG8TUl1dAfCnVqRw++6772r48OGqVq2aQkJCZLPZHOtsNhvhFgAAAC5RpHD7j3/8Qy+99JLGjh1b3PUAAAAARVak69yePn1affr0Ke5aAAAAgBtSpHDbp08fffnll8VdCwAAAHBDijQtoU6dOnrxxRe1adMmNWzYUJUqVXJa/9RTTxVLcQAAAEBh2CzLsgp7p8jIyGtv0GbTgQMHbqio0paWliY/Pz+lpqbK19fX1eUA5imrn2oHSgJXSwBKREHzWpFGbpOTk4tcGAAAAFBSijTnFgAAACiLijRyO3jw4Ouunzt3bpGKAQAAAG5EkcLt6dOnnZYvXryonTt36syZM2rbtm2xFAYAAAAUVpHC7dKlS3O1ZWdna/jw4apdu/YNFwUAAAAURbHNuXVzc1NcXJxef/314tokAAAAUCjF+oGy/fv369KlS8W5SQAAAKDAijQtIS4uzmnZsiwdOXJEn332mWJiYoqlMAAAAKCwihRuv//+e6dlNzc3BQYG6tVXX833SgoAAABASSlSuE1MTCzuOgAAAIAbVqRwm+P48ePas2ePJKlu3boKDAwslqIAAACAoijSB8rS09M1ePBghYaGqmXLlmrZsqXCwsI0ZMgQnT9/vrhrBAAAAAqkSOE2Li5O69ev16effqozZ87ozJkzWr58udavX69nnnmmuGsEAAAACqRI4fbjjz/WnDlz1LlzZ/n6+srX11ddunTRu+++q48++qjA29mwYYO6d++usLAw2Ww2LVu2zGm9ZVkaN26cQkND5enpqfbt22vv3r1OfU6dOqX+/fvL19dX/v7+GjJkiM6dO1eU3QIAAEA5V6Rwe/78eQUHB+dqDwoKKtS0hPT0dDVu3FgzZszIc/2UKVP05ptvavbs2dq8ebO8vLzUsWNHXbhwwdGnf//+2rVrl1atWqUVK1Zow4YNGjZsWOF3CgAAAOWezbIsq7B3ateunapWrap//etf8vDwkCT98ccfiomJ0alTp7R69erCF2KzaenSperZs6eky6O2YWFheuaZZ/Tss89KklJTUxUcHKz58+erb9+++vnnn1W/fn1t2bJFTZs2lSStXLlSXbp00W+//aawsLACPXZaWpr8/PyUmpoqX1/fQtcOIB8T/FxdAVB6JqS6ugLASAXNa0W6WsK0adPUqVMnVa9eXY0bN5Yk/fDDD7Lb7fryyy+LVvFVkpOTlZKSovbt2zva/Pz8FBUVpaSkJPXt21dJSUny9/d3BFtJat++vdzc3LR582Y98MADeW47IyNDGRkZjuW0tLRiqRkAAACuVaRw27BhQ+3du1cLFizQ7t27JUn9+vVT//795enpWSyFpaSkSFKu6Q/BwcGOdSkpKQoKCnJaX7FiRQUEBDj65CUhIUETJ04sljoBAABQdhQp3CYkJCg4OFhDhw51ap87d66OHz+usWPHFktxJSU+Pt7pK4TT0tIUHh7uwooAAABQHIr0gbK3335b9erVy9XeoEEDzZ49+4aLkqSQkBBJ0tGjR53ajx496lgXEhKiY8eOOa2/dOmSTp065eiTF7vd7rjKQ84NAAAA5V+Rwm1KSopCQ0NztQcGBurIkSM3XJQkRUZGKiQkRGvWrHG0paWlafPmzYqOjpYkRUdH68yZM9q6daujz9q1a5Wdna2oqKhiqQMAAADlR5GmJYSHh+ubb75RZGSkU/s333xT4CsUSNK5c+e0b98+x3JycrK2b9+ugIAA1ahRQ6NHj9Y//vEP3XzzzYqMjNSLL76osLAwxxUVbr31VnXq1ElDhw7V7NmzdfHiRY0cOVJ9+/YtVB0AAAAwQ5HC7dChQzV69GhdvHhRbdu2lSStWbNGzz33XKG+oey7775TmzZtHMs582BjYmI0f/58Pffcc0pPT9ewYcN05swZtWjRQitXrnRcfkySFixYoJEjR6pdu3Zyc3NTr1699OabbxZltwAAAFDOFek6t5Zl6fnnn9ebb76pzMxMSZKHh4fGjh2rcePGFXuRJY3r3AIljOvc4s+E69wCJaKgea1I4TbHuXPn9PPPP8vT01M333yz7HZ7UTflUoRboIQRbvFnQrgFSkSJfolDDm9vb9111103sgkAAACg2BTpagkAAABAWUS4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxR5sNtzZo1ZbPZct1iY2MlSa1bt8617sknn3Rx1QAAAHCFiq4uID9btmxRVlaWY3nnzp2677771KdPH0fb0KFDNWnSJMdy5cqVS7VGAAAAlA1lPtwGBgY6LU+ePFm1a9dWq1atHG2VK1dWSEhIaZcGAACAMqbMh9srZWZm6v3331dcXJxsNpujfcGCBXr//fcVEhKi7t2768UXX2T0Fn9OE/xcXQEAAC5VrsLtsmXLdObMGQ0cONDR9sgjjygiIkJhYWHasWOHxo4dqz179mjJkiXX3E5GRoYyMjIcy2lpaSVZNgAAAEpJuQq3c+bMUefOnRUWFuZoGzZsmOPnhg0bKjQ0VO3atdP+/ftVu3btPLeTkJCgiRMnlni9AAAAKF1l/moJOX799VetXr1ajz/++HX7RUVFSZL27dt3zT7x8fFKTU113A4dOlSstQIAAMA1ys3I7bx58xQUFKSuXbtet9/27dslSaGhodfsY7fbZbfbi7M8AAAAlAHlItxmZ2dr3rx5iomJUcWK/yt5//79Wrhwobp06aKqVatqx44dGjNmjFq2bKlGjRq5sGIAAAC4QrkIt6tXr9bBgwc1ePBgp3Z3d3etXr1a06ZNU3p6usLDw9WrVy/97W9/c1GlAAAAcKVyEW47dOggy7JytYeHh2v9+vUuqAgAAABlUbn5QBkAAACQH8ItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGqOjqAgAAMMoEP1dXkNuEVFdXAJQaRm4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYZTrcTpgwQTabzelWr149x/oLFy4oNjZWVatWlbe3t3r16qWjR4+6sGIAAAC4UpkOt5LUoEEDHTlyxHH7+uuvHevGjBmjTz/9VIsXL9b69et1+PBhPfjggy6sFgAAAK5U0dUF5KdixYoKCQnJ1Z6amqo5c+Zo4cKFatu2rSRp3rx5uvXWW7Vp0yY1b968tEsFAACAi5X5kdu9e/cqLCxMtWrVUv/+/XXw4EFJ0tatW3Xx4kW1b9/e0bdevXqqUaOGkpKSrrvNjIwMpaWlOd0AAABQ/pXpcBsVFaX58+dr5cqVmjVrlpKTk3Xvvffq7NmzSklJkbu7u/z9/Z3uExwcrJSUlOtuNyEhQX5+fo5beHh4Ce4FAAAASkuZnpbQuXNnx8+NGjVSVFSUIiIi9OGHH8rT07PI242Pj1dcXJxjOS0tjYALAABggDI9cns1f39/3XLLLdq3b59CQkKUmZmpM2fOOPU5evRonnN0r2S32+Xr6+t0AwAAQPlXrsLtuXPntH//foWGhqpJkyaqVKmS1qxZ41i/Z88eHTx4UNHR0S6sEgAAAK5SpqclPPvss+revbsiIiJ0+PBhjR8/XhUqVFC/fv3k5+enIUOGKC4uTgEBAfL19dWoUaMUHR3NlRIAAAD+pMp0uP3tt9/Ur18/nTx5UoGBgWrRooU2bdqkwMBASdLrr78uNzc39erVSxkZGerYsaNmzpzp4qoBAADgKjbLsixXF+FqaWlp8vPzU2pqKvNvUb5N8HN1BQDKogmprq4AuGEFzWvlas4tAAAAcD2EWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAG4RYAAADGINwCAADAGIRbAAAAGINwCwAAAGNUdHUBQLk1wc/VFQAAgKswcgsAAABjEG4BAABgDMItAAAAjMGcWwAATFdWPyMwIdXVFcBAjNwCAADAGIRbAAAAGINwCwAAAGOU6XCbkJCgu+66Sz4+PgoKClLPnj21Z88epz6tW7eWzWZzuj355JMuqhgAAACuVKbD7fr16xUbG6tNmzZp1apVunjxojp06KD09HSnfkOHDtWRI0cctylTprioYgAAALhSmb5awsqVK52W58+fr6CgIG3dulUtW7Z0tFeuXFkhISGlXR4AAADKmDI9cnu11NTLlwwJCAhwal+wYIGqVaum2267TfHx8Tp//rwrygMAAICLlemR2ytlZ2dr9OjRuueee3Tbbbc52h955BFFREQoLCxMO3bs0NixY7Vnzx4tWbLkmtvKyMhQRkaGYzktLa1EawcAAEDpKDfhNjY2Vjt37tTXX3/t1D5s2DDHzw0bNlRoaKjatWun/fv3q3bt2nluKyEhQRMnTizRegEAAFD6ysW0hJEjR2rFihVKTExU9erVr9s3KipKkrRv375r9omPj1dqaqrjdujQoWKtFwAAAK5RpkduLcvSqFGjtHTpUq1bt06RkZH53mf79u2SpNDQ0Gv2sdvtstvtxVUmAAAAyogyHW5jY2O1cOFCLV++XD4+PkpJSZEk+fn5ydPTU/v379fChQvVpUsXVa1aVTt27NCYMWPUsmVLNWrUyMXVAwAAoLSV6XA7a9YsSZe/qOFK8+bN08CBA+Xu7q7Vq1dr2rRpSk9PV3h4uHr16qW//e1vLqgWAAAArlamw61lWdddHx4ervXr15dSNQAAACjrysUHygAAAICCINwCAADAGIRbAAAAGINwCwAAAGMQbgEAAGAMwi0AAACMQbgFAACAMQi3AAAAMAbhFgAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjFHR1QUA+Zrg5+oKAABAOcHILQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDS4EBAADXKIuXepyQ6uoKcIMYuQUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABjEG4BAABgDMItAAAAjEG4BQAAgDEItwAAADAGX7/rKnzlIAAAZU9Z/Pss8Te6EBi5BQAAgDEYucX/lNX/VgEAAAqIkVsAAAAYw5iR2xkzZmjq1KlKSUlR48aNNX36dDVr1szVZQEAANy4svjuahmdB2zEyO0HH3yguLg4jR8/Xtu2bVPjxo3VsWNHHTt2zNWlAQAAoBQZEW5fe+01DR06VIMGDVL9+vU1e/ZsVa5cWXPnznV1aQAAAChF5X5aQmZmprZu3ar4+HhHm5ubm9q3b6+kpKQ875ORkaGMjAzHcmrq5WH1tLS0ki3WqQir9B4LAACguJVmbtL/cpplXT9Dlftwe+LECWVlZSk4ONipPTg4WLt3787zPgkJCZo4cWKu9vDw8BKpEQAAwDiTXTMP+OzZs/Lzu/Zjl/twWxTx8fGKi4tzLGdnZ+vUqVOqWrWqbDabCysrfWlpaQoPD9ehQ4fk6+vr6nLKJY5h8eA4Fg+OY/HgOBYPjuON4xj+j2VZOnv2rMLCwq7br9yH22rVqqlChQo6evSoU/vRo0cVEhKS533sdrvsdrtTm7+/f0mVWC74+vr+6X9pbhTHsHhwHIsHx7F4cByLB8fxxnEML7veiG2Ocv+BMnd3dzVp0kRr1qxxtGVnZ2vNmjWKjo52YWUAAAAobeV+5FaS4uLiFBMTo6ZNm6pZs2aaNm2a0tPTNWjQIFeXBgAAgFJkRLh9+OGHdfz4cY0bN04pKSm6/fbbtXLlylwfMkNudrtd48ePzzVNAwXHMSweHMfiwXEsHhzH4sFxvHEcw8KzWfldTwEAAAAoJ8r9nFsAAAAgB+EWAAAAxiDcAgAAwBiEWwAAABiDcAuHHj16qEaNGvLw8FBoaKgGDBigw4cPu7qscuWXX37RkCFDFBkZKU9PT9WuXVvjx49XZmamq0srV1566SXdfffdqly58p/+C1YKY8aMGapZs6Y8PDwUFRWlb7/91tUllTsbNmxQ9+7dFRYWJpvNpmXLlrm6pHInISFBd911l3x8fBQUFKSePXtqz549ri6r3Jk1a5YaNWrk+PKG6Ohoff75564uq1wg3MKhTZs2+vDDD7Vnzx59/PHH2r9/v3r37u3qssqV3bt3Kzs7W2+//bZ27dql119/XbNnz9YLL7zg6tLKlczMTPXp00fDhw93dSnlxgcffKC4uDiNHz9e27ZtU+PGjdWxY0cdO3bM1aWVK+np6WrcuLFmzJjh6lLKrfXr1ys2NlabNm3SqlWrdPHiRXXo0EHp6emuLq1cqV69uiZPnqytW7fqu+++U9u2bXX//fdr165dri6tzONSYLimTz75RD179lRGRoYqVark6nLKralTp2rWrFk6cOCAq0spd+bPn6/Ro0frzJkzri6lzIuKitJdd92lt956S9Llb2oMDw/XqFGj9Pzzz7u4uvLJZrNp6dKl6tmzp6tLKdeOHz+uoKAgrV+/Xi1btnR1OeVaQECApk6dqiFDhri6lDKNkVvk6dSpU1qwYIHuvvtugu0NSk1NVUBAgKvLgMEyMzO1detWtW/f3tHm5uam9u3bKykpyYWVAZdfAyXxOngDsrKytGjRIqWnpys6OtrV5ZR5hFs4GTt2rLy8vFS1alUdPHhQy5cvd3VJ5dq+ffs0ffp0PfHEE64uBQY7ceKEsrKycn0rY3BwsFJSUlxUFXD5HYTRo0frnnvu0W233ebqcsqdH3/8Ud7e3rLb7XryySe1dOlS1a9f39VllXmEW8M9//zzstls173t3r3b0f8vf/mLvv/+e3355ZeqUKGCHnvsMTFzpfDHUZJ+//13derUSX369NHQoUNdVHnZUZRjCKB8i42N1c6dO7Vo0SJXl1Iu1a1bV9u3b9fmzZs1fPhwxcTE6KeffnJ1WWUec24Nd/z4cZ08efK6fWrVqiV3d/dc7b/99pvCw8O1cePGP/3bIIU9jocPH1br1q3VvHlzzZ8/X25u/B9ZlHORObcFk5mZqcqVK+ujjz5ymh8aExOjM2fO8A5METHn9saMHDlSy5cv14YNGxQZGenqcozQvn171a5dW2+//barSynTKrq6AJSswMBABQYGFum+2dnZkqSMjIziLKlcKsxx/P3339WmTRs1adJE8+bNI9j+/27kXMT1ubu7q0mTJlqzZo0jiGVnZ2vNmjUaOXKka4vDn45lWRo1apSWLl2qdevWEWyLUXZ2Nn+TC4BwC0nS5s2btWXLFrVo0UJVqlTR/v379eKLL6p27dp/+lHbwvj999/VunVrRURE6JVXXtHx48cd60JCQlxYWfly8OBBnTp1SgcPHlRWVpa2b98uSapTp468vb1dW1wZFRcXp5iYGDVt2lTNmjXTtGnTlJ6erkGDBrm6tHLl3Llz2rdvn2M5OTlZ27dvV0BAgGrUqOHCysqP2NhYLVy4UMuXL5ePj49j3refn588PT1dXF35ER8fr86dO6tGjRo6e/asFi5cqHXr1umLL75wdWllnwVYlrVjxw6rTZs2VkBAgGW3262aNWtaTz75pPXbb7+5urRyZd68eZakPG8ouJiYmDyPYWJioqtLK9OmT59u1ahRw3J3d7eaNWtmbdq0ydUllTuJiYl5nnsxMTGuLq3cuNZr4Lx581xdWrkyePBgKyIiwnJ3d7cCAwOtdu3aWV9++aWryyoXmHMLAAAAYzAZEAAAAMYg3AIAAMAYhFsAAAAYg3ALAAAAYxBuAQAAYAzCLQAAAIxBuAUAAIAxCLcAAAAwBuEWAAAAxiDcAgAAwBiEWwAAABiDcAsAAABj/H/TPNJABGeeJwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(8, 6))\n", "ax.hist(pm.draw(x, draws=1_000), color=\"C1\", bins=15)\n", "ax.set(title=\"Samples from a normal distribution using pymc\", ylabel=\"count\");" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Yay! We learned how to sample from a `pymc` distribution!" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### What is going on behind the scenes?" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can now look into how this is done inside a {class}`~pymc.Model`." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normal_rv{0, (0, 0), floatX, False}.1 [id A] 'z'\n", " |RandomGeneratorSharedVariable() [id B]\n", " |TensorConstant{[]} [id C]\n", " |TensorConstant{11} [id D]\n", " |TensorConstant{(2,) of 0} [id E]\n", " |TensorConstant{[1. 2.]} [id F]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with pm.Model() as model:\n", " z = pm.Normal(name=\"z\", mu=np.array([0, 0]), sigma=np.array([1, 2]))\n", "\n", "pytensor.dprint(z)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We are just creating random variables like we saw before, but now registering them in a `pymc` model. To extract the list of random variables we can simply do:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[z ~ N(, )]" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.basic_RVs" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "normal_rv{0, (0, 0), floatX, False}.1 [id A] 'z'\n", " |RandomGeneratorSharedVariable() [id B]\n", " |TensorConstant{[]} [id C]\n", " |TensorConstant{11} [id D]\n", " |TensorConstant{(2,) of 0} [id E]\n", " |TensorConstant{[1. 2.]} [id F]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pytensor.dprint(model.basic_RVs[0])" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can try to sample via {meth}`~pytensor.graph.basic.Variable.eval` as above and it is no surprise that we are getting the same samples at each iteration." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sample 0: [-0.30775592 1.21469108]\n", "Sample 1: [-0.30775592 1.21469108]\n", "Sample 2: [-0.30775592 1.21469108]\n", "Sample 3: [-0.30775592 1.21469108]\n", "Sample 4: [-0.30775592 1.21469108]\n", "Sample 5: [-0.30775592 1.21469108]\n", "Sample 6: [-0.30775592 1.21469108]\n", "Sample 7: [-0.30775592 1.21469108]\n", "Sample 8: [-0.30775592 1.21469108]\n", "Sample 9: [-0.30775592 1.21469108]\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {z.eval()}\")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Again, the correct way of sampling is via {func}`~pymc.draw`. " ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sample 0: [-1.2390824 0.3744465]\n", "Sample 1: [0.76748461 0.95086347]\n", "Sample 2: [-1.11985098 -1.94744586]\n", "Sample 3: [-0.62003335 0.10075427]\n", "Sample 4: [-0.75744869 0.69140323]\n", "Sample 5: [-0.95472672 -1.0814984 ]\n", "Sample 6: [-0.81052179 -2.05414581]\n", "Sample 7: [0.37456894 1.76040717]\n", "Sample 8: [-0.61006854 -0.05034957]\n", "Sample 9: [1.19039658 1.10460999]\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {pm.draw(z)}\")" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAKqCAYAAACTnV4oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4iElEQVR4nO3deZhcZZ03/F91d7rT2ZoQOmxpsoCiCQjzJCGyDCaQYXlZnxnQR0FDZOJCwqqOCTMaUDE4cDlKkIC8CswrDKCIMIMsARM3tkCEETHBCJEYhCQEOntvdd4/fOixTQgJcN/VaT6f6zpXUqdP3b+7zqk69a27Tp1TKoqiCAAASKyq0h0AAOCdQfAEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIQvAE3nFKpVJcdNFFle7GVp1xxhkxbNiwSncD4G0leAJvyq9//es45ZRTYujQodG7d+/Yc8894+/+7u9i9uzZle5at7N06dIolUpx+eWXb/HvF110UZRKpVi1atVbqvP000/HRRddFEuXLn1L7QCkIngC2+3BBx+MMWPGxJNPPhlTpkyJK6+8Mv7xH/8xqqqq4pvf/Galu9cjXHvttbF48eLtus/TTz8dF198seAJdFs1le4AsOO55JJLoqGhIRYsWBA77bRTl7+tWLGiMp3qYXr16lXpLmy3TZs2RW1tbVRVGdMAtszeAdhuv//972PUqFGbhc6IiMGDB3e5fd1118URRxwRgwcPjrq6uhg5cmTMmTNns/sNGzYsjj/++Jg/f36MGTMm6uvrY//994/58+dHRMQPf/jD2H///aN3794xevTo+NWvftXl/meccUb069cvnn322Tj66KOjb9++sccee8SXvvSlKIriDR/T8uXL4+Mf/3jsuuuuUVdXF6NGjYrvfve7my03e/bsGDVqVPTp0ycGDhwYY8aMiZtuuukN299eWzrG8+abb47Ro0dH//79Y8CAAbH//vt3jjBff/31ceqpp0ZExIQJE6JUKkWpVOpcfxERV111VYwaNSrq6upijz32iKlTp8arr766We1vfetbMWLEiKivr4+DDjoofv7zn8f48eNj/PjxncvMnz8/SqVS3HzzzfEv//Ivseeee0afPn1izZo1sXr16vjsZz8b+++/f/Tr1y8GDBgQxx57bDz55JNd6rzWxq233hoXX3xx7LnnntG/f/845ZRTorm5OVpaWuK8886LwYMHR79+/WLy5MnR0tLytqxfoDKMeALbbejQofHQQw/FU089Ffvtt99Wl50zZ06MGjUqTjzxxKipqYn//M//jLPOOivK5XJMnTq1y7JLliyJj3zkI/HJT34yTj/99Lj88svjhBNOiKuvvjouvPDCOOussyIiYtasWfHBD34wFi9e3GV0raOjI4455ph4//vfH//6r/8a99xzT8ycOTPa29vjS1/60uv28aWXXor3v//9USqVYtq0adHY2Bh33313nHnmmbFmzZo477zzIuLPX3+fc845ccopp8S5554bmzZtiv/+7/+ORx55JD7ykY+84XrbsGHDFo/j3LBhwxved+7cufHhD384jjzyyPja174WERG//e1v45e//GWce+65cfjhh8c555wTV1xxRVx44YXx3ve+NyKi89+LLrooLr744pg4cWJ8+tOfjsWLF8ecOXNiwYIF8ctf/rJzhHXOnDkxbdq0+Nu//ds4//zzY+nSpXHyySfHwIEDY8iQIZv168tf/nLU1tbGZz/72WhpaYna2tp4+umn40c/+lGceuqpMXz48HjppZfimmuuiQ984APx9NNPxx577NGljVmzZkV9fX1Mnz49lixZErNnz45evXpFVVVVvPLKK3HRRRfFww8/HNdff30MHz48vvjFL77h+gK6qQJgO913331FdXV1UV1dXRx88MHFP/3TPxX33ntv0drautmyGzZs2Gze0UcfXYwYMaLLvKFDhxYRUTz44IOd8+69994iIor6+vriD3/4Q+f8a665poiIYt68eZ3zJk2aVEREcfbZZ3fOK5fLxXHHHVfU1tYWK1eu7JwfEcXMmTM7b5955pnF7rvvXqxatapLn/7P//k/RUNDQ+djOOmkk4pRo0a9wdrZ3HPPPVdExBtOf9nHSZMmFUOHDu28fe655xYDBgwo2tvbX7fO97///c3WS1EUxYoVK4ra2triqKOOKjo6OjrnX3nllUVEFN/97neLoiiKlpaWYtCgQcXYsWOLtra2zuWuv/76IiKKD3zgA53z5s2bV0REMWLEiM228aZNm7rUeW0d1NXVFV/60pc2a2O//fbr8tz58Ic/XJRKpeLYY4/t0sbBBx/cZZ0AOx5ftQPb7e/+7u/ioYceihNPPDGefPLJ+Nd//dc4+uijY88994w777yzy7L19fWd/29ubo5Vq1bFBz7wgXj22Wejubm5y7IjR46Mgw8+uPP2uHHjIiLiiCOOiL322muz+c8+++xmfZs2bVrn/18bwWxtbY37779/i4+lKIq47bbb4oQTToiiKGLVqlWd09FHHx3Nzc2xcOHCiIjYaaed4o9//GMsWLBgm9bTX/vEJz4Rc+fO3Wz66Ec/+ob33WmnnWL9+vUxd+7c7a57//33R2tra5x33nldRoinTJkSAwYMiLvuuisiIh577LF4+eWXY8qUKVFT8z9fiJ122mkxcODALbY9adKkLts4IqKurq6zTkdHR7z88svRr1+/2HfffTvX5V/62Mc+1uWY1nHjxkVRFPHxj3+8y3Ljxo2LZcuWRXt7+3auAaC78FU78KaMHTs2fvjDH0Zra2s8+eSTcfvtt8e//du/xSmnnBJPPPFEjBw5MiIifvnLX8bMmTPjoYce2uwr5ebm5mhoaOi8/ZfhMiI6/9bU1LTF+a+88kqX+VVVVTFixIgu89797ndHRLzuL71XrlwZr776anz729+Ob3/721tc5rUfTH3+85+P+++/Pw466KDYZ5994qijjoqPfOQjceihh27xfn/tXe96V0ycOHGz+b/4xS/e8L5nnXVW3HrrrXHsscfGnnvuGUcddVR88IMfjGOOOeYN7/uHP/whIiL23XffLvNra2tjxIgRnX9/7d999tmny3I1NTWve07R4cOHbzavXC7HN7/5zbjqqqviueeei46Ojs6/DRo0aLPlt2e7l8vlaG5u3mI7QPdnxBN4S2pra2Ps2LHx1a9+NebMmRNtbW3x/e9/PyL+/COkI488MlatWhVf//rX46677oq5c+fG+eefHxF/Dih/qbq6eos1Xm9+sQ0/Gnojr/Xh9NNP3+Jo5Ny5czuD5Xvf+95YvHhx3HzzzXHYYYfFbbfdFocddljMnDnzLffjjQwePDieeOKJuPPOO+PEE0+MefPmxbHHHhuTJk1KXntr/nq0MyLiq1/9alxwwQVx+OGHx/e+97249957Y+7cuTFq1KjNtnlEZbY7UBlGPIG3zZgxYyIi4k9/+lNERPznf/5ntLS0xJ133tllVGvevHlJ6pfL5Xj22Wc7RzkjIp555pmIiNcdsWtsbIz+/ftHR0fHFkcj/1rfvn3jQx/6UHzoQx+K1tbW+Pu///u45JJLYsaMGdG7d++35XG8ntra2jjhhBPihBNOiHK5HGeddVZcc8018YUvfCH22WefKJVKW7zf0KFDIyJi8eLFXUaEW1tb47nnnut83K8tt2TJkpgwYULncu3t7bF06dJ43/vet039/MEPfhATJkyI73znO13mv/rqq7HLLrts+wMGehwjnsB2mzdv3hZHnX784x9HxP98pfvaiNVfLtvc3BzXXXddsr5deeWVnf8viiKuvPLK6NWrVxx55JFbXL66ujr+4R/+IW677bZ46qmnNvv7ypUrO///8ssvd/lbbW1tjBw5MoqiiLa2trfpEWzZX9euqqrqDIKvnWKob9++ERGbnSJp4sSJUVtbG1dccUWXbfGd73wnmpub47jjjouIP39wGDRoUFx77bVdjqO88cYbNzusYWuqq6s3e358//vfj+XLl29zG0DPZMQT2G5nn312bNiwIf73//7f8Z73vCdaW1vjwQcfjFtuuSWGDRsWkydPjoiIo446qnOU7pOf/GSsW7curr322hg8eHDnqOjbqXfv3nHPPffEpEmTYty4cXH33XfHXXfdFRdeeGE0Nja+7v0uvfTSmDdvXowbNy6mTJkSI0eOjNWrV8fChQvj/vvvj9WrV3c+nt122y0OPfTQ2HXXXeO3v/1tXHnllXHcccdF//793/bH85f+8R//MVavXh1HHHFEDBkyJP7whz/E7Nmz48ADD+w8ZdKBBx4Y1dXV8bWvfS2am5ujrq6u8xyqM2bMiIsvvjiOOeaYOPHEE2Px4sVx1VVXxdixY+P000+PiD8H6YsuuijOPvvsOOKII+KDH/xgLF26NK6//vrYe++9X3dE9a8df/zx8aUvfSkmT54chxxySPz617+OG2+8cbPjb4F3HiOewHa7/PLLY8KECfHjH/84Lrjggrjgggvi0UcfjbPOOiseeeSRzhPL77vvvvGDH/wgSqVSfPazn42rr746PvGJT8S5556bpF/V1dVxzz33xIsvvhif+9znYsGCBTFz5sz48pe/vNX77brrrvHoo4/G5MmT44c//GFMmzYtvvnNb8bq1as7z5kZEZ3h+etf/3pMnTo1fvSjH8U555wT3/ve95I8nr90+umnR+/eveOqq66Ks846K2644Yb40Ic+FHfffXfnL8h32223uPrqq2PFihVx5plnxoc//OF4+umnI+LP5/G88sor4/nnn4/zzz8/br311vjEJz4R9913X5dflE+bNi2uuOKKeP755+Ozn/1s/PznP48777wzdtppp20+lODCCy+Mz3zmM3HvvffGueeeGwsXLoy77rprsx8LAe88pcJR2kAPcMYZZ8QPfvCDWLduXaW70uOUy+VobGyMv//7v49rr7220t0BdmBGPAHotGnTps2Oz/z3f//3WL16dZdLZgK8GY7xBKDTww8/HOeff36ceuqpMWjQoFi4cGF85zvfif3226/zWvAAb5bgCUCnYcOGRVNTU1xxxRWxevXq2HnnneNjH/tYXHrppVFbW1vp7gE7OMd4AgCQhWM8AQDIQvAEACCLbn2MZ7lcjhdeeCH69++/zScuBgAgn6IoYu3atbHHHnt0nlf49XTr4PnCCy844TAAwA5g2bJlMWTIkK0u062D52uXoDss/p+oiV5vsDQAb1opw5FXRTl9DSC79miLX8SPt+nSwd06eL729XpN9IqakuAJkEyO4BmCJ/RI//f8SNtyWKQfFwEAkIXgCQBAFoInAABZCJ4AAGSRNHh2dHTEF77whRg+fHjU19fH3nvvHV/+8pfDVToBAN55kv6q/Wtf+1rMmTMnbrjhhhg1alQ89thjMXny5GhoaIhzzjknZWkAALqZpMHzwQcfjJNOOimOO+64iIgYNmxY/Md//Ec8+uijKcsCANANJf2q/ZBDDokHHnggnnnmmYiIePLJJ+MXv/hFHHvssSnLAgDQDSUd8Zw+fXqsWbMm3vOe90R1dXV0dHTEJZdcEqeddtoWl29paYmWlpbO22vWrEnZPQAAMko64nnrrbfGjTfeGDfddFMsXLgwbrjhhrj88svjhhtu2OLys2bNioaGhs7JddoBAHqOUpHwJ+ZNTU0xffr0mDp1aue8r3zlK/G9730vFi1atNnyWxrxbGpqivFxkktmAqTkWu3Am9RetMX8uCOam5tjwIABW1026VftGzZsiKqqrjuz6urqKJe3vPOpq6uLurq6lF0CAKBCkgbPE044IS655JLYa6+9YtSoUfGrX/0qvv71r8fHP/7xlGUBAOiGkgbP2bNnxxe+8IU466yzYsWKFbHHHnvEJz/5yfjiF7+YsiwAAN1Q0mM836o1a9ZEQ0ODYzwBUnOMJ/Ambc8xnq7VDgBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGSR9ATyAOwgnGOze0l9XlXbmwox4gkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJBFTaU7ALmVanolr1G0tyWvEaXEnxuLctr2I9I/BrZdhu2d5bXX0ZG8RpbXRo4aUAH2+gAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWNZXuAORWtLelL1JK/5muVF2duELq9vNsi1JNr+Q1cig6OtIWyPCcBbCnAQAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALJIHjyXL18ep59+egwaNCjq6+tj//33j8ceeyx1WQAAupmkJ5B/5ZVX4tBDD40JEybE3XffHY2NjfG73/0uBg4cmLIsAADdUNLg+bWvfS2ampriuuuu65w3fPjwlCUBAOimkn7Vfuedd8aYMWPi1FNPjcGDB8ff/M3fxLXXXvu6y7e0tMSaNWu6TAAA9AxJg+ezzz4bc+bMiXe9611x7733xqc//ek455xz4oYbbtji8rNmzYqGhobOqampKWX3AADIqFQURZGq8dra2hgzZkw8+OCDnfPOOeecWLBgQTz00EObLd/S0hItLS2dt9esWRNNTU0xPk6KmlKvVN2Et18p/QkjStXVyWukVrS3Ja9RqukZ+46io6PSXXjLcjxns6ynopy+BuxA2ou2mB93RHNzcwwYMGCryyZ9d9x9991j5MiRXea9973vjeeff36Ly9fV1cWAAQO6TAAA9AxJg+ehhx4aixcv7jLvmWeeiaFDh6YsCwBAN5Q0eJ5//vnx8MMPx1e/+tVYsmRJ3HTTTfHtb387pk6dmrIsAADdUNLgOXbs2Lj99tvjP/7jP2K//faLL3/5y/GNb3wjTjvttJRlAQDohpKexzMi4vjjj4/jjz8+dRkAALo512oHACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCySn06JHiTD9cezyHCd5RzXpC5VldK2X1ubtP2IiI71O/71xyMiqvr1TV4j9XOqvHZt0vYjIkp9+iSvUWzYkLxGltfGunXJa0Al9JAkAQBAdyd4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWNZXuwDtGKUPGL8pJm6/qlf7pUqqtTV6jaG1NXqPc1p68RlV936Tt59gW1TU9YxdUqu+dvkhV2n1IVe+6pO1HRBRr1qavUS6S14iOjuQlqhK//nLso0rV1clrFBm2Rer31oiIUk2vpO0X7W1J298eRjwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyqKl0B94xinLyEtX9+iVtv2htTdp+RER5U0vyGlX1vZPXqGncJXmNKCd+TvXtk7b9iIhe6XdB7YPSP45eL61NXiM2bEzbfnX6cYjSgP7Ja1QPbEheo/ziiuQ1inKRukDa9iOiVJX+9V20p38cORTtbWkLlFK/vqsitvEpa8QTAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyCJb8Lz00kujVCrFeeedl6skAADdSJbguWDBgrjmmmvife97X45yAAB0Q8mD57p16+K0006La6+9NgYOHJi6HAAA3VTy4Dl16tQ47rjjYuLEiW+4bEtLS6xZs6bLBABAz5D0Qqo333xzLFy4MBYsWLBNy8+aNSsuvvjilF0CAKBCko14Llu2LM4999y48cYbo3fv3tt0nxkzZkRzc3PntGzZslTdAwAgs2Qjno8//nisWLEi/tf/+l+d8zo6OuJnP/tZXHnlldHS0hLV1dVd7lNXVxd1dXWpugQAQAUlC55HHnlk/PrXv+4yb/LkyfGe97wnPv/5z28WOgEA6NmSBc/+/fvHfvvt12Ve3759Y9CgQZvNBwCg53PlIgAAskj6q/a/Nn/+/JzlAADoRox4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWWU+n9KaVqv487cCqeqe/FGjH+g1J26/u2ydp+xERpVIpfY1+fZPXiPreyUuU+yauUZ1+W5Rr0++CNuyefltU7ZL+9d1rXUfS9osM27v2lU3Ja1Rtak9fY2P6x5Faed365DWKtvTbIkc2KGW40mLRkfb1HUW527S/Y6c5AAB2GIInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGRRU+kObJOiHBHldO2X0ufvoq09eY2eoFRbm7xGsdOA5DXK/euS12hprE/a/sadq5O2HxHRUVtKXmPj4OQlsui9KvH2SL8pohjeK3mNnX+7KXmNmlL6J1XV+rSPo5ThPalUnX4fUrS2Jq9RzlAjueQ5pyqi2OYlAQAgPcETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyKKm0h14pyj1Sr+qq3rXpS1QKqVtPyKiKn2Njob65DXWDUtfY+POaT83tvdJ2nxERGzYs0heozywNXmN2j5tyWtsqulI2n7r4oak7UdE1K9MXiJaBvZKXqPcK/2YTV1H2u1d2q0xafsREeVn/5C8Rqku8fteRERbe/oaqRXlbtO+EU8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgi6TBc9asWTF27Njo379/DB48OE4++eRYvHhxypIAAHRTSYPnT3/605g6dWo8/PDDMXfu3Ghra4ujjjoq1q9fn7IsAADdUNLL6dxzzz1dbl9//fUxePDgePzxx+Pwww9PWRoAgG4m6zGezc3NERGx88475ywLAEA3kO1a7eVyOc4777w49NBDY7/99tviMi0tLdHS0tJ5e82aNbm6BwBAYtlGPKdOnRpPPfVU3Hzzza+7zKxZs6KhoaFzampqytU9AAASyxI8p02bFv/1X/8V8+bNiyFDhrzucjNmzIjm5ubOadmyZTm6BwBABkm/ai+KIs4+++y4/fbbY/78+TF8+PCtLl9XVxd1dXUpuwQAQIUkDZ5Tp06Nm266Ke64447o379/vPjiixER0dDQEPX19SlLAwDQzST9qn3OnDnR3Nwc48ePj913371zuuWWW1KWBQCgG0r+VTsAAES4VjsAAJkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZJH0dEo7ilJ1dfIa5U0tyWvUDN4lbYHe6a8qVR40IHmNtoba5DWqN6U/ldi6prTtFzUZToc2eFPyEjsN2Ji8Rmtb+n3IGe9+JGn7V714ZNL2IyJaMqynDYPT16jeKf2YTc263mnbX9GctP2IiFKOKxFmeP/OkRGK9rbkNboLI54AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFjWV7kB3UHR0JK9R1Sv9qi42bEzafqkq/eeUqrWbktco79k3eY1X9q1OXiNKaZtvOvCFtAUi4tVN9clrnD7i0eQ1epXS70Nayr2Stv/RQ3+ZtP2IiNt+f0DyGuuLhuQ1Gn6fvESUyuW0BWrTPp8iIkq1tclrFC0tyWuUqhLvbCOiSF6h+zDiCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkEVNpTuwTUpVf55SNV9dnaztzhq1telr9K5LW6Bf37TtR0THwD7pa9T2jM9b7QPbk7b/wisNSduPiFh02P+XvEY5yslrPNe+IXmNteW0u+v/d9XhSduPiOhdm/Y5GxHxar8ieY1SR/IS0d437XtG9YuvJG0/IiKqSslLdKxP/9qr6pUhKiXMOHlURWzjS29Hf6QAAOwgBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyCJL8PzWt74Vw4YNi969e8e4cePi0UcfzVEWAIBuJHnwvOWWW+KCCy6ImTNnxsKFC+OAAw6Io48+OlasWJG6NAAA3Ujy4Pn1r389pkyZEpMnT46RI0fG1VdfHX369Invfve7qUsDANCNJA2era2t8fjjj8fEiRP/p2BVVUycODEeeuihlKUBAOhmkl6AdNWqVdHR0RG77rprl/m77rprLFq0aLPlW1paoqWlpfP2mjVrUnYPAICMutWv2mfNmhUNDQ2dU1NTU6W7BADA2yRp8Nxll12iuro6XnrppS7zX3rppdhtt902W37GjBnR3NzcOS1btixl9wAAyChp8KytrY3Ro0fHAw880DmvXC7HAw88EAcffPBmy9fV1cWAAQO6TAAA9AxJj/GMiLjgggti0qRJMWbMmDjooIPiG9/4Rqxfvz4mT56cujQAAN1I8uD5oQ99KFauXBlf/OIX48UXX4wDDzww7rnnns1+cAQAQM+WPHhGREybNi2mTZuWoxQAAN1Ut/pVOwAAPZfgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFllOp/SWFeWIKKdrviNZ0/9To7U1eY1SXV3aAhkeQ9Wm2uQ1ajaley511tiQ4TNduZS+RmKf+OMhyWscMmBJ8hrPtTQmr/Hshl2Stl9TSv+6WJ/h9Z1jOKW6tUheo6hO/PquzbAtWjK8Z/RO/L4XEeWNG5PXSK7UfcYZu09PAADo0QRPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMiiptId2Balml5RKvVK1n7R0ZGs7ZyK1ta0BV5uSdt+RJR6pdvOr2mvG5C8Rjn9w4g+z1cnbX9DfW3S9iMifvH8iOQ12oakXU8REWvaeievsXhlY9L2e9e2J20/IiJ+0z95ieoM72q9NqR/z6jZkHZ7lPumf85Wrd+QvEYOVbXp94Xltgyvv27CiCcAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZCF4AgCQheAJAEAWNZXuwLYo2tuiKKVrv6q2Nl3j/1e5rT15jarq6rTtNw5K2n5ERNSkfQwREX2f35C8Ru2auuQ1Xnl32udt7YL0r4tX35N+F/TTNe9OXqOqtiN5jY6NaddVyyvpt0XDH5OXiPqX02+L2ldak9eofiXxfupPK9K2HxFFe/r3vfLGjclrRCnDGF1RTl8jpe3ovxFPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIItkwXPp0qVx5plnxvDhw6O+vj723nvvmDlzZrS2pj//GQAA3U+yMwYvWrQoyuVyXHPNNbHPPvvEU089FVOmTIn169fH5ZdfnqosAADdVLLgecwxx8QxxxzTeXvEiBGxePHimDNnjuAJAPAOlPUYz+bm5th5551zlgQAoJvIdq32JUuWxOzZs7c62tnS0hItLS2dt9esWZOjawAAZLDdI57Tp0+PUqm01WnRokVd7rN8+fI45phj4tRTT40pU6a8btuzZs2KhoaGzqmpqWn7HxEAAN1SqSiKYnvusHLlynj55Ze3usyIESOitrY2IiJeeOGFGD9+fLz//e+P66+/PqqqXj/rbmnEs6mpKcbHSVFT6rU93dwuVf+3rymV29qT16jqXZe2/cZBSduPiIg+9clLdPTrnbxG205pt0VExCvvTvy8LaVtPiLi1feU0xfpm+G1V9uRvEbHxrRfUNW8kv4LsIbfpX9S1b+cflv0ebHljRd6i6pf2ZC2wJ9WpG0/IqI9/WuvY9265DWilOGoxCLDvjCh9qIt5scd0dzcHAMGDNjqstu9p2lsbIzGxsZtWnb58uUxYcKEGD16dFx33XVbDZ0REXV1dVFXl/4NGwCA/JJ9xF2+fHmMHz8+hg4dGpdffnmsXLmy82+77bZbqrIAAHRTyYLn3LlzY8mSJbFkyZIYMmRIl79t57f7AAD0AMkOXDjjjDOiKIotTgAAvPO4VjsAAFkIngAAZCF4AgCQheAJAEAWgicAAFkIngAAZJH+Gmlvh1JV0ktWlVtbk7XdkxRr0l+arJTh8qXVq9emr7FuU/IaA6MhafsvjU2/LRqeSf/Zt6MuwyVx013Rt1Ov9WnbH7A0/eUN1w5J/5bTe1X6/XnV+vQ1Sqkvs9wr/bboWJf4ScsOyYgnAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkEVNpTuwTYpyRJQr3Yu3ppQ+45c3bkzaflVHR9L2IyJiVfoSRUtL8hpVQ4ckr1G9qT1p+7s+lrT5iIjo1Zx+W7Q11CWv0dq/OnmNcm0pafu1r7YlbT8iYuDG9PvxXqs3JK9R2pR+XcXadWnbLxdp24+IUq/0EaPI8b5U7OD5o5sx4gkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJCF4AkAQBaCJwAAWQieAABkIXgCAJBFTaU7sE1KVX+eUinK6drOWSOxolwkr1FuXpO8RlXvuuQ14pXm5CV6bdiYtv0+9Unbj4goatPvgmpbOpLXqFlfm7xGkXiYoNeL6V97sXFT+hrlDPvaXumft8X6DWnb70j/uogcNdjhGPEEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIIkvwbGlpiQMPPDBKpVI88cQTOUoCANDNZAme//RP/xR77LFHjlIAAHRTyYPn3XffHffdd19cfvnlqUsBANCNJb3g7EsvvRRTpkyJH/3oR9GnT583XL6lpSVaWlo6b69Zk+HawQAAZJFsxLMoijjjjDPiU5/6VIwZM2ab7jNr1qxoaGjonJqamlJ1DwCAzLY7eE6fPj1KpdJWp0WLFsXs2bNj7dq1MWPGjG1ue8aMGdHc3Nw5LVu2bHu7BwBAN7XdX7V/5jOfiTPOOGOry4wYMSJ+8pOfxEMPPRR1dXVd/jZmzJg47bTT4oYbbtjsfnV1dZstDwBAz7DdwbOxsTEaGxvfcLkrrrgivvKVr3TefuGFF+Loo4+OW265JcaNG7e9ZQEA2MEl+3HRXnvt1eV2v379IiJi7733jiFDhqQqCwBAN+XKRQAAZJH0dEp/adiwYVEURa5yAAB0M0Y8AQDIQvAEACALwRMAgCwETwAAssj246J3vFKGjF+U0zbf3pa0/YjIsp6KtvbkNTpefiV5jeo9d0tbYOXqtO1HRKmqlL5Gfe/kNar+uDF5jaK1NW37GX78mfoxRESU+vRJXqN4eUPyGlFdnbT58sb0z9me8L4XET3ncXQTRjwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyEDwBAMhC8AQAIAvBEwCALARPAACyqKl0B7ZJUY6IcqV78dYUO3j/e5BSVanSXXhblP/0UtL2i3KRtP2IiFKvDLugdeuTlyjV7Bi70neCcobtXbS3Ja9Rqum1Q7cfEVF0dCSvkYX377eVEU8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsBE8AALIQPAEAyELwBAAgC8ETAIAsairdAXYgpZ7xOaXc2pq+SI51VZTTtp/hMRQdHclrVPWuS16jvHFT8hqlXml310Vbe9L2I/Js7+Svi4g8r432trQFesI+ih1Sz0gSAAB0e4InAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkEXS4HnXXXfFuHHjor6+PgYOHBgnn3xyynIAAHRjyS6Fcdttt8WUKVPiq1/9ahxxxBHR3t4eTz31VKpyAAB0c0mCZ3t7e5x77rlx2WWXxZlnntk5f+TIkSnKAQCwA0jyVfvChQtj+fLlUVVVFX/zN38Tu+++exx77LFvOOLZ0tISa9as6TIBANAzJAmezz77bEREXHTRRfEv//Iv8V//9V8xcODAGD9+fKxevfp17zdr1qxoaGjonJqamlJ0DwCACtiu4Dl9+vQolUpbnRYtWhTlcjkiIv75n/85/uEf/iFGjx4d1113XZRKpfj+97//uu3PmDEjmpubO6dly5a9tUcHAEC3sV3HeH7mM5+JM844Y6vLjBgxIv70pz9FRNdjOuvq6mLEiBHx/PPPv+596+rqoq6ubnu6BADADmK7gmdjY2M0Nja+4XKjR4+Ourq6WLx4cRx22GEREdHW1hZLly6NoUOHvrmeAgCwQ0vyq/YBAwbEpz71qZg5c2Y0NTXF0KFD47LLLouIiFNPPTVFSQAAurlk5/G87LLLoqamJj760Y/Gxo0bY9y4cfGTn/wkBg4cmKokAADdWKkoiqLSnXg9a9asiYaGhhgfJ0VNqVelu0Oph1xhtSinr5FjXaV+HD1ke1f1Tn/ceNHWnrxGqVeycYKIyPMYio6O5DW8vrdRT3gMdBvtRVvMjzuiubk5BgwYsNVle8Y7CwAA3Z7gCQBAFoInAABZCJ4AAGSR9mh1epYecqB4qSb9D9WK9rbkNXrKj39SK2/cWOkuvC2yPKdSy/Gc9brYNn6ERYV4hQIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGRRU+kOQG5Fe1ulu/D2KMpp2y9l+Fya+jHQvdje7yy2N1tgxBMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIQvAEACCLmkp3AOiminKle8BfKiUeJ7C9t13qbZFDju2dYz153u5wesCrBwCAHYHgCQBAFoInAABZCJ4AAGQheAIAkIXgCQBAFoInAABZCJ4AAGSRLHg+88wzcdJJJ8Uuu+wSAwYMiMMOOyzmzZuXqhwAAN1csuB5/PHHR3t7e/zkJz+Jxx9/PA444IA4/vjj48UXX0xVEgCAbixJ8Fy1alX87ne/i+nTp8f73ve+eNe73hWXXnppbNiwIZ566qkUJQEA6OaSBM9BgwbFvvvuG//+7/8e69evj/b29rjmmmti8ODBMXr06BQlAQDo5mpSNFoqleL++++Pk08+Ofr37x9VVVUxePDguOeee2LgwIGve7+WlpZoaWnpvL1mzZoU3QMAoAK2a8Rz+vTpUSqVtjotWrQoiqKIqVOnxuDBg+PnP/95PProo3HyySfHCSecEH/6059et/1Zs2ZFQ0ND59TU1PSWHyAAAN1DqSiKYlsXXrlyZbz88stbXWbEiBHx85//PI466qh45ZVXYsCAAZ1/e9e73hVnnnlmTJ8+fYv33dKIZ1NTU4yPk6Km1GtbuwnQ85QSn/2uKKdtvydJvS1yyLG9c6wnz9tuob1oi/lxRzQ3N3fJfVuyXV+1NzY2RmNj4xsut2HDhoiIqKrq+qSrqqqKcvn1nyR1dXVRV1e3PV0CAGAHkeTjyMEHHxwDBw6MSZMmxZNPPhnPPPNMfO5zn4vnnnsujjvuuBQlAQDo5pIEz1122SXuueeeWLduXRxxxBExZsyY+MUvfhF33HFHHHDAASlKAgDQzSX5VXtExJgxY+Lee+9N1TwAADuYHnCENAAAOwLBEwCALARPAACySHaM59vhtVOMtkdbxDafbRSgJ3Iez+6jB4zZZNnezuP5TtEebRHxP7lta7p18Fy7dm1ERPwiflzhngBUmA/f3YdtsW2sp3ectWvXRkNDw1aX2a4rF+VWLpfjhRdeiP79+0epVKp0dzbz2pWVli1b9oZn6mf7Wb9pWb9pWb9pWb9pWb9p9bT1WxRFrF27NvbYY4/NLh7017r1iGdVVVUMGTKk0t14QwMGDOgRT5zuyvpNy/pNy/pNy/pNy/pNqyet3zca6XxNDzhQBQCAHYHgCQBAFoLnW1BXVxczZ86Murq6SnelR7J+07J+07J+07J+07J+03onr99u/eMiAAB6DiOeAABkIXgCAJCF4AkAQBaCJwAAWQieb5MTTzwx9tprr+jdu3fsvvvu8dGPfjReeOGFSnerR1i6dGmceeaZMXz48Kivr4+99947Zs6cGa2trZXuWo9xySWXxCGHHBJ9+vSJnXbaqdLd6RG+9a1vxbBhw6J3794xbty4ePTRRyvdpR7hZz/7WZxwwgmxxx57RKlUih/96EeV7lKPMmvWrBg7dmz0798/Bg8eHCeffHIsXry40t3qMebMmRPve9/7Ok8cf/DBB8fdd99d6W5lJXi+TSZMmBC33nprLF68OG677bb4/e9/H6ecckqlu9UjLFq0KMrlclxzzTXxm9/8Jv7t3/4trr766rjwwgsr3bUeo7W1NU499dT49Kc/Xemu9Ai33HJLXHDBBTFz5sxYuHBhHHDAAXH00UfHihUrKt21Hd769evjgAMOiG9961uV7kqP9NOf/jSmTp0aDz/8cMydOzfa2triqKOOivXr11e6az3CkCFD4tJLL43HH388HnvssTjiiCPipJNOit/85jeV7lo2TqeUyJ133hknn3xytLS0RK9evSrdnR7nsssuizlz5sSzzz5b6a70KNdff32cd9558eqrr1a6Kzu0cePGxdixY+PKK6+MiIhyuRxNTU1x9tlnx/Tp0yvcu56jVCrF7bffHieffHKlu9JjrVy5MgYPHhw//elP4/DDD690d3qknXfeOS677LI488wzK92VLIx4JrB69eq48cYb45BDDhE6E2lubo6dd9650t2AzbS2tsbjjz8eEydO7JxXVVUVEydOjIceeqiCPYPt19zcHBFhf5tAR0dH3HzzzbF+/fo4+OCDK92dbATPt9HnP//56Nu3bwwaNCief/75uOOOOyrdpR5pyZIlMXv27PjkJz9Z6a7AZlatWhUdHR2x6667dpm/6667xosvvlihXsH2K5fLcd5558Whhx4a++23X6W702P8+te/jn79+kVdXV186lOfittvvz1GjhxZ6W5lI3huxfTp06NUKm11WrRoUefyn/vc5+JXv/pV3HfffVFdXR0f+9jHwpEMr297129ExPLly+OYY46JU089NaZMmVKhnu8Y3sz6BXjN1KlT46mnnoqbb7650l3pUfbdd9944okn4pFHHolPf/rTMWnSpHj66acr3a1sHOO5FStXroyXX355q8uMGDEiamtrN5v/xz/+MZqamuLBBx98Rw2hb4/tXb8vvPBCjB8/Pt7//vfH9ddfH1VVPjdtzZt5/jrG861rbW2NPn36xA9+8IMuxx5OmjQpXn31Vd+EvI0c45nOtGnT4o477oif/exnMXz48Ep3p0ebOHFi7L333nHNNddUuitZ1FS6A91ZY2NjNDY2vqn7lsvliIhoaWl5O7vUo2zP+l2+fHlMmDAhRo8eHdddd53QuQ3eyvOXN6+2tjZGjx4dDzzwQGcgKpfL8cADD8S0adMq2zl4A0VRxNlnnx233357zJ8/X+jMoFwuv6OyguD5NnjkkUdiwYIFcdhhh8XAgQPj97//fXzhC1+Ivffe22jn22D58uUxfvz4GDp0aFx++eWxcuXKzr/ttttuFexZz/H888/H6tWr4/nnn4+Ojo544oknIiJin332iX79+lW2czugCy64ICZNmhRjxoyJgw46KL7xjW/E+vXrY/LkyZXu2g5v3bp1sWTJks7bzz33XDzxxBOx8847x1577VXBnvUMU6dOjZtuuinuuOOO6N+/f+dxyQ0NDVFfX1/h3u34ZsyYEccee2zstddesXbt2rjpppti/vz5ce+991a6a/kUvGX//d//XUyYMKHYeeedi7q6umLYsGHFpz71qeKPf/xjpbvWI1x33XVFRGxx4u0xadKkLa7fefPmVbprO6zZs2cXe+21V1FbW1scdNBBxcMPP1zpLvUI8+bN2+JzddKkSZXuWo/wevva6667rtJd6xE+/vGPF0OHDi1qa2uLxsbG4sgjjyzuu+++SncrK8d4AgCQhQPlAADIQvAEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAshA8AQDIQvAEACALwRMAgCwETwAAsvj/AXkj2vArP/8bAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(8, 8))\n", "z_draws = pm.draw(vars=z, draws=10_000)\n", "ax.hist2d(x=z_draws[:, 0], y=z_draws[:, 1], bins=25)\n", "ax.set(title=\"Samples Histogram\");" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Enough with Random Variables, I want to see some (log)probabilities!" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Recall we have defined the following model above:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\n", " \\begin{array}{rcl}\n", " \\text{z} &\\sim & \\operatorname{N}(\\text{},~\\text{})\n", " \\end{array}\n", " $$" ], "text/plain": [ "z ~ N(, )" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "`pymc` is able to convert `RandomVariable`s to their respective probability functions. One simple way is to use {func}`~pymc.logp`, which takes as first input a RandomVariable, and as second input the value at which the logp is evaluated (we will discuss this in more detail later)." ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "z_value = pt.vector(name=\"z\")\n", "z_logp = pm.logp(rv=z, value=z_value)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "`z_logp` contains the PyTensor graph that represents the log-probability of the normal random variable `z`, evaluated at `z_value`." ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Check{sigma > 0} [id A] 'z_logprob'\n", " |Elemwise{sub,no_inplace} [id B]\n", " | |Elemwise{sub,no_inplace} [id C]\n", " | | |Elemwise{mul,no_inplace} [id D]\n", " | | | |InplaceDimShuffle{x} [id E]\n", " | | | | |TensorConstant{-0.5} [id F]\n", " | | | |Elemwise{pow,no_inplace} [id G]\n", " | | | |Elemwise{true_div,no_inplace} [id H]\n", " | | | | |Elemwise{sub,no_inplace} [id I]\n", " | | | | | |z [id J]\n", " | | | | | |TensorConstant{(2,) of 0} [id K]\n", " | | | | |TensorConstant{[1. 2.]} [id L]\n", " | | | |InplaceDimShuffle{x} [id M]\n", " | | | |TensorConstant{2} [id N]\n", " | | |InplaceDimShuffle{x} [id O]\n", " | | |Elemwise{log,no_inplace} [id P]\n", " | | |Elemwise{sqrt,no_inplace} [id Q]\n", " | | |TensorConstant{6.283185307179586} [id R]\n", " | |Elemwise{log,no_inplace} [id S]\n", " | |TensorConstant{[1. 2.]} [id L]\n", " |All [id T]\n", " |MakeVector{dtype='bool'} [id U]\n", " |All [id V]\n", " |Elemwise{gt,no_inplace} [id W]\n", " |TensorConstant{[1. 2.]} [id L]\n", " |InplaceDimShuffle{x} [id X]\n", " |TensorConstant{0} [id Y]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pytensor.dprint(z_logp)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ ":::{tip}\n", "There is also a handy `pymc` function to compute the log cumulative probability of a random variable {func}`~pymc.logcdf`." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Observe that, as explained at the beginning, there has been no computation yet. The actual computation is performed after compiling and passing the input. For illustration purposes alone, we will again use the handy {meth}`~pytensor.graph.basic.Variable.eval` method." ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-0.91893853, -1.61208571])" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "z_logp.eval({z_value: [0, 0]})" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "This is nothing else than evaluating the log probability of a normal distribution." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-0.91893853, -1.61208571])" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scipy.stats.norm.logpdf(x=np.array([0, 0]), loc=np.array([0, 0]), scale=np.array([1, 2]))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "`pymc` models provide some helpful routines to facilitating the conversion of `RandomVariable`s to probability functions. {meth}`~pymc.Model.logp`, for instance can be used to extract the joint probability of all variables in the model:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Check{sigma > 0} [id A] 'z_logprob'\n", " |Elemwise{sub,no_inplace} [id B]\n", " | |Elemwise{sub,no_inplace} [id C]\n", " | | |Elemwise{mul,no_inplace} [id D]\n", " | | | |InplaceDimShuffle{x} [id E]\n", " | | | | |TensorConstant{-0.5} [id F]\n", " | | | |Elemwise{pow,no_inplace} [id G]\n", " | | | |Elemwise{true_div,no_inplace} [id H]\n", " | | | | |Elemwise{sub,no_inplace} [id I]\n", " | | | | | |z [id J]\n", " | | | | | |TensorConstant{(2,) of 0} [id K]\n", " | | | | |TensorConstant{[1. 2.]} [id L]\n", " | | | |InplaceDimShuffle{x} [id M]\n", " | | | |TensorConstant{2} [id N]\n", " | | |InplaceDimShuffle{x} [id O]\n", " | | |Elemwise{log,no_inplace} [id P]\n", " | | |Elemwise{sqrt,no_inplace} [id Q]\n", " | | |TensorConstant{6.283185307179586} [id R]\n", " | |Elemwise{log,no_inplace} [id S]\n", " | |TensorConstant{[1. 2.]} [id L]\n", " |All [id T]\n", " |MakeVector{dtype='bool'} [id U]\n", " |All [id V]\n", " |Elemwise{gt,no_inplace} [id W]\n", " |TensorConstant{[1. 2.]} [id L]\n", " |InplaceDimShuffle{x} [id X]\n", " |TensorConstant{0} [id Y]\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pytensor.dprint(model.logp(sum=False))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Because we only have one variable, this function is equivalent to what we obtained by manually calling `pm.logp` before. We can also use a helper {meth}`~pymc.Model.compile_logp` to return an already compiled PyTensor function of the model logp." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "logp_function = model.compile_logp(sum=False)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "This function expects a \"point\" dictionary as input. We could create it ourselves, but just to illustrate another useful {class}`~pymc.Model` method, let's call {meth}`~pymc.Model.initial_point`, which returns the point that most samplers use when deciding where to start sampling." ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'z': array([0., 0.])}" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "point = model.initial_point()\n", "point" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[array([-0.91893853, -1.61208571])]" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "logp_function(point)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### What are value variables and why are they important?\n", "\n", "As we have seen above, a logp graph does not have random variables. Instead it's defined in terms of input (value) variables. When we want to sample, each random variable (RV) is replaced by a logp function evaluated at the respective input (value) variable. Let's see how this works through some examples. RV and value variables can be observed in these [`scipy`](https://github.com/scipy/scipy) operations:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rv = scipy.stats.norm(0, 1)\n", "\n", "# Equivalent to rv = pm.Normal(\"rv\", 0, 1)\n", "scipy.stats.norm(0, 1)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-1.41907251, -1.01111034, -0.16152042])" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Equivalent to rv_draw = pm.draw(rv, 3)\n", "rv.rvs(3)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-1.7001885332046727" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Equivalent to rv_logp = pm.logp(rv, 1.25)\n", "rv.logpdf(1.25)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's look at how these value variables behave in a slightly more complex model." ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "with pm.Model() as model_2:\n", " mu = pm.Normal(name=\"mu\", mu=0, sigma=2)\n", " sigma = pm.HalfNormal(name=\"sigma\", sigma=3)\n", " x = pm.Normal(name=\"x\", mu=mu, sigma=sigma)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Each model RV is related to a \"value variable\":" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{mu ~ N(0, 2): mu, sigma ~ N**+(0, 3): sigma_log__, x ~ N(mu, sigma): x}" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_2.rvs_to_values" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Observe that for sigma the associated value is in the *log* scale as in practice we require unbounded values for NUTS sampling." ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[mu, sigma_log__, x]" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_2.value_vars" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now that we know how to extract the model variables, we can compute the element-wise log-probability of the model for specific values." ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ -1.61208571, -11.32440364, 9.08106147])" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# extract values as pytensor.tensor.var.TensorVariable\n", "mu_value = model_2.rvs_to_values[mu]\n", "sigma_log_value = model_2.rvs_to_values[sigma]\n", "x_value = model_2.rvs_to_values[x]\n", "# element-wise log-probability of the model (we do not take te sum)\n", "logp_graph = pt.stack(model_2.logp(sum=False))\n", "# evaluate by passing concrete values\n", "logp_graph.eval({mu_value: 0, sigma_log_value: -10, x_value: 0})" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "This equivalent to:" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "mu_value -> -1.612085713764618\n", "sigma_log_value -> -11.324403641427345 \n", "x_value -> 9.081061466795328\n", "\n" ] } ], "source": [ "print(\n", " f\"\"\"\n", "mu_value -> {scipy.stats.norm.logpdf(x=0, loc=0, scale=2)}\n", "sigma_log_value -> {- 10 + scipy.stats.halfnorm.logpdf(x=np.exp(-10), loc=0, scale=3)} \n", "x_value -> {scipy.stats.norm.logpdf(x=0, loc=0, scale=np.exp(-10))}\n", "\"\"\"\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ ":::{Note}\n", "For `sigma_log_value` we add the $-10$ term for the `scipy` and `pytensor` to match because of the jacobian.\n", ":::" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "As we already saw, we can also use the method {meth}`~pymc.Model.compile_logp` to obtain a compiled pytensor function of the model logp, which takes a dictionary of `{value variable name : value}` as inputs:" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[array(-1.61208571), array(-11.32440364), array(9.08106147)]" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model_2.compile_logp(sum=False)({\"mu\": 0, \"sigma_log__\": -10, \"x\": 0})" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The {class}`~pymc.Model` class also has methods to extract the gradient ({meth}`~pymc.Model.dlogp`) and the hessian ({meth}`~pymc.Model.d2logp`) of the logp." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "If you want to go deeper into the internals of `pytensor` RandomVariables and `pymc` distributions please take a look into the [distribution developer guide](implementing-a-distribution)." ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: Tue Dec 06 2022\n", "\n", "Python implementation: CPython\n", "Python version : 3.11.0\n", "IPython version : 8.7.0\n", "\n", "pytensor: 2.8.10\n", "\n", "numpy : 1.23.4\n", "scipy : 1.9.3\n", "pymc : 4.4.0+207.g7c3068a1c\n", "pytensor : 2.8.10\n", "matplotlib: 3.6.2\n", "\n", "Watermark: 2.3.1\n", "\n" ] } ], "source": [ "%load_ext watermark\n", "%watermark -n -u -v -iv -w -p pytensor" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "colab": { "collapsed_sections": [ "dIYxBNT_5HgW", "JhmIBByY6T9h", "k7spOYvvTdZL", "qIMN2gHGzKof", "aMM0sJy6z7Ur", "WY6bUgqZztC3", "58gng2f17_P7", "C8Us4nEyhRdu", "wkZR0gDWRAgK", "wPw9kCvASOeJ", "WdZcUfvLUkwK", "EHH1hP0uzaWG", "I_Ph4o7ZW_XD", "jOI-E76fZ2bG", "Rq2W2fmobmFg" ], "name": "PyMC intro.ipynb", "provenance": [] }, "hide_input": false, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.16" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "vscode": { "interpreter": { "hash": "867ba48c05011db76db56a12fb95ccd32f7ac276df8f4ae698e0d475911a6ba0" } } }, "nbformat": 4, "nbformat_minor": 1 }