Created
December 24, 2021 13:06
-
-
Save twiecki/0219d3fa059776b90af0f08b5452e346 to your computer and use it in GitHub Desktop.
PyMC3 vs PyMC 4.0.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "metadata": {}, | |
| "id": "3e991cb4", | |
| "cell_type": "markdown", | |
| "source": "# PyMC 4.0 Release Announcement\n\nWe, the PyMC core development team, are incredibly excited to announce the release of a major rewrite of PyMC (formerly PyMC3): `4.0.0`. There are many exciting new updates that we will talk about in this and upcoming blog posts.\n\n## It's now called PyMC instead of PyMC3\nFirst, the biggest news: **PyMC3 has been renamed to PyMC with this release. PyMC3 version 3.x will stay under the current name to not break production systems.** There are a few reasons for this. The main one is that PyMC3 4.0 is quite confusing.\n\n## What about PyMC4?\nAs described in our previous post [\"The Future of PyMC3, or: Theano is Dead, Long Live Theano\"](https://pymc-devs.medium.com/the-future-of-pymc3-or-theano-is-dead-long-live-theano-d8005f8a0e9b), we have experimented with (but given up on) a PyMC4 that was based on TensorFlow and TensorFlow Probability.\n\nOne of the major reasons for this was (i) that TF was really difficult to work with, partly because of its eager execution mode, but, more importantly, (ii) we realized that `Theano` already was in a great spot to take `PyMC` to the next level. \n\nSo to summarize, here are the names and what they refer to:\n* PyMC3 3.x: The version you all love and use, version numbers start with 3.x.\n* PyMC 4.0: The successor to PyMC3, (mostly) API compatible, our main focus going forward, and what this blog post is about. This is the only version that will matter going forward. Version numbers start with 4.x.\n* PyMC4: The discontinued experiment we did with new API on the TensorFlow backend.\n\n## Theano -> Theano-PyMC -> Aesara\n\nSo after we rediscovered the power of `Theano`, we went ahead and forked `Theano` to `Theano-PyMC` and made `PyMC3` rely on this fork instead.\n\nDuring this forking process we discovered that we could take this project much further still, so we undertook a massive rewrite of `Theano-PyMC` (this charge was led by Brandon Willard), removing swaths of old and obscure code, and restructuring the entire library to be more developer friendly.\n\nThis rewrite motivated a rename of the package to [`Aesara`](https://github.com/aesara-devs/aesara) (Theano's daughter), and is now a project largely separate from PyMC, although PyMC is the main customer.\n\nMost imporantly, however, we added two new computational backends: `JAX` and `numba`. The way this works is that `aesara` is best understood as a computational graph library that allows you to build a computational graph out of array-operations (additions, multiplications, dot-products, indexing, for-loops). From that graph, you can do graph optimizations (`log(exp(x)) -> x`), symbolic rewrites (`N(0, 1) + a` -> `N(a, 1)`), as well as compile that graph, and here we have a ton of flexibility now and can add new compilation backends easily. Currently there is a `C`, `CUDA`, `JAX` (with GPU support), and `numba` backend. Future blog posts will go into more detail as to how to use these backends and their different trade-offs.\n\n## New Features\n\nMost of the changes in `v4` happened under the hood and lead to lots of code simplifications and setting the stage for future innovations that were not possible before.\n\nHowever, there are quite a few good reasons to check out this new version already.\n\n### Faster predictive sampling\n\nIn `v3` we were not really using the full power of `Theano`. While calling e.g. `x = pm.Normal('x')` created a `theano.TensorVariable` which you could manipulate further, like inputting it into another random variable. However, we were not really building a full-fledged graph. This created all kinds of issues but a major one was a very complicated\n\nThis is a major change in `4.0`." | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "54a6051e", | |
| "cell_type": "code", | |
| "source": "import numpy as np\nimport aesara\nimport aesara.tensor as at\nimport theano\nimport theano.tensor as tt\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport pandas as pd\n\nnp.set_printoptions(3)\n\nRANDOM_SEED = 8927\nrng = np.random.default_rng(RANDOM_SEED)\n\n#aesara.config.mode = \"Mode\" # Default is \"Mode\"\n\n#aesara.config.mode = \"JAX\" # Default is \"Mode\"\n\nimport pymc3 as pm3\nimport pymc as pm4\n\n#aesara.config.mode = \"NUMBA\" # Default is \"Mode\"", | |
| "execution_count": 43, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "data = pd.read_csv(pm4.get_data(\"radon.csv\"))\ncounty_names = data.county.unique()\n\ndata[\"log_radon\"] = data[\"log_radon\"].astype(theano.config.floatX)\n\ncounty_idxs, counties = pd.factorize(data.county)\ncoords = {\n \"county\": counties,\n \"obs_id\": np.arange(len(county_idxs)),\n}", | |
| "execution_count": 10, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "def build_model(pm):\n with pm.Model(coords=coords) as hierarchical_model:\n mu_a = pm.Normal(\"mu_a\", mu=0.0, sigma=10)\n sigma_a = pm.HalfNormal(\"sigma_a\", 5.0)\n mu_b = pm.Normal(\"mu_b\", mu=0.0, sigma=10)\n sigma_b = pm.HalfNormal(\"sigma_b\", 5.0)\n a = pm.Normal(\"a\", dims=\"county\") * sigma_a + mu_a\n b = pm.Normal(\"b\", dims=\"county\") * sigma_b + mu_b\n eps = pm.HalfCauchy(\"eps\", 5.0)\n radon_est = a[county_idx] + b[county_idx] * data.floor.values\n radon_like = pm.Normal(\n \"radon_like\", mu=radon_est, sigma=eps, observed=data.log_radon, dims=\"obs_id\"\n )\n \n idata = pm.sample(return_inferencedata=True)\n \n return hierarchical_model, idata\n\nmodel_pymc3, idata_pymc3 = build_model(pm3)\nmodel_pymc4, idata_pymc4 = build_model(pm4)", | |
| "execution_count": 3, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (4 chains in 4 jobs)\nNUTS: [eps, b, a, sigma_b, mu_b, sigma_a, mu_a]\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n", | |
| "name": "stderr" | |
| }, | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "text/plain": "<IPython.core.display.HTML object>", | |
| "text/html": "\n <div>\n <style>\n /* Turns off some styling */\n progress {\n /* gets rid of default border in Firefox and Opera. */\n border: none;\n /* Needs to be in here for Safari polyfill so background images work as expected. */\n background-size: auto;\n }\n .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n background: #F44336;\n }\n </style>\n <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n 100.00% [8000/8000 00:05<00:00 Sampling 4 chains, 183 divergences]\n </div>\n " | |
| }, | |
| "metadata": {} | |
| }, | |
| { | |
| "output_type": "stream", | |
| "text": "/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\n/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\n/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\n/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\nSampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 12 seconds.\nThere were 9 divergences after tuning. Increase `target_accept` or reparameterize.\nThere were 46 divergences after tuning. Increase `target_accept` or reparameterize.\nThere were 128 divergences after tuning. Increase `target_accept` or reparameterize.\nThe acceptance probability does not match the target. It is 0.6798815133749979, but should be close to 0.8. Try to increase the number of tuning steps.\nThe number of effective samples is smaller than 10% for some parameters.\nAuto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\n/Users/twiecki/projects/pymc/pymc/model.py:984: FutureWarning: `Model.initial_point` has been deprecated. Use `Model.recompute_initial_point(seed=None)`.\n warnings.warn(\nMultiprocess sampling (4 chains in 4 jobs)\nNUTS: [mu_a, sigma_a, mu_b, sigma_b, a, b, eps]\n/Users/twiecki/projects/pymc/pymc/model.py:984: FutureWarning: `Model.initial_point` has been deprecated. Use `Model.recompute_initial_point(seed=None)`.\n warnings.warn(\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYour Python environment has Theano(-PyMC) 1.1.2 installed, but you are importing PyMC 4.0.0b1 which uses Aesara as its backend.\nFor PyMC 4.0.0b1 to work as expected you should uninstall Theano(-PyMC).\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYou are importing PyMC 4.0.0b1, but your environment also has the legacy version PyMC3 3.11.4 installed.\nFor PyMC 4.0.0b1 to work as expected you should uninstall PyMC3.\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYour Python environment has Theano(-PyMC) 1.1.2 installed, but you are importing PyMC 4.0.0b1 which uses Aesara as its backend.\nFor PyMC 4.0.0b1 to work as expected you should uninstall Theano(-PyMC).\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYou are importing PyMC 4.0.0b1, but your environment also has the legacy version PyMC3 3.11.4 installed.\nFor PyMC 4.0.0b1 to work as expected you should uninstall PyMC3.\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYour Python environment has Theano(-PyMC) 1.1.2 installed, but you are importing PyMC 4.0.0b1 which uses Aesara as its backend.\nFor PyMC 4.0.0b1 to work as expected you should uninstall Theano(-PyMC).\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYou are importing PyMC 4.0.0b1, but your environment also has the legacy version PyMC3 3.11.4 installed.\nFor PyMC 4.0.0b1 to work as expected you should uninstall PyMC3.\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nWARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYour Python environment has Theano(-PyMC) 1.1.2 installed, but you are importing PyMC 4.0.0b1 which uses Aesara as its backend.\nFor PyMC 4.0.0b1 to work as expected you should uninstall Theano(-PyMC).\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\nYou are importing PyMC 4.0.0b1, but your environment also has the legacy version PyMC3 3.11.4 installed.\nFor PyMC 4.0.0b1 to work as expected you should uninstall PyMC3.\nSee https://github.com/pymc-devs/pymc/wiki for update instructions.\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n", | |
| "name": "stderr" | |
| }, | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "text/plain": "<IPython.core.display.HTML object>", | |
| "text/html": "\n <div>\n <style>\n /* Turns off some styling */\n progress {\n /* gets rid of default border in Firefox and Opera. */\n border: none;\n /* Needs to be in here for Safari polyfill so background images work as expected. */\n background-size: auto;\n }\n .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n background: #F44336;\n }\n </style>\n <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n 100.00% [8000/8000 00:09<00:00 Sampling 4 chains, 4 divergences]\n </div>\n " | |
| }, | |
| "metadata": {} | |
| }, | |
| { | |
| "output_type": "stream", | |
| "text": "/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\n/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\n/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\n/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n return _boost._beta_ppf(q, a, b)\nSampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 19 seconds.\nThere were 4 divergences after tuning. Increase `target_accept` or reparameterize.\nThe number of effective samples is smaller than 25% for some parameters.\n", | |
| "name": "stderr" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "be9adfd8", | |
| "cell_type": "code", | |
| "source": "%%time\n\nwith model_pymc3:\n pm3.sample_posterior_predictive(idata_pymc3)", | |
| "execution_count": 8, | |
| "outputs": [ | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "text/plain": "<IPython.core.display.HTML object>", | |
| "text/html": "\n <div>\n <style>\n /* Turns off some styling */\n progress {\n /* gets rid of default border in Firefox and Opera. */\n border: none;\n /* Needs to be in here for Safari polyfill so background images work as expected. */\n background-size: auto;\n }\n .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n background: #F44336;\n }\n </style>\n <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n 100.00% [4000/4000 00:04<00:00]\n </div>\n " | |
| }, | |
| "metadata": {} | |
| }, | |
| { | |
| "output_type": "stream", | |
| "text": "CPU times: user 4.45 s, sys: 974 ms, total: 5.43 s\nWall time: 5.43 s\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "9ef50e85", | |
| "cell_type": "code", | |
| "source": "%%time\n\nwith model_pymc4:\n pm4.sample_posterior_predictive(idata_pymc4)", | |
| "execution_count": 9, | |
| "outputs": [ | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "text/plain": "<IPython.core.display.HTML object>", | |
| "text/html": "\n <div>\n <style>\n /* Turns off some styling */\n progress {\n /* gets rid of default border in Firefox and Opera. */\n border: none;\n /* Needs to be in here for Safari polyfill so background images work as expected. */\n background-size: auto;\n }\n .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n background: #F44336;\n }\n </style>\n <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n 100.00% [4000/4000 00:00<00:00]\n </div>\n " | |
| }, | |
| "metadata": {} | |
| }, | |
| { | |
| "output_type": "stream", | |
| "text": "CPU times: user 1.45 s, sys: 28.5 ms, total: 1.47 s\nWall time: 1.46 s\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## Random Variables behave like Tensors" | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### Clipping" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "f8c21f0a", | |
| "cell_type": "code", | |
| "source": "with pm3.Model():\n x3 = pm3.Normal(\"x\")\n \nwith pm4.Model():\n x4 = pm4.Normal(\"x\", shape=100)", | |
| "execution_count": 39, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "96bf2513", | |
| "cell_type": "code", | |
| "source": "type(x3)", | |
| "execution_count": 29, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "execution_count": 29, | |
| "data": { | |
| "text/plain": "pymc3.model.FreeRV" | |
| }, | |
| "metadata": {} | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "4974298a", | |
| "cell_type": "code", | |
| "source": "type(x4)", | |
| "execution_count": 30, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "execution_count": 30, | |
| "data": { | |
| "text/plain": "aesara.tensor.var.TensorVariable" | |
| }, | |
| "metadata": {} | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "cc309fae", | |
| "cell_type": "code", | |
| "source": "tt.clip(x3, 0, np.inf)", | |
| "execution_count": 31, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "execution_count": 31, | |
| "data": { | |
| "text/plain": "Elemwise{clip,no_inplace}.0" | |
| }, | |
| "metadata": {} | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "78504a5e", | |
| "cell_type": "code", | |
| "source": "sns.histplot(at.clip(x4, 0, np.inf).eval())", | |
| "execution_count": 47, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "execution_count": 47, | |
| "data": { | |
| "text/plain": "<AxesSubplot:ylabel='Count'>" | |
| }, | |
| "metadata": {} | |
| }, | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "text/plain": "<Figure size 432x288 with 1 Axes>", | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD4CAYAAADrRI2NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAOSklEQVR4nO3dfayedX3H8feHtoiZOCU9dF1pV43EyJZMzNEhLAvKWJhbVmYAMRs2C1vdg4vMxdm5ZP7bJYtxD2baqLFujMF8GJWhjlUeYmSMU4YK6VydUeho2ooPQFycxe/+OBdyaHvOuTg913333L/3K7lzX4/n+v74hU+v/O7r/t2pKiRJ7Tht3AVIkkbL4Jekxhj8ktQYg1+SGmPwS1JjVo+7gD7Wrl1bmzdvHncZkrSi7N279xtVNXXs9hUR/Js3b2ZmZmbcZUjSipLk6yfa7lCPJDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1ZuKDf8PGTSQZy2vDxk3jbr4kHWdFTNlwMh458DBveP/nx3LtG9984ViuK0kLmfg7fknSMxn8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMG/bH1JF8DHgeeBI5W1XSSs4Abgc3A14CrqupbQ9YhSXraKO74X1NVL6+q6W59O7Cnqs4F9nTrkqQRGcdQzxZgV7e8C7h8DDVIUrOGDv4C/iXJ3iTbum3rquogQPd+9olOTLItyUySmSNHjgxcpiS1Y9AxfuCiqnokydnAbUn+s++JVbUT2AkwPT1dQxUoSa0Z9I6/qh7p3g8DnwBeBRxKsh6gez88ZA2SpGcaLPiT/EiSM59aBn4BeADYDWztDtsK3DxUDZKk4w051LMO+ESSp67z91X16ST3AjcluRZ4CLhywBokSccYLPir6qvAT59g+6PAJUNdV5K0ML+5K0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1ZvDgT7IqyX8kuaVbPyvJbUn2d+8vHLoGSdLTRnHH/1Zg35z17cCeqjoX2NOtS5JGZNDgT3IO8EvAB+Zs3gLs6pZ3AZcPWYMk6ZmGvuN/D/BHwA/mbFtXVQcBuvezT3Rikm1JZpLMHDlyZOAyJakdgwV/kl8GDlfV3qWcX1U7q2q6qqanpqaWuTpJatfqAf/2RcCvJHkdcAbw/CR/BxxKsr6qDiZZDxwesAZJ0jEGu+Ovqj+uqnOqajNwNfDZqvp1YDewtTtsK3DzUDVIko43juf4dwCXJtkPXNqtS5JGZMihnh+qqjuAO7rlR4FLRnFdSdLx/OauJDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1JhewZ/koj7bJEmnvr53/H/Vc5sk6RS3eqGdSV4NXAhMJXnbnF3PB1YNWZgkaRgLBj9wOvC87rgz52x/DLhiqKIkScNZMPir6k7gziQfrqqvj6gmSdKAFrvjf8pzkuwENs89p6peO0RRkqTh9A3+fwTeB3wAeHK4ciRJQ+sb/Eer6m8GrUSSNBJ9H+f8ZJLfTbI+yVlPvQatTJI0iL53/Fu797fP2VbAi+c7IckZwF3Ac7rrfLSq3tX9g3Ejs58XfA24qqq+9ezKliQtVa/gr6oXLeFvfw94bVU9kWQN8LkknwJeD+ypqh1JtgPbgXcs4e9LkpagV/AnedOJtlfVR+Y7p6oKeKJbXdO9CtgCXNxt3wXcgcEvSSPTd6jnlXOWzwAuAe4D5g1+gCSrgL3AS4D3VtU9SdZV1UGAqjqY5Ox5zt0GbAPYtGlTzzIlSYvpO9Tz+3PXk/wo8Lc9znsSeHmSFwCfSPJTfQurqp3AToDp6enqe54kaWFLnZb5u8C5fQ+uqm8zO6RzGXAoyXqA7v3wEmuQJC1B3zH+TzI7Pg+zk7O9DLhpkXOmgO9X1beTPBf4eeDPgN3MPiW0o3u/eWmlS5KWou8Y/5/PWT4KfL2qDixyznpgVzfOfxpwU1XdkuRu4KYk1wIPAVc+26IlSUvXd4z/ziTrePpD3v09zvkicP4Jtj/K7IfDkqQx6PsLXFcB/87s3flVwD1JnJZZklagvkM9fwK8sqoOww/H7/8V+OhQhUmShtH3qZ7Tngr9zqPP4lxJ0imk7x3/p5N8BrihW38DcOswJUmShrTYb+6+BFhXVW9P8nrgZ4EAdwPXj6A+SdIyW2y45j3A4wBV9fGqeltV/QGzd/vvGbY0SdIQFgv+zd1jmc9QVTPMTqssSVphFgv+MxbY99zlLESSNBqLBf+9SX7r2I3dt273DlOSJGlIiz3Vcx2zs2r+Gk8H/TRwOvCrA9YlSRrIgsFfVYeAC5O8BnhqSuV/rqrPDl6ZJGkQfefquR24feBaJEkj4LdvJakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JasxgwZ9kY5Lbk+xL8mCSt3bbz0pyW5L93fsLh6pBknS8Ie/4jwJ/WFUvAy4Afi/JecB2YE9VnQvs6dYlSSMyWPBX1cGquq9bfhzYB2wAtgC7usN2AZcPVYMk6XgjGeNPshk4H7gHWFdVB2H2Hwfg7HnO2ZZkJsnMkSNHRlGmJDVh8OBP8jzgY8B1VfVY3/OqamdVTVfV9NTU1HAFSlJjBg3+JGuYDf3rq+rj3eZDSdZ3+9cDh4esQZL0TEM+1RPgg8C+qnr3nF27ga3d8lbg5qFqkCQdb/WAf/si4BrgS0nu77a9E9gB3JTkWuAh4MoBa5AkHWOw4K+qzwGZZ/clQ11XkrQwv7krSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMM/gm0YeMmkoz8tWHjpnE3XVIPq8ddgJbfIwce5g3v//zIr3vjmy8c+TUlPXve8UtSYwx+SWrMYMGf5ENJDid5YM62s5LclmR/9/7Coa5/Sjht9VjG2iVpIUOO8X8Y+GvgI3O2bQf2VNWOJNu79XcMWMN4/eCoY+2STjmD3fFX1V3AN4/ZvAXY1S3vAi4f6vqSpBMb9Rj/uqo6CNC9nz3fgUm2JZlJMnPkyJGRFShJk+6U/XC3qnZW1XRVTU9NTY27HEmaGKMO/kNJ1gN074dHfH1Jat6og383sLVb3grcPOLrS1Lzhnyc8wbgbuClSQ4kuRbYAVyaZD9wabcuSRqhwR7nrKo3zrPrkqGuKUla3Cn74a4kaRgGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8mwoaNm8bya2cbNm4ad9OlZ23IX+CSRuaRAw/7a2dST97xS1JjDH5JaozBL0mNMfglqTEGvyQ1xqd6pJNx2mqSjOXSq9Y8hye//72RX/fHz9nI/zz80Mivq+Vj8Esn4wdHx/IYKcw+SuojrFoKh3okqTEGvyQ1xqEeLZ8xjndL6s/g1/IZ83i3pH4c6pGkxhj8ktQYg1+SGmPwS1JjDH5JaoxP9Uh6dhqcpmJc14VhpsgYS/AnuQz4C2AV8IGq2jGOOiQtQaPTVEzSo8ojH+pJsgp4L/CLwHnAG5OcN+o6JKlV4xjjfxXwlar6alX9H/APwJYx1CFJTUpVjfaCyRXAZVX1m936NcDPVNVbjjluG7CtW30p8OUlXnIt8I0lnrtS2eY22OY2nEybf6Kqpo7dOI4x/hN9KnTcvz5VtRPYedIXS2aqavpk/85KYpvbYJvbMESbxzHUcwDYOGf9HOCRMdQhSU0aR/DfC5yb5EVJTgeuBnaPoQ5JatLIh3qq6miStwCfYfZxzg9V1YMDXvKkh4tWINvcBtvchmVv88g/3JUkjZdTNkhSYwx+SWrMxAR/ksuSfDnJV5JsP8H+JPnLbv8Xk7xiHHUupx5tvjjJd5Lc373+dBx1LpckH0pyOMkD8+yfxD5erM0T1ccASTYmuT3JviQPJnnrCY6ZqL7u2ebl6+uqWvEvZj8k/m/gxcDpwBeA84455nXAp5j9HsEFwD3jrnsEbb4YuGXctS5jm38OeAXwwDz7J6qPe7Z5ovq4a9N64BXd8pnAfzXw/3OfNi9bX0/KHX+faSC2AB+pWf8GvCDJ+lEXuoyam/qiqu4CvrnAIZPWx33aPHGq6mBV3dctPw7sAzYcc9hE9XXPNi+bSQn+DcDDc9YPcPx/tD7HrCR92/PqJF9I8qkkPzma0sZm0vq4r4nt4ySbgfOBe47ZNbF9vUCbYZn6elLm4+8zDUSvqSJWkD7tuY/ZuTqeSPI64J+Ac4cubIwmrY/7mNg+TvI84GPAdVX12LG7T3DKiu/rRdq8bH09KXf8faaBmLSpIhZtT1U9VlVPdMu3AmuSrB1diSM3aX28qEnt4yRrmA3A66vq4yc4ZOL6erE2L2dfT0rw95kGYjfwpu5pgAuA71TVwVEXuowWbXOSH0v3U0lJXsVsfz868kpHZ9L6eFGT2Mddez4I7Kuqd89z2ET1dZ82L2dfT8RQT80zDUSS3+72vw+4ldknAb4CfBf4jXHVuxx6tvkK4HeSHAX+F7i6uscDVqIkNzD7ZMPaJAeAdwFrYDL7GHq1eaL6uHMRcA3wpST3d9veCWyCie3rPm1etr52ygZJasykDPVIknoy+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1Jj/h+X4vS93g1mHAAAAABJRU5ErkJggg==\n" | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| } | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### Stacking & Indexing" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "6a42b3b7", | |
| "cell_type": "code", | |
| "source": "with pm4.Model():\n x = pm4.Uniform(\"x\", lower=-1, upper=0)\n y = pm4.Uniform(\"y\", lower=0, upper=1)\n xy = at.stack([x, y])\n index = pm4.Bernoulli(\"index\", p=0.5)\n\nfor _ in range(5):\n print(\"Sampled value = {:.2f}\".format(xy[index].eval()))", | |
| "execution_count": 14, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "Sampled value = -0.63\nSampled value = -0.79\nSampled value = 0.06\nSampled value = -0.47\nSampled value = -0.94\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### For loops with scan" | |
| }, | |
| { | |
| "metadata": {}, | |
| "id": "374a426b", | |
| "cell_type": "markdown", | |
| "source": "## Dynamic Shape Support" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "id": "0daf28f1", | |
| "cell_type": "code", | |
| "source": "with pm4.Model() as m:\n x = pm4.Poisson('x', 3)\n z = pm4.Normal('z', shape=x)\n \nfor _ in range(5):\n print(\"Sampled value z = {}\".format(z.eval()))", | |
| "execution_count": 19, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "Sampled value z = [0.792]\nSampled value z = [-0.536 0.037 -0.445]\nSampled value z = [ 1.316 0.373 -0.017 0.108 0.045]\nSampled value z = [-0.86]\nSampled value z = [ 1.085 1.435 -0.225 0.085]\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## Conditioning & Causal Inference" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# 1. define joint P(mu, sigma, x)\nwith pm4.Model() as m:\n mu = pm4.Normal(\"mu\", 0, 1)\n sigma = pm4.HalfNormal(\"sigma\", 1)\n x = pm4.Normal(\"obs\", mu, sigma)\n\n# 2. Generate N samples from P(x|mu=0, sigma=0.1) \nN = 10\ncond_x = aesara.function([mu, sigma], x)\n_x = np.array([cond_x(0, 0.1) for _ in range(N)])\n\nprint(\"10 samples from P(x | mu=0, sigma=0.1) = {}\".format(_x))", | |
| "execution_count": 27, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "text": "10 samples from P(x | mu=0, sigma=0.1) = [-0.005 -0.238 0.031 0.045 0.136 -0.112 -0.232 -0.077 -0.033 0.004]\n", | |
| "name": "stdout" | |
| } | |
| ] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## Better NUTS initialization" | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## What's not yet working\n\nCurrently we are still missing a few features from `v3`. These should be hard to add and will be there for the `4.0` stable release. The list includes:\n\n* Timeseries distributions\n* Mixture Distributions" | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "name": "pymc4b1", | |
| "display_name": "pymc4b1", | |
| "language": "python" | |
| }, | |
| "language_info": { | |
| "name": "python", | |
| "version": "3.8.12", | |
| "mimetype": "text/x-python", | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "pygments_lexer": "ipython3", | |
| "nbconvert_exporter": "python", | |
| "file_extension": ".py" | |
| }, | |
| "gist": { | |
| "id": "", | |
| "data": { | |
| "description": "PyMC3 vs PyMC 4.0.ipynb", | |
| "public": true | |
| } | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 5 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment