{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "(pymc_aesara)=\n", "\n", "# PyMC and Aesara\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 Aesara graphs. The purpose is not to give a detailed description of all [`aesara`](https://github.com/aesara-devs/aesara)'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." ] }, { "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 aesara\n", "import aesara.tensor as at\n", "import pymc as pm\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy.stats" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction to Aesara\n", "\n", "We start by looking into `aesara`. According to their documentation\n", "\n", "> Aesara is a Python library that allows one to define, optimize, and efficiently evaluate mathematical expressions involving multi-dimensional arrays." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![aesara logo](https://raw.githubusercontent.com/aesara-devs/aesara/main/doc/images/aesara_logo_2400.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A simple example\n", "\n", "To begin, we define some aesara 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, (None,))\n", "y name = y\n", "\n" ] } ], "source": [ "x = at.scalar(name=\"x\")\n", "y = at.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", ")" ] }, { "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\"" ] }, { "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 = at.log(z)\n", "w.name = \"log(x + y)\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use the {func}`~aesara.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": [ "aesara.dprint(w)" ] }, { "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}`~aesara.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 = aesara.function(inputs=[x, y], outputs=w)" ] }, { "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])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{tip}\n", "Sometimes we just want to debug, we can use {meth}`~aesara.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]})" ] }, { "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]})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Aesara is clever!\n", "\n", "One of the most important features of `aesara` 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 = at.scalar(name=\"a\")\n", "b = at.scalar(name=\"b\")\n", "\n", "c = a / b\n", "c.name = \"a / b\"\n", "\n", "aesara.dprint(c)" ] }, { "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", "aesara.dprint(d)" ] }, { "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 = aesara.function(inputs=[a, b], outputs=d)\n", "\n", "aesara.dprint(g)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What is in an Aesara graph?\n", "\n", "The following diagram shows the basic structure of an `aesara` graph.\n", "\n", "![aesara graph](https://raw.githubusercontent.com/aesara-devs/aesara/main/doc/tutorial/apply.png)" ] }, { "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, (None,))\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", ")" ] }, { "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, (None,))\n", " > Op is Elemwise{log,no_inplace}\n", " > Input 0 is x + y\n", "---\n", "Checking variable x + y of type TensorType(float64, (None,))\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, (None,))\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\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this is very similar to the output of {func}`~aesara.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": [ "aesara.dprint(w)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Graph manipulation 101\n", "\n", "Another interesting feature of Aesara 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(aesara.graph.graph_inputs(graphs=[w]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a simple example, let's add an {func}`~aesara.tensor.exp` before the {func}`~aesara.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 = at.exp(parent_of_w) # modify the parent of w\n", "new_parent_of_w.name = \"exp(x + y)\"" ] }, { "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": [ "aesara.dprint(w)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To modify the graph we need to use the {func}`~aesara.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 = aesara.clone_replace(output=[w], replace={parent_of_w: new_parent_of_w})[0]\n", "new_w.name = \"log(exp(x + y))\"\n", "aesara.dprint(new_w)" ] }, { "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]})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the new graph is just the identity function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{note}\n", "Again, note that `aesara` 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 = aesara.function(inputs=[x, y], outputs=new_w)\n", "\n", "aesara.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])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Aesara RandomVariables\n", "\n", "Now that we have seen aesara'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": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAF1CAYAAAAeOhj3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgAElEQVR4nO3deZhkdX3v8feHRVxRCA0CA44LEhFNTOYSk3gNEb0SMWDyxFyMmnFJiIkavTGGQZMg6iQYs+jVeBOibC4Q4hJRNAFRJMYIDorKojLKCAMjM4CERUUHv/ePc1qKpnumpmeqqvs379fz9FN1ljrne35V1Z865/zqVKoKSZK0uO0w6QIkSdLWM9AlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiaiCSvS/KebbzMJDklyXeSXLwtl709SfKCJJ/ZgvnXJHlqf/81Sd65DWu5Pckj+vunJnnjNlz2PyT5s221vK21rdtO25+dJl2AxivJk4C/Ah4L3AVcCbyyqj4/0cK2jScBTwOWVNUdky5me1RVfzHMfEkuAN5TVZsMsKp64LaoK8kLgN+pqicNLPsl22LZ28qwbSfNxUDfjiTZFfgo8PvAWcB9gP8J3DnJurahhwFr5grzJDtV1cYx1zRyLW5Xi9skjZqH3LcvjwaoqjOq6q6q+l5VnVtVXwZI8sgkn0xyU5Ibk7w3yUOmH9wfWn11ki8nuSPJu5LsleTjSW5L8okku/XzLk1SSY5Jcn2SdUleNVdhSZ6Y5LNJbknypSSHDkx7QZJv9uu4OslzZ3n8i4F3Aj/fH6Y9IcmhSdYmOTbJt4FTkuyS5C19Tdf393fplzE9/58kWd/X/Kwkz0jy9SQ3J3nNJrbhiCRfTHJrkmuTvG4T806v61UD63rhwPQHJzk9yYYk30ryp0l2GGiP/0zyd0luBl7XH45+R/9c3N5Pf2i/fd9J8tUkTxhY/ook3+jb9IokvzZXrbPU/vy+ppuSvHbGtB+fSkly3yTv6ee7Jcnn+9fLSroPkm/va317P38leWmSq4CrBsY9amAVeyQ5r6/700ke1s83/XrbaaCWC5L8TpLHAP/A3a+NW/rp9ziEn+R3k6zun+ezk+wzMK2SvCTJVX17/n2SzNE+M5d7aJK1A8PHJrmu34avJTlslrab3p7lSa5J93587cAy7pfktL6WK/vX7FrmsKn6M+P018y27Nvxjenen7cn+UiSn0j3/+HW/nldOmNdf5juPXtjkjcn2SHde+/mJI8bmHfPJN9LMjVX7doCVeXfdvIH7ArcBJwG/Aqw24zpj6I7ZL0LMAVcCLxlYPoa4HPAXsC+wHrgC8AT+sd8Eji+n3cpUMAZwAOAxwEbgKf2019Hd8iVflk3Ac+g+5D5tH54qn/srcCB/bx7A4+dY/teAHxmYPhQYCPwpr6++wGv77dhz375nwXeMGP+Pwd2Bn63r/l9wIPoTlN8H3jEHOs/tN/OHYDHAzcAz9rEvBv7enbut/27088JcDrw4X69S4GvAy8e2M6NwMvpjrLdDzgVuBH4WeC+/XNxNfDbwI7AG4FPDaz/2cA+fa3/G7gD2Hu2dpxR90HA7cCT+zb9276W2Z7X3wM+Aty/r+FngV37aRfQHQIfXHYB5wG7A/cbGPeo/v6pwG0D637rdJ3c/XrbaWB5P17HbNvUL++N/f2n9O33M/2y3wZcOKO2jwIPAfane10cPkcb/Xi5A8/12v7+gcC1wD4DdT9ylrab3p5/6p/fn6I7kvaYfvqJwKeB3YAlwJen1zFHTXPWP7je2dqyb8fVwCOBBwNX0L0en0r3+jsdOGXGuj7VP4/79/NOPw/vAN40MO8rgI9M+n9jK3/uoW9HqupWuvPM0/8oNvR7Inv101dX1XlVdWdVbaD7Z/1LMxbztqq6oaquA/4DuKiqvlhVdwIfogv3QSdU1R1V9RXgFOA5s5T2POBjVfWxqvpRVZ0HrKILOYAfAQcnuV9Vrauqy7dgs39E9yHjzqr6HvBc4PVVtb7fxhOA5w/M/0NgZVX9EDgT2AN4a1Xd1q/3crqwvpequqCqvtJvw5fpPszMbL9BP+xr+WFVfYwuKA9MsiNdyB7Xr3cN8Dcz6ry+qt5WVRv77QL4UFVdUlXfp3suvl9Vp1fVXcA/M/DcVNW/VNX1fa3/TLdHfMhm2hLgN4CPVtWF/XP+Z3RtPNf2/QRdIN/V13brZpb/l1V188A2zXTOwLpfS7fXvd8QdW/Oc4GTq+oL/bKP65e9dGCeE6vqlqq6hi6wfnoe67mL7gPDQUl2rqo1VfWNTcx/QnVH0r4EfIku2AF+E/iLqvpOVa0F/u8Q696a+k+pqm9U1X8DHwe+UVWfqO60yL9w7/f9m/rn8RrgLdz9vj8N+K30R5voXtPv3oI6tAkG+namqq6sqhdU1RLgYLq9tLfAjw9/ndkfDrwVeA9doA26YeD+92YZntmJ6dqB+9/q1zfTw4Bn94dlb+kPiT6Jbo/xDrpwewmwLsk5SX5y+C1mQx9w0/bp65irppv6AJzeHtj8NgKQ5OeSfCrdYfL/7mue2X6Dbqp7nif+br/sPej6N8ysc9+B4cF2nTb0c5Pkt5NcOtDeB2+m1mn7DK67f35ummPedwP/DpyZ7vTGXyXZeTPLn227Zp1eVbcDNzP7a2pL3eN10S/7Ju7Z5t8euD/9XG2RqloNvJJur3h9/37bVP1zrfMezwObb7dNLWsY2+R9X1UX0R0N+qX+ffwo4OwtqEObYKBvx6rqq3SHBw/uR/0l3d7746tqV7o951nPE26Bwb2n/YHrZ5nnWuDdVfWQgb8HVNWJfZ3/XlVPozvc/lW6owvDmvlzgtfTfYDYXE3z8T66f077VdWD6c7bzqf9bqTbu51Z53UDw/P+mcT+vPM/AS8DfqKqHgJcxnC1rmPgOU1yf7q98HvpjzycUFUHAb8APJPuFMCm6t/cdg2u+4F0h3WvpwsJ6A7vT3voFiz3Hq+LJA+g267r5nzE3O7YRB1U1fuq623/sL6uN81jHevoDrVP25qjFJusd5429b4/je5/y/OB98/4wK2tYKBvR5L8ZLpOWEv64f3oDoV9rp/lQXSHfW9Jsi/w6m2w2j9Lcv8kjwVeSHfod6b3AL+a5OlJdkzXmerQJEvSdaI6sv8He2df312zLGNYZwB/mmQqyR5058u31ffhHwTcXFXfT3II8FvzWUh/hOAsYGWSB/UB/EfbsM4H0AXJBoB0nfEO3uQj7vZ+4JlJnpTkPnR9AGb9P5Lkl5M8rj+FcCvdh5Tp5+4G4BHzqP0ZA+t+A90pn2v70yfXAc/rX0MvojvnO+0GYEn/uNm8D3hhkp9O10nyL/plr5lHjZf2de6e5KF0e+QAJDkwyVP6dXyfbu92Pq/ns4DjkuzWv1dfNo9lDNb75CT7J3kw3emGrfXqvrb96M6TD77v3w38Gl2on74N1qWegb59uQ34OeCiJHfQBfllwHTv8xPoOgX9N3AO8MFtsM5P03WoOR/466o6d+YMVXUtcBTwGrqQuZbuw8QO/d+r6D7h30x3TvoPtqKeN9Kdn/8y8BW6Tn3b6mIlfwC8PsltdB8UztqKZb2cbs/pm8Bn6ALn5K2uEKiqK+jOyf8XXdA9DvjPIR97OfDSvp51wHeAuXpXP5TuA8CtdNc7+DR3fyh5K/AbfY/rYc7/TnsfcDzda+Fn6c59T/tdutfNTXQdGD87MO2TdP0fvp3kxlm263y6/gAf6LfrkcDRW1DXoHfTne9eA5zLPcNsF7oObTfSHQLfk+51v6VeT9fuVwOfoGvneX39tO+z8s9074lL6DrPba0P98u6lO5/ybsG1reW7n1XdP1wtI2kat5H7qQ59Z2JrgZ2Lr9PLI1Ukt8Hjq6qTXXCHFctBRzQ9xeYa56T6Tp2/un4KmufF5aRpEUmyd50pyz+CziA7ijW2yda1JD6D/u/zr17xmsrechdkhaf+wD/SHca7ZN0h7jfMdGKhpDkDXSn+d5cVVdPup7WeMhdkqQGuIcuSVIDDHRJkhqwqDvF7bHHHrV06dJJlyFJ0thccsklN1bVvX7QZlEH+tKlS1m1atWky5AkaWySfGu28R5ylySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhqwqH9tTdLWW7rinLGta82JR4xtXdL2xj10SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGjCzQk5ycZH2Sy2aMf3mSryW5PMlfDYw/LsnqftrTR1WXJEktGuWPs5wKvB04fXpEkl8GjgIeX1V3JtmzH38QcDTwWGAf4BNJHl1Vd42wPkmSmjGyPfSquhC4ecbo3wdOrKo7+3nW9+OPAs6sqjur6mpgNXDIqGqTJKk14z6H/mjgfya5KMmnk/yPfvy+wLUD863tx91LkmOSrEqyasOGDSMuV5KkxWHcgb4TsBvwRODVwFlJAmSWeWu2BVTVSVW1rKqWTU1Nja5SSZIWkXEH+lrgg9W5GPgRsEc/fr+B+ZYA14+5NkmSFq1xB/q/Ak8BSPJo4D7AjcDZwNFJdknycOAA4OIx1yZJ0qI1sl7uSc4ADgX2SLIWOB44GTi5/yrbD4DlVVXA5UnOAq4ANgIvtYe7JEnDG1mgV9Vz5pj0vDnmXwmsHFU9kiS1zCvFSZLUgFFeWEaS7mHpinPGtq41Jx4xtnVJC4F76JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBuw06QIk3dvSFedMugRJi4x76JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDVgZIGe5OQk65NcNsu0P05SSfYYGHdcktVJvpbk6aOqS5KkFo1yD/1U4PCZI5PsBzwNuGZg3EHA0cBj+8e8I8mOI6xNkqSmjCzQq+pC4OZZJv0d8CdADYw7Cjizqu6sqquB1cAho6pNkqTWjPUcepIjgeuq6kszJu0LXDswvLYfN9syjkmyKsmqDRs2jKhSSZIWl7EFepL7A68F/ny2ybOMq1nGUVUnVdWyqlo2NTW1LUuUJGnRGufPpz4SeDjwpSQAS4AvJDmEbo98v4F5lwDXj7E2SZIWtbHtoVfVV6pqz6paWlVL6UL8Z6rq28DZwNFJdknycOAA4OJx1SZJ0mI3yq+tnQH8F3BgkrVJXjzXvFV1OXAWcAXwb8BLq+quUdUmSVJrRnbIvaqes5npS2cMrwRWjqoeSZJa5pXiJElqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhowskBPcnKS9UkuGxj35iRfTfLlJB9K8pCBacclWZ3ka0mePqq6JElq0Sj30E8FDp8x7jzg4Kp6PPB14DiAJAcBRwOP7R/zjiQ7jrA2SZKaMrJAr6oLgZtnjDu3qjb2g58DlvT3jwLOrKo7q+pqYDVwyKhqkySpNZM8h/4i4OP9/X2Bawemre3HSZKkIew0iZUmeS2wEXjv9KhZZqs5HnsMcAzA/vvvP5L6JC1+S1ecM7Z1rTnxiLGtS5rL2PfQkywHngk8t6qmQ3stsN/AbEuA62d7fFWdVFXLqmrZ1NTUaIuVJGmRGGugJzkcOBY4sqq+OzDpbODoJLskeThwAHDxOGuTJGkxG9kh9yRnAIcCeyRZCxxP16t9F+C8JACfq6qXVNXlSc4CrqA7FP/SqrprVLVJktSakQV6VT1nltHv2sT8K4GVo6pHkqSWeaU4SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEj+7U1qTVLV5wz6RIkaU7uoUuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWrAyAI9yclJ1ie5bGDc7knOS3JVf7vbwLTjkqxO8rUkTx9VXZIktWiUe+inAofPGLcCOL+qDgDO74dJchBwNPDY/jHvSLLjCGuTJKkpIwv0qroQuHnG6KOA0/r7pwHPGhh/ZlXdWVVXA6uBQ0ZVmyRJrRn3OfS9qmodQH+7Zz9+X+DagfnW9uMkSdIQFkqnuMwyrmadMTkmyaokqzZs2DDisiRJWhzGHeg3JNkboL9d349fC+w3MN8S4PrZFlBVJ1XVsqpaNjU1NdJiJUlaLMYd6GcDy/v7y4EPD4w/OskuSR4OHABcPObaJElatHYa1YKTnAEcCuyRZC1wPHAicFaSFwPXAM8GqKrLk5wFXAFsBF5aVXeNqjZJklozskCvqufMMemwOeZfCawcVT2SJLVsoXSKkyRJW8FAlySpAQa6JEkNMNAlSWqAgS5JUgOGCvQk5w8zTpIkTcYmv7aW5L7A/em+S74bd1+idVdgnxHXJkmShrS576H/HvBKuvC+hLsD/Vbg70dXliRJ2hKbDPSqeivw1iQvr6q3jakmSZK0hYa6UlxVvS3JLwBLBx9TVaePqC5JkrQFhgr0JO8GHglcCkxfY70AA12SpAVg2Gu5LwMOqqpZf6NckiRN1rDfQ78MeOgoC5EkSfM37B76HsAVSS4G7pweWVVHjqQqSZK0RYYN9NeNsghJkrR1hu3l/ulRFyJJkuZv2F7ut9H1age4D7AzcEdV7TqqwiRJ0vCG3UN/0OBwkmcBh4yiIEmStOXm9WtrVfWvwFO2bSmSJGm+hj3k/usDgzvQfS/d76RLkrRADNvL/VcH7m8E1gBHbfNqJEnSvAx7Dv2Foy5Emo+lK86ZdAmStCAMdQ49yZIkH0qyPskNST6QZMmoi5MkScMZtlPcKcDZdL+Lvi/wkX6cJElaAIYN9KmqOqWqNvZ/pwJTI6xLkiRtgWED/cYkz0uyY//3POCmURYmSZKGN2ygvwj4TeDbwDrgNwA7ykmStEAM+7W1NwDLq+o7AEl2B/6aLuglSdKEDbuH/vjpMAeoqpuBJ8x3pUn+T5LLk1yW5Iwk902ye5LzklzV3+423+VLkrS9GTbQdxgM2H4Pfdi9+3tIsi/wh8CyqjoY2BE4GlgBnF9VBwDn98OSJGkIw4by3wCfTfJ+uku+/iawcivXe78kPwTuD1wPHAcc2k8/DbgAOHYr1iFJ0nZj2CvFnZ5kFd0PsgT49aq6Yj4rrKrrkvw1cA3wPeDcqjo3yV5Vta6fZ12SPeezfEmStkdDHzbvA3xeIT6oP3R/FPBw4BbgX/qvwQ37+GOAYwD233//rS1HkqQmzOvnU7fSU4Grq2pDVf0Q+CDwC8ANSfYG6G/Xz/bgqjqpqpZV1bKpKa9tI0kSTCbQrwGemOT+SQIcBlxJd2nZ5f08y4EPT6A2SZIWpXn1VN8aVXVR37nuC3Q/xfpF4CTggcBZSV5MF/rPHndtkiQtVmMPdICqOh44fsboO+n21iVJ0haaxCF3SZK0jRnokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ2YSKAneUiS9yf5apIrk/x8kt2TnJfkqv52t0nUJknSYjSpPfS3Av9WVT8J/BRwJbACOL+qDgDO74clSdIQdhr3CpPsCjwZeAFAVf0A+EGSo4BD+9lOAy4Ajh13fZK0pZauOGds61pz4hFjW5cWl0nsoT8C2ACckuSLSd6Z5AHAXlW1DqC/3XMCtUmStChNItB3An4G+H9V9QTgDrbg8HqSY5KsSrJqw4YNo6pRkqRFZRKBvhZYW1UX9cPvpwv4G5LsDdDfrp/twVV1UlUtq6plU1NTYylYkqSFbuyBXlXfBq5NcmA/6jDgCuBsYHk/bjnw4XHXJknSYjX2TnG9lwPvTXIf4JvAC+k+XJyV5MXANcCzJ1SbJEmLzkQCvaouBZbNMumwMZciSVITvFKcJEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktSAnSZdgNqzdMU5ky5BkrY77qFLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWrAxAI9yY5Jvpjko/3w7knOS3JVf7vbpGqTJGmxmeQe+iuAKweGVwDnV9UBwPn9sCRJGsJELiyTZAlwBLAS+KN+9FHAof3904ALgGPHXZskLWTjvHDTmhOPGNu6tPUmtYf+FuBPgB8NjNurqtYB9Ld7TqAuSZIWpbEHepJnAuur6pJ5Pv6YJKuSrNqwYcM2rk6SpMVpEnvovwgcmWQNcCbwlCTvAW5IsjdAf7t+tgdX1UlVtayqlk1NTY2rZkmSFrSxB3pVHVdVS6pqKXA08Mmqeh5wNrC8n2058OFx1yZJ0mK1kL6HfiLwtCRXAU/rhyVJ0hAm+vOpVXUBXW92quom4LBJ1iNJ0mK1kPbQJUnSPBnokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqwE6TLkDjsXTFOZMuQZI0Qu6hS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqwNgDPcl+ST6V5Moklyd5RT9+9yTnJbmqv91t3LVJkrRYTWIPfSPwqqp6DPBE4KVJDgJWAOdX1QHA+f2wJEkawtgDvarWVdUX+vu3AVcC+wJHAaf1s50GPGvctUmStFhN9Bx6kqXAE4CLgL2qah10oQ/sOcHSJElaVCYW6EkeCHwAeGVV3boFjzsmyaokqzZs2DC6AiVJWkQmEuhJdqYL8/dW1Qf70Tck2bufvjewfrbHVtVJVbWsqpZNTU2Np2BJkha4SfRyD/Au4Mqq+tuBSWcDy/v7y4EPj7s2SZIWq0n82tovAs8HvpLk0n7ca4ATgbOSvBi4Bnj2BGqTJGlRGnugV9VngMwx+bBx1iJJUiu8UpwkSQ0w0CVJasAkzqFLkhaBpSvOGdu61px4xNjW1Sr30CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAn0+doHH+NKEkLWT+VOvWcw9dkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgL3cZ7DnuSRpMXIPXZKkBhjokiQ1wECXJKkBBrokSQ1YcJ3ikhwOvBXYEXhnVZ044ZIkSQ1p9TKzC2oPPcmOwN8DvwIcBDwnyUGTrUqSpIVvQQU6cAiwuqq+WVU/AM4EjppwTZIkLXgLLdD3Ba4dGF7bj5MkSZuw0M6hZ5ZxdY8ZkmOAY/rB25N8beRVwR7AjWNYz0JmG9gGYBuAbQC2AQzZBnnTSNb9sNlGLrRAXwvsNzC8BLh+cIaqOgk4aZxFJVlVVcvGuc6FxjawDcA2ANsAbANYmG2w0A65fx44IMnDk9wHOBo4e8I1SZK04C2oPfSq2pjkZcC/031t7eSqunzCZUmStOAtqEAHqKqPAR+bdB0zjPUQ/wJlG9gGYBuAbQC2ASzANkhVbX4uSZK0oC20c+iSJGkeDPQhJHlDki8nuTTJuUn2mXRN45bkzUm+2rfDh5I8ZNI1jVuSZye5PMmPkiyo3q2jluTwJF9LsjrJiknXMwlJTk6yPsllk65lUpLsl+RTSa7s3wuvmHRN45bkvkkuTvKlvg1OmHRN0zzkPoQku1bVrf39PwQOqqqXTLissUryv4BP9h0X3wRQVcdOuKyxSvIY4EfAPwJ/XFWrJlzSWPSXZP468DS6r5Z+HnhOVV0x0cLGLMmTgduB06vq4EnXMwlJ9gb2rqovJHkQcAnwrO3ptZAkwAOq6vYkOwOfAV5RVZ+bcGnuoQ9jOsx7D2DGxW62B1V1blVt7Ac/R3eNgO1KVV1ZVeO4kNFC4yWZgaq6ELh50nVMUlWtq6ov9PdvA65kO7uaZ3Vu7wd37v8WRCYY6ENKsjLJtcBzgT+fdD0T9iLg45MuQmPjJZl1L0mWAk8ALppwKWOXZMcklwLrgfOqakG0gYHeS/KJJJfN8ncUQFW9tqr2A94LvGyy1Y7G5tqgn+e1wEa6dmjOMG2wHdrsJZm1fUnyQOADwCtnHMHcLlTVXVX103RHKg9JsiBOwSy476FPSlU9dchZ3wecAxw/wnImYnNtkGQ58EzgsGq088UWvA62J5u9JLO2H/154w8A762qD066nkmqqluSXAAcDky8s6R76ENIcsDA4JHAVydVy6QkORw4Fjiyqr476Xo0Vl6SWcCPO4S9C7iyqv520vVMQpKp6W/5JLkf8FQWSCbYy30IST4AHEjXw/lbwEuq6rrJVjVeSVYDuwA39aM+tx329P814G3AFHALcGlVPX2iRY1JkmcAb+HuSzKvnGxF45fkDOBQul/ZugE4vqreNdGixizJk4D/AL5C9/8Q4DX9FT63C0keD5xG917YATirql4/2ao6BrokSQ3wkLskSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAf8fx1ZFpDaNFggAAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" }, "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\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's try to do it in Aesara." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": "TensorType(float64, ())" }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = at.random.normal(loc=0, scale=1, name=\"y\")\n", "y.type" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we show the graph using {func}`~aesara.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": [ "aesara.dprint(y)" ] }, { "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`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We *could* sample by calling {meth}`~aesara.graph.basic.Variable.eval`. on the random variable." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": "array(-0.8043385)" }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.eval()" ] }, { "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: -0.804338501335673\n", "Sample 1: -0.804338501335673\n", "Sample 2: -0.804338501335673\n", "Sample 3: -0.804338501335673\n", "Sample 4: -0.804338501335673\n", "Sample 5: -0.804338501335673\n", "Sample 6: -0.804338501335673\n", "Sample 7: -0.804338501335673\n", "Sample 8: -0.804338501335673\n", "Sample 9: -0.804338501335673\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {y.eval()}\")" ] }, { "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## PyMC\n", "\n", "![pymc logo](https://raw.githubusercontent.com/pymc-devs/pymc/main/docs/logos/PyMC.png)" ] }, { "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", "aesara.dprint(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that `x` is just a normal `RandomVariable` and which is the same as `y` except for the `rng`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can try to generate samples by calling {meth}`~aesara.graph.basic.Variable.eval` as above." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sample 0: 0.9562139156242165\n", "Sample 1: 0.9562139156242165\n", "Sample 2: 0.9562139156242165\n", "Sample 3: 0.9562139156242165\n", "Sample 4: 0.9562139156242165\n", "Sample 5: 0.9562139156242165\n", "Sample 6: 0.9562139156242165\n", "Sample 7: 0.9562139156242165\n", "Sample 8: 0.9562139156242165\n", "Sample 9: 0.9562139156242165\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {x.eval()}\")" ] }, { "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": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAF1CAYAAAAeOhj3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAefklEQVR4nO3deZhkdX3v8feHRVRkDaOyOoCEBNHAzTzk3kQNcYnEa9yeaCBqcIlIIok+MbmuCYNIonG/GjUoCKigRMLVRHIDQYUYgzogIquCjDIwDsOMrCrXGb73j3MairZ7pmemq6r71+/X89TTVb+zfc+pqv7U+Z1Tp1JVSJKk+W2rcRcgSZK2nIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wEDXnJBkaZJPzvI8k+TjSX6U5OuzOe+FJMlLk3xlE8ZfnuRp/f03JfnYLNZyd5L9+vunJXnbLM77I0n+arbmt6Vme9upfduMuwCNV5InAn8HPA5YD1wDvLaqvjHWwmbHE4GnA3tV1T3jLmYhqqq/mcl4Sb4MfLKqNhhgVfWI2agryUuBP6qqJw7M+9jZmPdsmem2kyYY6AtYkh2BfwH+GDgbeAjwJODecdY1ix4DLJ8uzJNsU1XrRlzT0LW4Xi2ukzTb7HJf2H4RoKrOqqr1VfWTqjq/qq4ASLJ/ki8mWZPktiSfSrLzxMR91+pfJrkiyT1JTknyqCT/muSuJP+eZJd+3MVJKskxSW5JsjLJ66YrLMl/T/LVJLcn+VaSwweGvTTJ9/pl3JjkRVNM/wrgY8D/6LtpT0hyeJIVSV6f5IfAx5Nsl+R9fU239Pe36+cxMf7/SnJrX/NzkzwzyXeSrE3ypg2sw/9M8s0kdya5KcnSDYw7sazXDSzrZQPDd0pyRpLVSb6f5C1JthrYHv+Z5L1J1gJL++7oD/XPxd398Ef36/ejJNcmOXRg/m9IckO/Ta9O8rzpap2i9pf0Na1J8uZJw+4/lJLkoUk+2Y93e5Jv9K+Xk+g+SH6wr/WD/fiV5NVJvgt8d6DtsQOL2C3JBX3dFyV5TD/exOttm4Favpzkj5L8MvARHnht3N4Pf1AXfpJXJrm+f54/n2SPgWGV5Ngk3+23598nyTTbZ/J8D0+yYuDx65Pc3K/DdUmeOsW2m1ifo5P8IN378c0D83hYktP7Wq7pX7MrmEY/rz9L9z66Lck7k2yV7v2wNsnjB8Z9ZJKfJFmUTXxPJNk63aGDidfWpUn2nq4ubaGq8rZAb8COwBrgdOB3gF0mDX8sXZf1dsAi4GLgfQPDlwOXAI8C9gRuBS4DDu2n+SJwfD/uYqCAs4DtgccDq4Gn9cOX0nW50s9rDfBMug+dT+8fL+qnvRM4sB93d+Bx06zfS4GvDDw+HFgHvKOv72HAW/t1eGQ//68CJ04a/6+BbYFX9jWfCexAd5jip8B+0yz/8H49twKeAKwCnruBcdf19Wzbr/uPJ54T4Azgc/1yFwPfAV4xsJ7rgD+l63V7GHAacBvwq8BD++fiRuAPga2BtwFfGlj+C4A9+lp/H7gH2H2q7Tip7oOAu4En99v0PX0tUz2vrwL+GXh4X8OvAjv2w75M1wU+OO8CLgB2BR420PbY/v5pwF0Dy37/RJ088HrbZmB+9y9jqnXq5/e2/v5T+u333/p5fwC4eFJt/wLsDOxD97o4YpptdP98B57rFf39A4GbgD0G6t5/im03sT4f7Z/fX6HrSfvlfvjbgYuAXYC9gCsmljFNTQV8qd+2+9C9nia2zYeAdwyM+xrgnzfnPQH8JfDtfj3T1/0L4/7f1+rNPfQFrKrupDvOPPGPYnW/J/Kofvj1VXVBVd1bVavp/ln/5qTZfKCqVlXVzcB/AF+rqm9W1b3AuXThPuiEqrqnqr4NfBw4aorSXgycV1XnVdV9VXUBsIwu5ADuAw5O8rCqWllVV23Cat9H9yHj3qr6CfAi4K1VdWu/jicALxkY/2fASVX1M+DTwG7A+6vqrn65V9GF9c+pqi9X1bf7dbiC7sPM5O036Gd9LT+rqvPogvLAJFvThewb++UuB949qc5bquoDVbWuXy+Ac6vq0qr6Kd1z8dOqOqOq1gOfYeC5qap/rKpb+lo/Q7dHfNhGtiXA7wH/UlUX98/5X9Ft4+nW7xfoAnl9X9udG5n/31bV2oF1muwLA8t+M91e92zsAb4IOLWqLuvn/cZ+3osHxnl7Vd1eVT+gC8dDNmM56+k+MByUZNuqWl5VN2xg/BOq60n7FvAtuoAEeCHwN1X1o6paAfzvGSz7Hf22/QHwPh54L54O/MFEDxDd6+wTA9Ntynvij4C3VNV11flWVa2ZQW3aDAb6AldV11TVS6tqL+Bgur2098H9XW2f7rsD7wQ+SffmHbRq4P5Ppng8+SSmmwbuf79f3mSPAV7Qd8ve3neJPpFuj/EeunA7FliZ5AtJfmnma8zqPuAm7NHXMV1Na/oAnFgf2Pg6ApDk15J8KV03+R19zZO336A19eDjxD/u570b3fkNk+vcc+Dx4HadMOPnJskfJrl8YHsfvJFaJ+wxuOz++ZnuH/YngH8DPp3u8MbfJdl2I/Ofar2mHF5VdwNrmfo1take9Lro572GB2/zHw7cn3iuNklVXQ+8lm5v/Nb+/bah+qdb5oOeBza+3SaPc//rvqq+RtdD85v9e+uxwOcHxt2U98TewIY+oGgWGei6X1VdS9c9eHDf9Ld0e+9PqKod6facpzxOuAkG9572AW6ZYpybgE9U1c4Dt+2r6u19nf9WVU+n626/lq53YaYm/7zgLXQfIDZW0+Y4k+4f4d5VtRPdcdvN2X630e0VTa7z5oHHm/2zif1x548Cx9F1h+4MXMnMal3JwHOa5OF0e+E/p+95OKGqDgJ+HXgW3SGADdW/sfUaXPYj6LqQb6ELJOi69yc8ehPm+6DXRZLt6dbr5mmnmN49G6iDqjqzurPtH9PX9Y7NWMZKuq72CTPppdjQe/F0uvf7S4DPTvoQvCluAvbfzGm1iQz0BSzJL6U7CWuv/vHedN1ul/Sj7EDX7Xt7kj3pjodtqb9K8vAkjwNeRtf1O9kngd9N8oz+pJqH9ifj7JXuJKpn9/9g7+3rWz/FPGbqLOAt/Qk/u9EdG5yt78PvAKytqp8mOQz4g82ZSb83dDZwUpId+gD+81msc3u6IFkNkO5kvIM3OMUDPgs8K8kTkzyE7hyAKf+vJPmtJI/vDyHcSfchZeK5WwXstxm1P3Ng2SfSHfK5qT98cjPw4v419HIeHCyrgL366aZyJvCyJIekO0nyb/p5L9+MGi/v69w1yaPp9sgBSHJgkqf0y/gp3d7t5ryezwbemGSX/r163Aym+ct+/L3pjpMPvhc/ATyPLtTP2Ix6JnwMODHJAek8IcmUH/i05Qz0he0u4NeAryW5hy7IrwQmzj4/ge6koDuALwD/NAvLvAi4HrgQeFdVnT95hKq6CXgO8Ca6kLmJ7sPEVv3tdXR7E2vpjkn/yRbU8za64/NX0J28c1nfNhv+BHhrkrvoPiicvQXz+lO6Pb3vAV+hC5xTt7hCoKqupjsm/190Qfd44D9nOO1VwKv7elYCPwKmO7v60XQfAO6ku97BRTzwoeT9wO/1Z2nP5PjvhDOB4+leC79Kd+x7wivpXjdr6E7W+urAsC/SHev9YZLbplivC+nOBzinX6/9gSM3oa5Bn6A73r0cOJ8HB+d2dCe03UbXnf5Iutf9pnor3Xa/Efh3uu28sa+ffg64lO4DxxeAUyYG9MfhL6P7oPcfm1HPhPfQve7Pp3veT6E7qU9DkKrN7qmTZqw/mehGYNvy+8TSUCX5Y+DIqpryJMwkBRzQH8Ofbh6n0p1s+ZYhlalZ5oVlJGmeS7I73SGL/wIOoOvF+uAWzG8x8Hx+/lsqmsPscpek+e8hwD/QHUb7Il13+oc2Z0ZJTqQ79PbOqrpx1irU0NnlLklSA9xDlySpAQa6JEkNmNcnxe222261ePHicZchSdLIXHrppbdV1aLJ7fM60BcvXsyyZcvGXYYkSSOT5PtTtdvlLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDVgXv/amtSspTuNcFl3jG5ZkobGPXRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSA4YW6ElOTXJrkisH2j6T5PL+tjzJ5X374iQ/GRj2kWHVJUlSi4Z5pbjTgA8CZ0w0VNXvT9xP8m5g8BJVN1TVIUOsR5KkZg0t0Kvq4iSLpxqWJMALgacMa/mSJC0k4zqG/iRgVVV9d6Bt3yTfTHJRkidNN2GSY5IsS7Js9erVw69UkqR5YFyBfhRw1sDjlcA+VXUo8OfAmUl2nGrCqjq5qpZU1ZJFixaNoFRJkua+kQd6km2A5wOfmWirqnurak1//1LgBuAXR12bJEnz1Tj20J8GXFtVKyYakixKsnV/fz/gAOB7Y6hNkqR5aZhfWzsL+C/gwCQrkryiH3QkD+5uB3gycEWSbwGfBY6tqrXDqk2SpNYM8yz3o6Zpf+kUbecA5wyrFkmSWueV4iRJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIasM24C5A0Zkt3GuGy7hjdsqQFxj10SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAUML9CSnJrk1yZUDbUuT3Jzk8v72zIFhb0xyfZLrkjxjWHVJktSiYe6hnwYcMUX7e6vqkP52HkCSg4Ajgcf103woydZDrE2SpKYMLdCr6mJg7QxHfw7w6aq6t6puBK4HDhtWbZIktWYcx9CPS3JF3yW/S9+2J3DTwDgr+jZJkjQDo75S3IeBE4Hq/74beDmQKcatqWaQ5BjgGIB99tlnOFVKGg6vSicNzUj30KtqVVWtr6r7gI/yQLf6CmDvgVH3Am6ZZh4nV9WSqlqyaNGi4RYsSdI8MdJAT7L7wMPnARNnwH8eODLJdkn2BQ4Avj7K2iRJms+G1uWe5CzgcGC3JCuA44HDkxxC152+HHgVQFVdleRs4GpgHfDqqlo/rNokSWrN0AK9qo6aovmUDYx/EnDSsOqRJKllXilOkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDVgm3EXIM0bS3cadwWSNC330CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAUML9CSnJrk1yZUDbe9Mcm2SK5Kcm2Tnvn1xkp8kuby/fWRYdUmS1KJh7qGfBhwxqe0C4OCqegLwHeCNA8NuqKpD+tuxQ6xLkqTmDC3Qq+piYO2ktvOral3/8BJgr2EtX5KkhWScx9BfDvzrwON9k3wzyUVJnjTdREmOSbIsybLVq1cPv0pJkuaBsQR6kjcD64BP9U0rgX2q6lDgz4Ezk+w41bRVdXJVLamqJYsWLRpNwZIkzXEjD/QkRwPPAl5UVQVQVfdW1Zr+/qXADcAvjro2SZLmq5EGepIjgNcDz66qHw+0L0qydX9/P+AA4HujrE2SpPlsaD+fmuQs4HBgtyQrgOPpzmrfDrggCcAl/RntTwbemmQdsB44tqrWTjljSZL0c4YW6FV11BTNp0wz7jnAOcOqRZKk1nmlOEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDdhm3AVI0lAs3WmEy7pjdMuSpuEeuiRJDTDQJUlqwIwCPcmFM2mTJEnjscFAT/LQJLsCuyXZJcmu/W0xsMdGpj01ya1Jrhxo2zXJBUm+2//dZWDYG5Ncn+S6JM/YwvWSJGlB2dge+quAS4Ff6v9O3D4H/P1Gpj0NOGJS2xuAC6vqAODC/jFJDgKOBB7XT/OhJFvPeC0kSVrgNhjoVfX+qtoX+Iuq2q+q9u1vv1JVH9zItBcDayc1Pwc4vb9/OvDcgfZPV9W9VXUjcD1w2CauiyRJC9aMvrZWVR9I8uvA4sFpquqMTVzeo6pqZT/tyiSP7Nv3BC4ZGG9F3/ZzkhwDHAOwzz77bOLiJUlq04wCPckngP2By4H1fXMBmxro0y5iiraaasSqOhk4GWDJkiVTjiNJ0kIz0wvLLAEOqqotDdBVSXbv9853B27t21cAew+MtxdwyxYuS5KkBWOm30O/Enj0LCzv88DR/f2j6U6um2g/Msl2SfYFDgC+PgvLkyRpQZjpHvpuwNVJvg7cO9FYVc+eboIkZwGH033lbQVwPPB24OwkrwB+ALygn89VSc4GrgbWAa+uqvVTzliSJP2cmQb60k2dcVUdNc2gp04z/knASZu6HEmSNPOz3C8adiGSJGnzzfQs97t44KzzhwDbAvdU1Y7DKkySJM3cTPfQdxh8nOS5eOEXSZLmjM36tbWq+j/AU2a3FEmStLlm2uX+/IGHW9F9L92LukiSNEfM9Cz33x24vw5YTnf9dUmSNAfM9Bj6y4ZdiCRJ2nwzOoaeZK8k5/a/b74qyTlJ9hp2cZIkaWZmelLcx+kuz7oH3a+g/XPfJkmS5oCZBvqiqvp4Va3rb6cBi4ZYlyRJ2gQzDfTbkrw4ydb97cXAmmEWJkmSZm6mgf5y4IXAD4GVwO8BnignSdIcMdOvrZ0IHF1VPwJIsivwLrqglyRJYzbTPfQnTIQ5QFWtBQ4dTkmSJGlTzTTQt0qyy8SDfg99pnv3kiRpyGYayu8Gvprks3SXfH0h/na5JElzxkyvFHdGkmV0P8gS4PlVdfVQK5MkSTM2427zPsANcUmS5qDN+vlUSZI0txjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGjDj30OfLUkOBD4z0LQf8NfAzsArgdV9+5uq6rzRVidJ0vw08kCvquuAQwCSbA3cDJwLvAx4b1W9a9Q1SZI03427y/2pwA1V9f0x1yFJ0rw27kA/Ejhr4PFxSa5IcmqSXcZVlCRJ883YAj3JQ4BnA//YN30Y2J+uO34l8O5ppjsmybIky1avXj3VKJIkLTjj3EP/HeCyqloFUFWrqmp9Vd0HfBQ4bKqJqurkqlpSVUsWLVo0wnIlSZq7xhnoRzHQ3Z5k94FhzwOuHHlFkiTNUyM/yx0gycOBpwOvGmj+uySHAAUsnzRMkiRtwFgCvap+DPzCpLaXjKMWSZJaMO6z3CVJ0iww0CVJaoCBLklSA8ZyDF2aNUt3GncF0mhfh0vvGN2yNK+4hy5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQHbjGOhSZYDdwHrgXVVtSTJrsBngMXAcuCFVfWjcdQnSdJ8M8499N+qqkOqakn/+A3AhVV1AHBh/1iSJM3AXOpyfw5wen//dOC54ytFkqT5ZVyBXsD5SS5Nckzf9qiqWgnQ/33kVBMmOSbJsiTLVq9ePaJyJUma28ZyDB34jaq6JckjgQuSXDvTCavqZOBkgCVLltSwCpQkaT4Zyx56Vd3S/70VOBc4DFiVZHeA/u+t46hNkqT5aOSBnmT7JDtM3Ad+G7gS+DxwdD/a0cDnRl2bJEnz1Ti63B8FnJtkYvlnVtX/TfIN4OwkrwB+ALxgDLVJkjQvjTzQq+p7wK9M0b4GeOqo65EkqQVz6WtrkiRpMxnokiQ1wECXJKkBBrokSQ0w0CVJasC4rhQnSdocS3ca4bLuGN2ytMXcQ5ckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAfz5Vs2+UP+8oSQLcQ5ckqQkGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaMPJAT7J3ki8luSbJVUle07cvTXJzksv72zNHXZskSfPVOH4+dR3wuqq6LMkOwKVJLuiHvbeq3jWGmiRJmtdGHuhVtRJY2d+/K8k1wJ6jrkOSpJaM9Rh6ksXAocDX+qbjklyR5NQku0wzzTFJliVZtnr16lGVKknSnDa2QE/yCOAc4LVVdSfwYWB/4BC6Pfh3TzVdVZ1cVUuqasmiRYtGVa4kSXPaWAI9ybZ0Yf6pqvongKpaVVXrq+o+4KPAYeOoTZKk+WgcZ7kHOAW4pqreM9C++8BozwOuHHVtkiTNV+M4y/03gJcA305yed/2JuCoJIcABSwHXjWG2iRJmpfGcZb7V4BMMei8UdeyoCzdadwVSJKGyCvFSZLUAANdkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqwDi+hy5Jmg9G+XXXpXeMblmNcg9dkqQGGOiSJDXAQJckqQEGuiRJDTDQJUlqgIEuSVIDDHRJkhpgoEuS1AADXZKkBhjokiQ1wECXJKkBBrokSQ0w0CVJaoCBLklSAwx0SZIaYKBLktQAA12SpAYY6JIkNcBAlySpAQa6JEkN2GbcBSxoS3cadwWStPCM8n/v0jtGtij30CVJaoCBLklSA+xylySNn4cgt5h76JIkNcBAlySpAXMu0JMckeS6JNcnecO465EkaT6YU4GeZGvg74HfAQ4Cjkpy0HirkiRp7ptrJ8UdBlxfVd8DSPJp4DnA1SOrwBMzJEnz0JzaQwf2BG4aeLyib5MkSRsw1/bQM0VbPWiE5BjgmP7h3UmuG3pVm2434LZxFzEGrvfC4novLK735jhhqljbYo+ZqnGuBfoKYO+Bx3sBtwyOUFUnAyePsqhNlWRZVS0Zdx2j5novLK73wuJ6z31zrcv9G8ABSfZN8hDgSODzY65JkqQ5b07toVfVuiTHAf8GbA2cWlVXjbksSZLmvDkV6ABVdR5w3rjr2EJz+pDAELneC4vrvbC43nNcqmrjY0mSpDltrh1DlyRJm8FAH4IkJya5IsnlSc5Psse4axqFJO9Mcm2/7ucm2XncNY1CkhckuSrJfUnmxdmwW2KhXp45yalJbk1y5bhrGaUkeyf5UpJr+tf5a8Zd0ygkeWiSryf5Vr/eJ4y7po2xy30IkuxYVXf29/8MOKiqjh1zWUOX5LeBL/YnN74DoKpeP+ayhi7JLwP3Af8A/EVVLRtzSUPTX575O8DT6b5m+g3gqKoa3dUcxyTJk4G7gTOq6uBx1zMqSXYHdq+qy5LsAFwKPLf15zxJgO2r6u4k2wJfAV5TVZeMubRpuYc+BBNh3tueSRfHaVVVnV9V6/qHl9BdR6B5VXVNVc3FCxwNw/2XZ66q/wdMXJ65eVV1MbB23HWMWlWtrKrL+vt3AdewAK7gWZ27+4fb9rc5/b/cQB+SJCcluQl4EfDX465nDF4O/Ou4i9Cs8/LMC1iSxcChwNfGXMpIJNk6yeXArcAFVTWn19tA30xJ/j3JlVPcngNQVW+uqr2BTwHHjbfa2bOx9e7HeTOwjm7dmzCT9V4gNnp5ZrUpySOAc4DXTuqFbFZVra+qQ+h6Gw9LMqcPtcy576HPF1X1tBmOeibwBeD4IZYzMhtb7yRHA88CnloNnaCxCc936zZ6eWa1pz+GfA7wqar6p3HXM2pVdXuSLwNHAHP2pEj30IcgyQEDD58NXDuuWkYpyRHA64FnV9WPx12PhsLLMy8w/clhpwDXVNV7xl3PqCRZNPFNnSQPA57GHP9f7lnuQ5DkHOBAujOfvw8cW1U3j7eq4UtyPbAdsKZvumSBnN3/POADwCLgduDyqnrGWIsaoiTPBN7HA5dnPmm8FY1GkrOAw+l+fWsVcHxVnTLWokYgyROB/wC+Tfc/DeBN/VU9m5XkCcDpdK/zrYCzq+qt461qwwx0SZIaYJe7JEkNMNAlSWqAgS5JUgMMdEmSGmCgS5LUAANdkqQGGOiSJDXAQJckqQH/H+Rs+go/ymb0AAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" }, "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\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yay! We learned how to sample from a `pymc` distribution!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What is going on behind the scenes?" ] }, { "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", "aesara.dprint(z)" ] }, { "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": [ "aesara.dprint(model.basic_RVs[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can try to sample via {meth}`~aesara.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.14882183 4.14157407]\n", "Sample 1: [-0.14882183 4.14157407]\n", "Sample 2: [-0.14882183 4.14157407]\n", "Sample 3: [-0.14882183 4.14157407]\n", "Sample 4: [-0.14882183 4.14157407]\n", "Sample 5: [-0.14882183 4.14157407]\n", "Sample 6: [-0.14882183 4.14157407]\n", "Sample 7: [-0.14882183 4.14157407]\n", "Sample 8: [-0.14882183 4.14157407]\n", "Sample 9: [-0.14882183 4.14157407]\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {z.eval()}\")" ] }, { "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: [0.43562858 0.2937141 ]\n", "Sample 1: [-0.82913473 -0.45162817]\n", "Sample 2: [-0.39186432 0.84863972]\n", "Sample 3: [-0.9840986 -0.29480585]\n", "Sample 4: [-0.13096991 -0.10632378]\n", "Sample 5: [-0.94007078 3.080421 ]\n", "Sample 6: [-0.32518418 1.52276117]\n", "Sample 7: [-0.21348089 3.09445457]\n", "Sample 8: [-0.62399986 -2.13572382]\n", "Sample 9: [-0.21707549 1.44738996]\n" ] } ], "source": [ "for i in range(10):\n", " print(f\"Sample {i}: {pm.draw(z)}\")" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAHiCAYAAAA06c+jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeAklEQVR4nO3de5SkdX3n8c+nq3u65woMwwjMjAwqXpGoZ8Br1ohgUAme3dxM1GiSPXPWLFlzDkoUdhM8m2Q5cY9RV7M5E0U3QiSuQsLhEi7rbbOrCCJXB1zWADMMl5kB5j59q+/+UYXbNj0zDM+359tT9X6dM4fpqeJT3+ep7vrU81RX/RwRAgAANQaqBwAAoJ9RxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYmENsX2T70oLbvcD25w/17QKgiAFJku032f7ftrfZfsL2/7J9avVcGfZV7rbD9oskKSL+LCL+9bPI+pbtA14PwLM3WD0AUM32EklXS/qgpK9Kmifp5yWNVs7Vb2wPRsRE9RzAocYRMSC9WJIi4isRMRkReyLihoi4U5Jsv9D2N2xvtb3F9mW2j3z6f7b9gO2P2L7T9i7bX7D9PNvX2d5h+ybbR3Wvu7p7JLrW9ibbj9g+b1+D2X5d90j9Kdt32P6FKZd9wPZPurfxz7bf81x3wNSjZtsjti/tbu9Ttm/pbs+fqvME5bO2d9r+bPf6b+heZ1v3v2+Yknui7e9M2Q+fm3I7T++L37X9kKRvdP/9v9t+tJv3HduvmJL3Jdt/2d23O7tnLo61/SnbT9q+1/arn+t+ACpQxID0Y0mTtv+b7bc/XZpTWNJ/knS8pJdJWiXpomnX+WVJZ6pT6r8k6TpJF0haps7P2b+bdv23SDpJ0tskfdT2GdOHsr1C0jWS/kTSUkkflvR128fYXijpM5LeHhGLJb1B0u0HveUze7+kI9TZzqMl/RtJeyLiQkn/U9K5EbEoIs61vbQ742e61/2kpGtsH93N+ltJ3+9edpGk981we29WZ7/+Yvfr69TZN8sl3SbpsmnX/zVJ/16dfTsq6bvd6y2T9LXuDMBhgyJG34uI7ZLeJCkk/bWkzbavsv287uX3R8SNETEaEZvVeaB/87SY/xIRj0XEw+qU1c0R8cOIGJV0paTpR2kfj4hdEXGXpC9K+o0ZRnuvpGsj4tqIaEfEjZJulfSO7uVtSSfbnh8Rj0TEPfvZzF/rHt3+9M9+rjuuTnG+qHuG4AfdfTSTd0r6PxHx5YiYiIivSLpX0i/Zfr6kUyX9UUSMRcQ/SbpqhoyLuvtijyRFxCURsaO77y6S9HO2j5hy/Su7M+1VZ9/ujYi/iYhJSX+nZ+5rYE6jiAFJEbE+Ij4QESslnazO0e+nJMn2ctuX237Y9nZJl6pz9DXVY1P+vmeGrxdNu/6GKX9/sHt7050g6VenleebJB0XEbsk/bo6R6uP2L7G9kv3s4lfjYgjp/7Zz3W/LOl6SZd3T5//ue2hfVz3+O78Uz0oaUX3siciYveUyzbomX76b7Zbti+2/X+7+/qB7kVT9/fB7mtgTqOIgWki4l5JX1KnkKXOaemQdEpELFHnSNUNb2bVlL8/X9KmGa6zQdKXpxXowoi4uDvn9RFxpqTj1DkK/euGM6mbOx4RH4+Il6tzyvtsSb/19MXTrr5JnScMUz1f0sOSHpG01PaCKZet0jNNzfxNSe+SdIY6p8dXd/+96f4G5iyKGH3P9kttn2d7ZffrVeqcKv5e9yqLJe2U9FT3dduPJNzsf7C9oPuLSL+tzinV6S5V5xTvL3aPFEds/4Ltld1fnjqn+1rxaHe+yYS5ZPsttl9puyVpuzqnqp/OfkzSC6Zc/VpJL7b9m7YHbf+6pJdLujoiHlTnVPpFtufZfr06r5/vz+Lu9myVtEDSn2VsEzCXUcSAtEPSayXdbHuXOgV8t6Snf5v545JeI2mbOr+YdEXCbX5b0v2S/oek/xwRN0y/QkRsUOfo8AJJm9U5Qv6IOj+3A935Nkl6Qp3XrH8vYS5JOladX3raLml9d9an34f8aUm/0v0N5c9ExFZ1jpjPU6c8z5d0dkRs6V7/PZJe373sT9R5wrG/t4X9jTqnth+W9CP9/ydDQM9yxPQzTQBmi+3Vkv5Z0lA/vmfW9t9Jujci/rh6FmCu4IgYwKyxfao778MesH2WOkf4f188FjCn8MlaAGbTseqcyj9a0kZJH4yIH9aOBMwtnJoGAKAQp6YBAChEEQMAUKjkNeJ5Ho4RLay4aQAADrkdenJLRBwz02UlRTyihXqt31px0wAOFSedcIt2Tg5Q6Kb42vSPgv0pTk0DAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABQqWY8YQB/IWke4F9c17sVtwnPGETEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFBqsHAJDESc+ro52TM8d4cCglJybGE0KS9jH3eU/giBgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAocHqAYB+l7VgvQeckpMl2pGSk7VdnjcvJae9t3lGTE42D0HPSDkitn2k7a/Zvtf2etuvz8gFAKDXZR0Rf1rSP0bEr9ieJ2lBUi4AAD2tcRHbXiLpX0j6gCRFxJiksaa5AAD0g4xT0y+QtFnSF23/0PbnbS9MyAUAoOdlFPGgpNdI+q8R8WpJuyR9dPqVbK+1favtW8c1mnCzAAAc/jKKeKOkjRFxc/frr6lTzD8jItZFxJqIWDOk4YSbBQDg8Ne4iCPiUUkbbL+k+09vlfSjprkAAPSDrN+a/n1Jl3V/Y/onkn47KRcAgJ6WUsQRcbukNRlZAAD0Ez7iEgCAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAoaxP1kKvctJztWjn5CTN41YrJSdDTE6m5Hgo5zPcPZj0sDAxkZOTtH9iLGd1Vg813z8xMZ4wifJ+HgaHUnKyvpfTZD3uzDKOiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUSloBHD0ra2HtpAXMBxIWZZekaEdKTsYi8QOLFiZMIilrUfZWKyXGSTlzTcZ9nrVvIuk+j/GJlJy0x4s+wxExAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEI5q6yjZ3lwqHqEn5W0oPrAcM63vhfMT8lJMZG0SPxxy1JyPJq02Hwr53gh5iU93G3d1jjC80cSBpEmN2xKyWkdsSQlZ3Lb9pQcRTsn5zDBETEAAIUoYgAAClHEAAAUoogBAChEEQMAUCitiG23bP/Q9tVZmQAA9LrMI+IPSVqfmAcAQM9LKWLbKyW9U9LnM/IAAOgXWUfEn5J0vqT+ehc2AAANNS5i22dLejwifnCA6621favtW8c12vRmAQDoCRlHxG+UdI7tByRdLul025dOv1JErIuINRGxZkjDCTcLAMDhr3ERR8THImJlRKyW9G5J34iI9zaeDACAPsD7iAEAKJS6+lJEfEvStzIzAQDoZRwRAwBQiCIGAKAQRQwAQKHU14jRe2JiPCWntWhRSs7krt0pOYOrjk7JicULmofYzTMkeXQiJWd82cKUnJ0rct6mOLx9MiVnYCxScoYHmx+/DDz2ZMIk0uDzV6TktB/fkpIzMJJzn7f35nzWhFutxhlZj4H7wxExAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEKD1QMAB2Nw2dKcoNGxlJhYdkTzjKHmi5dL0s6XHpWSM3pkzvPzvUl31e7lOQ9TI09GSs7YksWNMxYl3ecDu0ZzcrQsJSe2PpmS41bO/onJyZSc2cYRMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCOStuI4/n1nOjgZHhnCA7JydJ+/ichdBHj1nQOCOS7vKtJ+cETc6PlJz2cFJOzhrxGtuas39Gtjb/Xm4PLUyYRBrZMi8n59GcfePxiZycycmUnJgYT8mZbXPrUR8AgD5DEQMAUIgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQiCIGAKBQ4yK2vcr2N22vt32P7Q9lDAYAQD/IWPRhQtJ5EXGb7cWSfmD7xoj4UUI2AAA9rfERcUQ8EhG3df++Q9J6SSua5gIA0A9SXyO2vVrSqyXdnJkLAECvSluP2PYiSV+X9AcRsX2Gy9dKWitJI2q+hisAAL0g5YjY9pA6JXxZRFwx03UiYl1ErImINUNKWmweAIDDXOMjYtuW9AVJ6yPik81HOox57rwbbGAk58mOB3NOmnjJ4pQczR9JiRk7OueszN6jm++fx0+LhEkkL9+TknPmSetTcu7btjwl54FHlqXk6Imc7509z2ueMTDu5iGSBve0UnLSjI6lxHg46fFrfCIlJ8X4vi/KaI43SnqfpNNt3979846EXAAAel7jp/MR8U+Scp7eAQDQZ+bOuVQAAPoQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUyln1vYqTnkdEe27lJHArZ8HwGB1NybEWp+S0l8xPyfFkpOTsPK75CqAxkDPLymOeTMkZcs738cUvvCIl54tHvCkl577jlqfkPPT40Y0zRnfnfB8veiQlRrtPyPn5XPjUjpQcTyY9lu7clZMzyzgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoNBg9QCNRLt6gp/l5s9rBkaGEwaR2ntHU3JaRx2RkqOxsZyciJSYbS/M2c+ebJ7RWr6neYikNy+/PyVn5bwnUnJOmZcSoyHn/Jy/87i7U3L+8tE3N87YfVzONi27K+fnYeH6LSk52rM3JSbGJ1JyPNS84rJm2R+OiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAACqUUse2zbN9n+37bH83IBACgHzReNdl2S9LnJJ0paaOkW2xfFRE/2v//mPAcIHIW187iVqtxRtqC2ANOyZl8cltKTuv456XkeHwyJWdod86C6ttf0Hw/Dw/n3OcnDm9Oyfn2Uy9JyfmtJQ+l5Hz6+FtScn7v4dem5Bx91M7GGdvvn58wibTz2JyTmkPbjkjJmTc+npKjbTtSYjIeT2My5zFnfzLuxdMk3R8RP4mIMUmXS3pXQi4AAD0vo4hXSNow5euN3X8DAAAH0PjUtKSZzs0947yf7bWS1krSiBYk3CwAAIe/jCPijZJWTfl6paRN068UEesiYk1ErBnScMLNAgBw+Mso4lsknWT7RNvzJL1b0lUJuQAA9LzGp6YjYsL2uZKul9SSdElE3NN4MgAA+kDGa8SKiGslXZuRBQBAP+GTtQAAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAolPLJWs9JtMtuerbERMKi2M56btRKSfFQ0rfI8Nxa6GNwT873X2tP8/2864mc1ci+v+PElJz5rbGUnJv2LE7J2dGen5Iz2s75Xt68eUnjjIEjc77/lnw3Z9H6wR2jKTlZYvfu6hF+yq2cx1Lt5y7niBgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCg9UD9BTzvGafRkdTYjw2lpIzMH5ESs687c0zYuNQ8xBJ/zj4ipScV5ywKSXn/Id+OSXnRcu2pOTsHJ+XkjP84HDzjCcSBpHUHmqn5AzsGU/J0fhESky0IydnImm7ZhnNAQBAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQiCIGAKDQYPUAPSUSFun23Hpu5OHmi6BLUmzfmZLjZUtTcgbGcxYe92TzjKxF4icW5Cx8f/foqpScI4/dnpJz190npOTM29rKyUnYrEWbEh4rJC28/8mUHO3Zm5OTxEM51TQw0vzxa3JnzmPX/sytR30AAPoMRQwAQCGKGACAQhQxAACFGhWx7U/Yvtf2nbavtH1k0lwAAPSFpkfEN0o6OSJOkfRjSR9rPhIAAP2jURFHxA0RMdH98nuSVjYfCQCA/pH5GvHvSLouMQ8AgJ53wHdN275J0rEzXHRhRPxD9zoXSpqQdNl+ctZKWitJI1rwnIYFAKDXHLCII+KM/V1u+/2Szpb01ojY58cVRcQ6SeskaYmX5nysEQAAh7lGnyNm+yxJfyjpzRGxO2ckAAD6R9PXiD8rabGkG23fbvuvEmYCAKBvNDoijogXZQ0CAEA/4pO1AAAoRBEDAFCIIgYAoBBFDABAoUa/rIVZEO2UGA/k3LWxe269K83bcrZrcNeSlJzlt403zth+4vyESaT5W1JitOVVrZSc0Q1LU3KOfiQlRuMLc3IG9zbP8OTc+iiF2LajeoSfEeMTB77Ss9DesyclZ7ZxRAwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgUM4q65hz2mNjKTkeHMrJGXBKTtaC4YO335+So9UrGkcsvXl7wiDSrpctS8lZ8e2UGO1Z1krJGdoVKTnD21JitPDhvY0zBh99qvkgkjQwt46lYvfulJy0x4uMECft4/0MM7fuRQAA+gxFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQaLB6AMySpMWsY2I8JyclRWoND6fkeP5ISo4e29o4IpYvTRhEWnjHppScyeNy5lnwwGhKTnthzn3e2rozJSdGhpqHTEw0z5CkicmUGLdaKTlZP+ft8aT9kyHas34THBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUCiliG1/2HbYXpaRBwBAv2hcxLZXSTpT0kPNxwEAoL9kHBH/haTzlfcxowAA9I1GRWz7HEkPR8QdSfMAANBXDrj6ku2bJB07w0UXSrpA0tuezQ3ZXitprSSNaMFBjAgAQO86YBFHxBkz/bvtV0o6UdIdtiVppaTbbJ8WEY/OkLNO0jpJWuKlnMYGAEAN1iOOiLskLX/6a9sPSFoTEVsS5gIAoC/wPmIAAAo95yPi6SJidVYWEkQ7J8c5z9XcaqXkxOhoSk57z96UnIH5I40zvOWp5oNIil27U3JanZeamtubc1+1NqfESBOTKTEpe2cgZx+3d+7KyUm6r7J+zvsNR8QAABSiiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAACg1WD4A5Lto5MRM5ORocSolx0sLsMdp8QfWYmEiYRPJg0o9zRE7M7j05OVn7J2nR+vbe5vd51vdfe2wsJSdL2s95FiccayY9Bu4PR8QAABSiiAEAKEQRAwBQiCIGAKAQRQwAQCGKGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUGqweADgYMTGek5OSIskJz2XHJ5pnSBoYmkzJiUcfT8lpJ21Xa9HClJzJnbtScgaGmj9stsfGEiZRzvefJEU7JydLr27XPnBEDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABSiiAEAKNS4iG3/vu37bN9j+88zhgIAoF80WljT9lskvUvSKRExant5zlgAAPSHpitcf1DSxRExKkkRkbOiOHpP1kLfWbIWDM/ISdo37fGJlBy3Wik5WSZ37soJSrrPU/Zzny18f9B6dbv2oel3w4sl/bztm21/2/apGUMBANAvDnhEbPsmScfOcNGF3f//KEmvk3SqpK/afkFExAw5ayWtlaQRLWgyMwAAPeOARRwRZ+zrMtsflHRFt3i/b7staZmkzTPkrJO0TpKWeOkzihoAgH7U9NT030s6XZJsv1jSPElbGmYCANA3mv6y1iWSLrF9t6QxSe+f6bQ0AACYWaMijogxSe9NmgUAgL4zx95TAgBAf6GIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBACjU9LOmgWenVxf6zljgPWvfJC02H5OTKTlz7j5P2j9zbrtw2OOIGACAQhQxAACFKGIAAApRxAAAFKKIAQAoRBEDAFCIIgYAoBBFDABAIYoYAIBCFDEAAIUoYgAAClHEAAAUoogBAChEEQMAUIgiBgCgEEUMAEAhihgAgEKOiEN/o/ZmSQ8e4ptdJmnLIb7NCmxnb+mH7eyHbZTYzl5zsNt5QkQcM9MFJUVcwfatEbGmeo7Zxnb2ln7Yzn7YRont7DWZ28mpaQAAClHEAAAU6qciXlc9wCHCdvaWftjOfthGie3sNWnb2TevEQMAMBf10xExAABzTl8Wse0P2w7by6pnmQ22/6PtO23fbvsG28dXzzQbbH/C9r3dbb3S9pHVM2Wz/au277Hdtt1zv4lq+yzb99m+3/ZHq+eZDbYvsf247burZ5kttlfZ/qbt9d3v1w9VzzQbbI/Y/r7tO7rb+fGM3L4rYturJJ0p6aHqWWbRJyLilIh4laSrJf1R8Tyz5UZJJ0fEKZJ+LOljxfPMhrsl/StJ36keJJvtlqTPSXq7pJdL+g3bL6+dalZ8SdJZ1UPMsglJ50XEyyS9TtK/7dH7clTS6RHxc5JeJeks269rGtp3RSzpLySdL6lnXxyPiO1TvlyoHt3WiLghIia6X35P0srKeWZDRKyPiPuq55glp0m6PyJ+EhFjki6X9K7imdJFxHckPVE9x2yKiEci4rbu33dIWi9pRe1U+aJjZ/fLoe6fxo+vfVXEts+R9HBE3FE9y2yz/ae2N0h6j3r3iHiq35F0XfUQOCgrJG2Y8vVG9eCDd7+xvVrSqyXdXDzKrLDdsn27pMcl3RgRjbdzsPFUc4ztmyQdO8NFF0q6QNLbDu1Es2N/2xkR/xARF0q60PbHJJ0r6Y8P6YBJDrSd3etcqM6pscsO5WxZns029ijP8G89efamX9heJOnrkv5g2pm5nhERk5Je1f2dlCttnxwRjV7/77kijogzZvp326+UdKKkO2xLndOYt9k+LSIePYQjptjXds7gbyVdo8O0iA+0nbbfL+lsSW+Nw/S9eAdxX/aajZJWTfl6paRNRbOgIdtD6pTwZRFxRfU8sy0inrL9LXVe/29UxH1zajoi7oqI5RGxOiJWq/Mg8JrDsYQPxPZJU748R9K9VbPMJttnSfpDSedExO7qeXDQbpF0ku0Tbc+T9G5JVxXPhOfAnaObL0haHxGfrJ5nttg+5ul3Z9ieL+kMJTy+9k0R95mLbd9t+051TsX35FsJJH1W0mJJN3bfqvVX1QNls/0vbW+U9HpJ19i+vnqmLN1ftDtX0vXq/HLPVyPintqp8tn+iqTvSnqJ7Y22f7d6plnwRknvk3R692fxdtvvqB5qFhwn6Zvdx9Zb1HmN+OqmoXyyFgAAhTgiBgCgEEUMAEAhihgAgEIUMQAAhShiAAAKUcQAABSiiAEAKEQRAwBQ6P8Bkf7MQQFmhvQAAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" }, "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\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Enough with Random Variables, I want to see some (log)probabilities!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall we have defined the following model above:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": "z ~ N(, )", "text/latex": "$$\n \\begin{array}{rcl}\n \\text{z} &\\sim & \\operatorname{N}(\\text{},~\\text{})\n \\end{array}\n $$" }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model" ] }, { "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 = at.vector(name=\"z\")\n", "z_logp = pm.logp(rv=z, value=z_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`z_logp` contains the Aesara 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", " |Elemwise{gt,no_inplace} [id U]\n", " |TensorConstant{[1. 2.]} [id L]\n", " |InplaceDimShuffle{x} [id V]\n", " |TensorConstant{0.0} [id W]\n" ] }, { "data": { "text/plain": "" }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "aesara.dprint(z_logp)" ] }, { "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`." ] }, { "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}`~aesara.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]})" ] }, { "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]))" ] }, { "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": [ "Elemwise{mul,no_inplace} [id A]\n", " |Check{sigma > 0} [id B] 'z_logprob'\n", " | |Elemwise{sub,no_inplace} [id C]\n", " | | |Elemwise{sub,no_inplace} [id D]\n", " | | | |Elemwise{mul,no_inplace} [id E]\n", " | | | | |InplaceDimShuffle{x} [id F]\n", " | | | | | |TensorConstant{-0.5} [id G]\n", " | | | | |Elemwise{pow,no_inplace} [id H]\n", " | | | | |Elemwise{true_div,no_inplace} [id I]\n", " | | | | | |Elemwise{sub,no_inplace} [id J]\n", " | | | | | | |z [id K]\n", " | | | | | | |TensorConstant{(2,) of 0} [id L]\n", " | | | | | |TensorConstant{[1. 2.]} [id M]\n", " | | | | |InplaceDimShuffle{x} [id N]\n", " | | | | |TensorConstant{2} [id O]\n", " | | | |InplaceDimShuffle{x} [id P]\n", " | | | |Elemwise{log,no_inplace} [id Q]\n", " | | | |Elemwise{sqrt,no_inplace} [id R]\n", " | | | |TensorConstant{6.283185307179586} [id S]\n", " | | |Elemwise{log,no_inplace} [id T]\n", " | | |TensorConstant{[1. 2.]} [id M]\n", " | |All [id U]\n", " | |Elemwise{gt,no_inplace} [id V]\n", " | |TensorConstant{[1. 2.]} [id M]\n", " | |InplaceDimShuffle{x} [id W]\n", " | |TensorConstant{0.0} [id X]\n", " |InplaceDimShuffle{x} [id Y]\n", " |TensorConstant{1.0} [id Z]\n" ] }, { "data": { "text/plain": "" }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "aesara.dprint(model.logp(sum=False))" ] }, { "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 Aesara function of the model logp." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "logp_function = model.compile_logp(sum=False)" ] }, { "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)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What are value variables and why are they important?\n", "\n", "As he 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([ 0.14809949, -1.18242809, 1.34229652])" }, "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)" ] }, { "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)" ] }, { "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" ] }, { "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" ] }, { "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 aesara.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 = at.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})" ] }, { "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", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{Note}\n", "For `sigma_log_value` we add the $-10$ term for the `scipy` and `aesara` to match because of the jacobian.\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we already saw, we can also use the method {meth}`~pymc.Model.compile_logp` to obtain a compiled aesara 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})" ] }, { "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want to go deeper into the internals of `aesara` RandomVariables and `pymc` distributions please take a look into the [distribution developer guide](implementing-a-distribution)." ] } ], "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": { "name": "pymc", "language": "python", "display_name": "pymc" }, "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.13" }, "vscode": { "interpreter": { "hash": "867ba48c05011db76db56a12fb95ccd32f7ac276df8f4ae698e0d475911a6ba0" } } }, "nbformat": 4, "nbformat_minor": 1 }