Skip to content

Instantly share code, notes, and snippets.

@twiecki

twiecki/PyMC3 vs PyMC 4.0.ipynb

Last active January 6, 2022 16:01
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save twiecki/badf23b78fd1237d16f4c1909ddbe3e3 to your computer and use it in GitHub Desktop.
Save twiecki/badf23b78fd1237d16f4c1909ddbe3e3 to your computer and use it in GitHub Desktop.
Download ZIP
PyMC3 vs PyMC 4.0.ipynb
Display the source blob
Display the rendered blob
Raw
PyMC3 vs PyMC 4.0.ipynb
{
"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. PyMC3 version 3.x will stay under the current name to not break production systems but future versions will use PyMC everywhere.** While there were 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:\n1. that TF was really difficult to work with, partly because of its eager execution mode, but, more importantly, and\n2. 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, and our main focus going forward. The version that 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. You can just forget about that one.\n\nIn general, you should refer to the new version then as \"PyMC 4.0\".\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 in Greek mythology), and is now a project largely separate from `PyMC`, although `PyMC` is the main customer and there is considerable overlap in the development teams.\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). With this graph representation, we can do various things:\n* graph optimizations like `log(exp(x)) -> x`\n* symbolic rewrites like `N(0, 1) + a` -> `N(a, 1)`\n* compilation of that graph to various computational backends.\n\nPreviously, `theano` supported Python and C as computational backends. But with `aesara` it is now possible, and in fact quite easy, to add new computational backends. We have currently added a `JAX` backend that comes with GPU support (see [this blog post](https://martiningram.github.io/mcmc-comparison/) for some impressive speed-ups using GPUs for sampling). We're also in the process of adding a `numba` backend.\n\n## What's new\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."
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## Using JAX for faster sampling\n\nBy far the most exciting new feature are the new computational backends and the associated speed-ups. While we believe that the `numba` backend shows a lot of promise, the JAX backend is currently more mature and can already be used for NUTS sampling.\n\nAs mentioned above, `aesara` provides a representation of the model logp graph in form of various `aesara` `Ops` (operators) which represent the computations to be be performed. For example `exp(x + y)` would be an `Add` `Op` with two input arguments `x` and `y`, which are inputted into an `exp` `Op`.\n\nThis doesn't say anything about how we actually *execute* this graph. Before, we would transpile this graph to C-code which could then compile, load into Python as a C-extension, and then execute. But now, we can just transpile this graph to `JAX` (or `numba`, or whatever) instead.\n\nWhile this by itself is already pretty exciting, because `JAX` is capable of a whole bunch of low-level optimizations which lead to faster logp functions, we still have our samplers written in Python, so there is still some call-overhead and things we can't optimize. \n\nTo get rid of this, we can link this up with samplers also written in `JAX`. That way, the model logp evaluation *and* the sampler are one big JAX graph that gets optimized and executed without any Python call-overhead. We currently support a NUTS implementation provided by `numpyro` as well as `blackjax`.\n\nWe have not yet done extensive benchmarks, but what we have so far is very promising. Here is a small example of how much faster this is already on a fairly small and simple model: the hierarchical linear regression of the famous Radon example:"
},
{
"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 arviz as az\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport pandas as pd\nRANDOM_SEED = 8927\nrng = np.random.default_rng(RANDOM_SEED)\nnp.set_printoptions(2)\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": 2,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Load in data:"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "data = pd.read_csv(\"/Users/twiecki/personal/homepage/WhileMyMCMCGentlySamples/content/radon.csv\")",
"execution_count": 8,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"id": "eea8dba8",
"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_idx, counties = pd.factorize(data.county)\ncoords = {\n \"county\": counties,\n \"obs_id\": np.arange(len(county_idx)),\n}",
"execution_count": 10,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Define a function to create the model. Note that we provide `pm`, our PyMC library, as an argument here. This allows us to create this model in `pymc3` or `pymc` 4.0. Here you can also see that most models that work in `pymc3` also work in `pymc` 4.0 without any code change, you only need to change your imports."
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "def build_model(pm, **kwargs):\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.HalfNormal(\"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 return hierarchical_model",
"execution_count": 29,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Create and sample model in `pymc3`, nothing special:"
},
{
"metadata": {
"trusted": true,
"scrolled": true
},
"cell_type": "code",
"source": "%%time\nmodel_pymc3 = build_model(pm3)\nwith model_pymc3:\n idata_pymc3 = pm3.sample(return_inferencedata=True)",
"execution_count": 32,
"outputs": [
{
"output_type": "stream",
"text": "Auto-assigning NUTS sampler...\nINFO:pymc3:Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nINFO:pymc3:Initializing NUTS using jitter+adapt_diag...\nMultiprocess sampling (4 chains in 4 jobs)\nINFO:pymc3:Multiprocess sampling (4 chains in 4 jobs)\nNUTS: [eps, b, a, sigma_b, mu_b, sigma_a, mu_a]\nINFO:pymc3:NUTS: [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:04<00:00 Sampling 4 chains, 58 divergences]\n </div>\n "
},
"metadata": {}
},
{
"output_type": "stream",
"text": "/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.9/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.9/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.9/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.9/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 10 seconds.\nINFO:pymc3:Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 10 seconds.\nThere were 50 divergences after tuning. Increase `target_accept` or reparameterize.\nERROR:pymc3:There were 50 divergences after tuning. Increase `target_accept` or reparameterize.\nThere were 6 divergences after tuning. Increase `target_accept` or reparameterize.\nERROR:pymc3:There were 6 divergences after tuning. Increase `target_accept` or reparameterize.\nThere were 2 divergences after tuning. Increase `target_accept` or reparameterize.\nERROR:pymc3:There were 2 divergences after tuning. Increase `target_accept` or reparameterize.\nThe number of effective samples is smaller than 25% for some parameters.\nINFO:pymc3:The number of effective samples is smaller than 25% for some parameters.\n",
"name": "stderr"
},
{
"output_type": "stream",
"text": "CPU times: user 4.23 s, sys: 244 ms, total: 4.47 s\nWall time: 12.7 s\n",
"name": "stdout"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Create and sample model in `pymc` 4.0, also nothing special:"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%time\nmodel_pymc4 = build_model(pm4)\nwith model_pymc4:\n idata_pymc4 = pm4.sample()",
"execution_count": 33,
"outputs": [
{
"output_type": "stream",
"text": "Auto-assigning NUTS sampler...\nINFO:pymc:Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag...\nINFO:pymc:Initializing 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)\nINFO:pymc:Multiprocess sampling (4 chains in 4 jobs)\nNUTS: [mu_a, sigma_a, mu_b, sigma_b, a, b, eps]\nINFO:pymc:NUTS: [mu_a, sigma_a, mu_b, sigma_b, a, b, eps]\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:06<00:00 Sampling 4 chains, 204 divergences]\n </div>\n "
},
"metadata": {}
},
{
"output_type": "stream",
"text": "/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.9/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.9/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.9/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.9/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 14 seconds.\nINFO:pymc:Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 14 seconds.\nThere were 147 divergences after tuning. Increase `target_accept` or reparameterize.\nERROR:pymc:There were 147 divergences after tuning. Increase `target_accept` or reparameterize.\nThe acceptance probability does not match the target. It is 0.673, but should be close to 0.8. Try to increase the number of tuning steps.\nWARNING:pymc:The acceptance probability does not match the target. It is 0.673, but should be close to 0.8. Try to increase the number of tuning steps.\nThere were 55 divergences after tuning. Increase `target_accept` or reparameterize.\nERROR:pymc:There were 55 divergences after tuning. Increase `target_accept` or reparameterize.\nThe acceptance probability does not match the target. It is 0.7123, but should be close to 0.8. Try to increase the number of tuning steps.\nWARNING:pymc:The acceptance probability does not match the target. It is 0.7123, but should be close to 0.8. Try to increase the number of tuning steps.\nThere were 2 divergences after tuning. Increase `target_accept` or reparameterize.\nERROR:pymc:There were 2 divergences after tuning. Increase `target_accept` or reparameterize.\nThe number of effective samples is smaller than 10% for some parameters.\nWARNING:pymc:The number of effective samples is smaller than 10% for some parameters.\n",
"name": "stderr"
},
{
"output_type": "stream",
"text": "CPU times: user 4.54 s, sys: 254 ms, total: 4.79 s\nWall time: 16.4 s\n",
"name": "stdout"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "Now, lets use a JAX sampler instead. Here we use the one provided by `numpyro`:"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import pymc.sampling_jax",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true,
"scrolled": true
},
"cell_type": "code",
"source": "%%time\nwith model_pymc4:\n idata = pm4.sampling_jax.sample_numpyro_nuts(progress_bar=False)",
"execution_count": 35,
"outputs": [
{
"output_type": "stream",
"text": "Compiling...\n",
"name": "stdout"
},
{
"output_type": "stream",
"text": "/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(\n",
"name": "stderr"
},
{
"output_type": "stream",
"text": "Compilation time = 0 days 00:00:00.267406\nSampling...\nSampling time = 0 days 00:00:05.143948\nTransforming variables...\nTransformation time = 0 days 00:00:00.017413\nCPU times: user 7.99 s, sys: 47.2 ms, total: 8.04 s\nWall time: 5.53 s\n",
"name": "stdout"
},
{
"output_type": "stream",
"text": "/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.9/site-packages/aesara/graph/opt.py:279: UserWarning: Supervisor is not added. Please build a FunctionGraphvia aesara.compile.function.types.std_graph()or add the Supervisor class manually.\n sub_prof = optimizer.optimize(fgraph)\n/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.9/site-packages/aesara/graph/opt.py:279: UserWarning: Supervisor is not added. Please build a FunctionGraphvia aesara.compile.function.types.std_graph()or add the Supervisor class manually.\n sub_prof = optimizer.optimize(fgraph)\n",
"name": "stderr"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "16.4 / 5.53",
"execution_count": 36,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 36,
"data": {
"text/plain": "2.9656419529837246"
},
"metadata": {}
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "That's a speed-up of 3x -- for a single-line code change! And this is just running things on the CPU, we can just as easily run this on the GPU where we saw even more impressive speed-ups (especially as we scale the data).\n\nAgain, for a more proper benchmark that also compares this to Stan, see [this blog post](https://martiningram.github.io/mcmc-comparison/).\n\n### The next step: Samplers written in `aesara`\n\nWhile this current approach is already quite exciting, we can take this one step further. The setup we showed above has the model logp graph (represented in `aesara` compiled to `JAX`. This `JAX` function we then linked to a sampler written in `JAX`.\n\nThis is suboptimal for two reasons:\n1. For new backends, like `numba`, we need to rewrite the sampler also in `numba`.\n2. While we get low-level optimizations from `JAX` on the logp+sampler JAX-graph, we do not get any high-level optimizations, which is what `aesara` is great at.\n\nWith [`aehmc`](https://www.github.com/aesara-devs/aehmc) and [`aemcmc`](https://www.github.com/aesara-devs/aemcmc) the `aesara` devs are developing a library of samplers *written in `aesara`*. That way, our model logp, consisting out of `aesara` `Ops` can be combined with the sampler logic, also consisting out of `aesara` `Ops`, into one big `aesara` graph.\n\nOn that graph, `aesara` can do high-level optimizations to get a more efficient graph representation, and then compile it to whatever backend we want: `JAX`, `numba`, `C`, or whatever other backend we add in the future.\n\nIf you think this is interesting, definitely check out these packages and consider contributing, this is where the next round of innovation will come from!"
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## Better integration into `aesara`\n\nThe next feature we are excited about is a better integration of `PyMC` and `aesara`.\n\nIn `PyMC3 3.x`, the random variables (RVs) created by e.g. calling `x = pm.Normal('x')` were not truly `theano` `Ops` which you would integrate nicely with the rest of `theano`. This created a lot of issues, limitations, and complexities in the library.\n\n`Aesara` now provides a proper `RandomVariable` Op which perfectly integrates with the rest of the other `Ops`. \n\nThis is a major change in `4.0` and lead to huge swaths of code getting removed and greatly simplified. In many ways, this change is much more exciting than the different computational backends, but the effects are not quite as visible to the user.\n\nThere are a few cases, however, where you can see the benefits."
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### Faster posterior predictive sampling"
},
{
"metadata": {
"trusted": false
},
"id": "be9adfd8",
"cell_type": "code",
"source": "%%time\n\nwith model_pymc3:\n pm3.sample_posterior_predictive(idata_pymc3)",
"execution_count": 86,
"outputs": [
{
"data": {
"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 ",
"text/plain": "<IPython.core.display.HTML object>"
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": "CPU times: user 5.04 s, sys: 1.12 s, total: 6.16 s\nWall time: 6.17 s\n"
}
]
},
{
"metadata": {
"trusted": false
},
"id": "9ef50e85",
"cell_type": "code",
"source": "%%time\n\nwith model_pymc4:\n pm4.sample_posterior_predictive(idata_pymc4)",
"execution_count": 87,
"outputs": [
{
"data": {
"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 ",
"text/plain": "<IPython.core.display.HTML object>"
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": "CPU times: user 1.59 s, sys: 32 ms, total: 1.62 s\nWall time: 1.6 s\n"
}
]
},
{
"metadata": {},
"id": "15599440",
"cell_type": "markdown",
"source": "On this model, we get a speed-up of almost 4x. This will vary on the model but in general, on models where predictive sampling was already fast, it will probably not be much faster, but on models where predictive sampling was really slow, you can expect a big speed-up.\n\nThe reason for this is that predictive sampling is now happening as part of the `aesara` graph. Before, we were walking through the random variables in Python which was not only slow, but also very error-prone, so a lot of dev time was spent fixing bugs and rewriting this complicated piece of code. In `PyMC` 4.0, all that complexity is gone."
},
{
"metadata": {},
"id": "3d7ea093",
"cell_type": "markdown",
"source": "## Work with RVs just like with Tensors\n\nIn PyMC3, RVs as returned by e.g. `pm.Normal(\"x\")` behaved somewhat like a Tensor variable, but not *quite*. In PyMC 4, RVs are first-class Tensor variables that can be operated on much more freely."
},
{
"metadata": {
"trusted": false
},
"id": "f8c21f0a",
"cell_type": "code",
"source": "with pm3.Model():\n x3 = pm3.Normal(\"x\")\n \nwith pm4.Model():\n x4 = pm4.Normal(\"x\")",
"execution_count": 13,
"outputs": []
},
{
"metadata": {
"trusted": false
},
"id": "96bf2513",
"cell_type": "code",
"source": "type(x3)",
"execution_count": 14,
"outputs": [
{
"data": {
"text/plain": "pymc3.model.FreeRV"
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": false
},
"id": "4974298a",
"cell_type": "code",
"source": "type(x4)",
"execution_count": 15,
"outputs": [
{
"data": {
"text/plain": "aesara.tensor.var.TensorVariable"
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {},
"id": "db22a81f",
"cell_type": "markdown",
"source": "Through the power of [`aeppl`](https://github.com/aesara-devs/aeppl) (a new low-level library that provides core building blocks for probabilistic programming languages on top of `aesara`), PyMC 4.0 allows you to do even more operations directly on the RV as was possible before.\n\nFor example, we can `clip()` a RV using `aesara`'s `clip()` function. Also, calling `.eval()` on a RV samples a random draw from the RV, this is also new in PyMC 4.0 and makes things more consistent and allows easy interactions. "
},
{
"metadata": {
"trusted": false
},
"id": "78504a5e",
"cell_type": "code",
"source": "at.clip(x4, 0, np.inf).eval()",
"execution_count": 24,
"outputs": [
{
"data": {
"text/plain": "array(1.39867239)"
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {
"trusted": false
},
"id": "eea963c2",
"cell_type": "code",
"source": "trunc_norm = [at.clip(x4, 0, np.inf).eval() for _ in range(1000)]\nsns.histplot(np.asarray(trunc_norm))",
"execution_count": 36,
"outputs": [
{
"data": {
"text/plain": "<AxesSubplot:ylabel='Count'>"
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAARgUlEQVR4nO3da6wdV3nG8f8T54YKlERxUuNLHYRbkSBxkUkDqRAQ1LiACEUEjApYVVpTGiooFZBQqYgPlvKhQrQIChZFmHIJboHGpFwaDAEhIMHh7oQ0LoHk1FZsgiAgqlCHtx/OpOx4ncv4Mmd27P9P2toza9bs/Wbl2I/XrL3npKqQJGnSSWMXIEmaPoaDJKlhOEiSGoaDJKlhOEiSGiePXcDROOuss2rt2rVjlyFJDyk333zzj6pq+UJ9HtLhsHbtWnbt2jV2GZL0kJLkh4v18bKSJKlhOEiSGoaDJKlhOEiSGoaDJKlhOEiSGoaDJKlhOEiSGoaDJKlxQofDytVrSLLoY+XqNWOXKklL6iF9+4yjtXfmLl7y7i8v2u8jr3zaElQjSdPjhJ45SJLmZjhIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpYThIkhqGgySpMWg4JPlBku8k+WaSXV3bmUmuT3J793zGRP+rkuxJcluSS4asTZI0v6WYOTyzqp5YVeu7/SuBnVW1DtjZ7ZPkPGAjcD6wAXhnkmVLUJ8k6RBjXFa6FNjWbW8DXjDRfk1V3VdVdwB7gAuWvjxJ0tDhUMB/JLk5yeau7Zyq2gfQPZ/dta8E7po4d6Zre5Akm5PsSrLrwIEDA5YuSSeukwd+/Yuqam+Ss4Hrk3xvgb6Zo62ahqqtwFaA9evXN8clSUdv0JlDVe3tnvcDH2f2MtHdSVYAdM/7u+4zwOqJ01cBe4esT5I0t8HCIclvJHnEA9vAHwDfBXYAm7pum4Bru+0dwMYkpyU5F1gH3DRUfZKk+Q15Wekc4ONJHnifD1XVp5N8Ddie5HLgTuAygKranWQ7cAtwELiiqu4fsD5J0jwGC4eq+j7whDna7wEunuecLcCWoWqSJPXjN6QlSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUMBwkSQ3DQZLUGDwckixL8o0k13X7Zya5Psnt3fMZE32vSrInyW1JLhm6NknS3JZi5vAa4NaJ/SuBnVW1DtjZ7ZPkPGAjcD6wAXhnkmVLUJ8k6RCDhkOSVcBzgfdMNF8KbOu2twEvmGi/pqruq6o7gD3ABUPWJ0ma29Azh7cBbwB+NdF2TlXtA+iez+7aVwJ3TfSb6doeJMnmJLuS7Dpw4MAgRUvSiW6wcEjyPGB/Vd3c95Q52qppqNpaVeurav3y5cuPqkZJ0txOHvC1LwKen+Q5wOnAI5N8ALg7yYqq2pdkBbC/6z8DrJ44fxWwd8D6JEnzGGzmUFVXVdWqqlrL7ELz56rqZcAOYFPXbRNwbbe9A9iY5LQk5wLrgJuGqk+SNL8hZw7zuRrYnuRy4E7gMoCq2p1kO3ALcBC4oqruH6E+STrhLUk4VNUNwA3d9j3AxfP02wJsWYqaJEnz8xvSkqSG4SBJahgOkqSG4SBJahgOkqSG4SBJahgOkqSG4SBJahgOkqSG4SBJahgOkqSG4SBJahgOkqSG4SBJahgOkqSG4SBJahgOkqSG4SBJahgOkqRGr3BIclGfNknS8aHvzOHtPdskSceBkxc6mOSpwNOA5UleN3HokcCyIQuTJI1nwXAATgUe3vV7xET7vcCLhipKkjSuBcOhqr4AfCHJ+6rqh0tUkyRpZIvNHB5wWpKtwNrJc6rqWUMUJUkaV99w+BfgXcB7gPuHK0eSNA36hsPBqvrHQSuRJE2Nvh9l/USSv0iyIsmZDzwWOiHJ6UluSvKtJLuTvKVrPzPJ9Ulu757PmDjnqiR7ktyW5JKj+O+SJB2FvjOHTd3z6yfaCnjMAufcBzyrqn6e5BTgS0k+BbwQ2FlVVye5ErgSeGOS84CNwPnAo4HPJvmdqvIyliQtsV7hUFXnHu4LV1UBP+92T+keBVwKPKNr3wbcALyxa7+mqu4D7kiyB7gA+Mrhvrck6ej0Cockr5irvarev8h5y4CbgccC76iqG5OcU1X7uvP3JTm7674S+OrE6TNd26GvuRnYDLBmzZo+5UuSDlPfy0pPmdg+HbgY+DqwYDh0l4SemORRwMeTPH6B7pnrJeZ4za3AVoD169c3xyVJR6/vZaW/nNxP8pvAP/d9k6r6SZIbgA3A3UlWdLOGFcD+rtsMsHritFXA3r7vIUk6do70lt2/ANYt1CHJ8m7GQJKHAc8Gvgfs4NcL3JuAa7vtHcDGJKclObd7/ZuOsD5J0lHou+bwCX59iWcZ8Dhg+yKnrQC2desOJwHbq+q6JF8Btie5HLgTuAygqnYn2Q7cAhwErvCTSpI0jr5rDn83sX0Q+GFVzSx0QlV9G3jSHO33MLtmMdc5W4AtPWuSJA2k12Wl7gZ832P2zqxnAL8csihJ0rj6/ia4FzN7/f8y4MXAjUm8ZbckHaf6Xlb6G+ApVbUfZhebgc8C/zpUYZKk8fT9tNJJDwRD557DOFeS9BDTd+bw6SSfAT7c7b8E+OQwJUmSxrbY75B+LHBOVb0+yQuB32f2m8xfAT64BPVJkkaw2KWhtwE/A6iqj1XV66rqr5idNbxt2NIkSWNZLBzWdt9XeJCq2sXsrwyVJB2HFguH0xc49rBjWYgkaXosFg5fS/JnhzZ2t764eZiSJEljW+zTSq9l9lbbf8yvw2A9cCrwRwPWJUka0YLhUFV3A09L8kzggd/F8O9V9bnBK5Mkjabv73P4PPD5gWuRJE0Jv+UsSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoMFg5JVif5fJJbk+xO8pqu/cwk1ye5vXs+Y+Kcq5LsSXJbkkuGqk2StLAhZw4Hgb+uqscBFwJXJDkPuBLYWVXrgJ3dPt2xjcD5wAbgnUmWDVifJGkeg4VDVe2rqq932z8DbgVWApcC27pu24AXdNuXAtdU1X1VdQewB7hgqPokSfNbkjWHJGuBJwE3AudU1T6YDRDg7K7bSuCuidNmurZDX2tzkl1Jdh04cGDQuiXpRDV4OCR5OPBR4LVVde9CXedoq6ahamtVra+q9cuXLz9WZUqSJgwaDklOYTYYPlhVH+ua706yoju+Atjftc8AqydOXwXsHbI+SdLchvy0UoB/Am6tqrdOHNoBbOq2NwHXTrRvTHJaknOBdcBNQ9UnSZrfyQO+9kXAy4HvJPlm1/Ym4Gpge5LLgTuBywCqaneS7cAtzH7S6Yqqun/A+iRJ8xgsHKrqS8y9jgBw8TznbAG2DFWTJKkfvyEtSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoMFg5J3ptkf5LvTrSdmeT6JLd3z2dMHLsqyZ4ktyW5ZKi6JEmLG3Lm8D5gwyFtVwI7q2odsLPbJ8l5wEbg/O6cdyZZNmBtkqQFDBYOVfVF4MeHNF8KbOu2twEvmGi/pqruq6o7gD3ABUPVJkla2FKvOZxTVfsAuuezu/aVwF0T/Wa6tkaSzUl2Jdl14MCBQYuVpBPVtCxIZ462mqtjVW2tqvVVtX758uUDlyVJJ6alDoe7k6wA6J73d+0zwOqJfquAvUtcmySps9ThsAPY1G1vAq6daN+Y5LQk5wLrgJuWuDZJUufkoV44yYeBZwBnJZkB3gxcDWxPcjlwJ3AZQFXtTrIduAU4CFxRVfcPVZskaWGDhUNVvXSeQxfP038LsGWoeiRJ/U3LgrQkaYoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDpKkhuEgSWoYDsfQytVrSLLoY+XqNWOXKkkLGuz2GSeivTN38ZJ3f3nRfh955dOWoBpJOnLOHCRJDcNBktQwHCRJDcOhj5NO7rXQPBYXwiUday5I9/Grg1O90OxCuKRjzZnDGHrORPyXvqSxOHMYw5TPRCTJcJhm3QxDkpaa4TDNnGFIGolrDpKkhuFwIum5EJ6Ek0893UVz6QTmZaUTSc/LVDB7qepYXtJauXoNe2fuWrTfo1et5r/vurPXa0oajuGgo3MYi+ZjrJ8YStKRMRx0dKZ80dwvCEpHxjUHTZexviB4GOsxrrPoRODMQdOl70zkVU8/tt8BOcz1GOl4ZzjooWnMy1k911mWnXIa9//vfYv2c71D02jqwiHJBuDvgWXAe6rq6pFLkh7sMILJ9Q49VE3VmkOSZcA7gD8EzgNemuS8cauSBjbSOkvfW733/c6L6zHHl2mbOVwA7Kmq7wMkuQa4FLhl1KqkIR3jdZa+l7Og/8eLe6/HHOMaj3W/h8IlvGn5+HWqarAXP1xJXgRsqKo/7fZfDvxeVb16os9mYHO3+7vAbUfxlmcBPzqK84dkbUdmmmuD6a7P2o7MNNcGc9f321W1fKGTpm3mMNc/OR6UXlW1Fdh6TN4s2VVV64/Fax1r1nZkprk2mO76rO3ITHNtcOT1TdWaAzADrJ7YXwXsHakWSTphTVs4fA1Yl+TcJKcCG4EdI9ckSSecqbqsVFUHk7wa+AyzH2V9b1XtHvAtj8nlqYFY25GZ5tpguuuztiMzzbXBEdY3VQvSkqTpMG2XlSRJU8BwkCQ1jvtwSLIhyW1J9iS5co7jSfIP3fFvJ3nyFNX2jCQ/TfLN7vG3S1jbe5PsT/LdeY6POW6L1TbmuK1O8vkktybZneQ1c/QZc+z61DfK+CU5PclNSb7V1faWOfqMMnY9axvt5657/2VJvpHkujmOHf64VdVx+2B2Ufu/gMcApwLfAs47pM9zgE8x+x2LC4Ebp6i2ZwDXjTR2TweeDHx3nuOjjFvP2sYctxXAk7vtRwD/OS0/c4dR3yjj143Hw7vtU4AbgQunYex61jbaz133/q8DPjRXDUcybsf7zOH/b8dRVb8EHrgdx6RLgffXrK8Cj0qyYkpqG01VfRH48QJdxhq3PrWNpqr2VdXXu+2fAbcCKw/pNubY9alvFN14/LzbPaV7HPqJmVHGrmdto0myCngu8J55uhz2uB3v4bASmLxJyQztH4Q+fYbQ932f2k1lP5Xk/CWoq6+xxq2v0cctyVrgScz+K3PSVIzdAvXBSOPXXRr5JrAfuL6qpmbsetQG4/3cvQ14A/CreY4f9rgd7+Gw6O04evYZQp/3/Tqz90B5AvB24N+GLuowjDVufYw+bkkeDnwUeG1V3Xvo4TlOWdKxW6S+0cavqu6vqicye3eEC5I8/pAuo41dj9pGGbckzwP2V9XNC3Wbo23BcTvew6HP7TjGumXHou9bVfc+MJWtqk8CpyQ5awlq62Nqb3Uy9rglOYXZv3g/WFUfm6PLqGO3WH1jj1/3vj8BbgA2HHJo9J+7+WobcdwuAp6f5AfMXp5+VpIPHNLnsMfteA+HPrfj2AG8olvNvxD4aVXtm4bakvxWMnv/4yQXMPv/654lqK2PscZtUWOOW/e+/wTcWlVvnafbaGPXp76xxi/J8iSP6rYfBjwb+N4h3UYZuz61jTVuVXVVVa2qqrXM/j3yuap62SHdDnvcpur2GcdazXM7jiR/3h1/F/BJZlfy9wC/AP5kimp7EfCqJAeB/wE2VvfRg6El+TCzn744K8kM8GZmF+FGHbeetY02bsz+K+7lwHe669MAbwLWTNQ32tj1rG+s8VsBbMvsL/06CdheVddNw5/XnrWN+XPXONpx8/YZkqTG8X5ZSZJ0BAwHSVLDcJAkNQwHSVLDcJAkNQwHSVLDcJAkNf4PAZSH8x61OVUAAAAASUVORK5CYII=\n",
"text/plain": "<Figure size 432x288 with 1 Axes>"
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
]
},
{
"metadata": {},
"id": "36b741ff",
"cell_type": "markdown",
"source": "As you can see, negative values are clipped to be 0. And you can use this, just like any other transform, directly in your model."
},
{
"metadata": {},
"id": "dba592c9",
"cell_type": "markdown",
"source": "But there are other things you can do as well, like `stack()` RVs, and then index into them with a binary RV."
},
{
"metadata": {
"trusted": false
},
"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": 69,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": "Sampled value = -0.84\nSampled value = 0.55\nSampled value = 0.70\nSampled value = -0.66\nSampled value = -0.34\n"
}
]
},
{
"metadata": {},
"id": "94b12cb2",
"cell_type": "markdown",
"source": "As you can see, depending on whether `index` is `0` or `1` we either sample from the negative or positive uniform. This also supports fancy indexing, so you can manually create complicated mixture distribution using a `Categorical` like this:"
},
{
"metadata": {
"trusted": false
},
"id": "77432dbe",
"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 z = pm4.Uniform(\"z\", lower=1, upper=2)\n xyz = at.stack([x, y, z])\n index = pm4.Categorical(\"index\", [.3, .3], shape=3)\n\nfor _ in range(5):\n print(\"Sampled value = {}\".format(xyz[index].eval()))",
"execution_count": 55,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": "Sampled value = [1.35927217 0.20451265 1.35927217]\nSampled value = [ 1.01280942 -0.1336057 0.45787874]\nSampled value = [1.15299912 0.59138617 1.15299912]\nSampled value = [1.54220526 1.54220526 1.54220526]\nSampled value = [-0.83659795 1.43210203 -0.83659795]\n"
}
]
},
{
"metadata": {},
"id": "75cc3855",
"cell_type": "markdown",
"source": "## Support for for-loops in `scan`\n\nSometimes you would like to create RVs in for-loops. This is very common in timeseries models. Wouldn't it be cool if we could write a Gaussian random walk like this:"
},
{
"metadata": {
"trusted": false
},
"id": "8750fb7b",
"cell_type": "code",
"source": "with pm4.Model():\n y = [0]\n for t in range(5):\n y_t = pm4.Normal(f\"y_{t}\", sigma=.1)\n y.append(y[-1] + y_t)\n \n random_walk = at.stack(y)\n\nrandom_walk.eval()",
"execution_count": 91,
"outputs": [
{
"data": {
"text/plain": "array([ 0. , 0.01, 0.09, -0. , 0.12, 0.16])"
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
]
},
{
"metadata": {},
"id": "3948624d",
"cell_type": "markdown",
"source": "While that is possible, it is incredibly inefficient because creation and keeping track of individual RVs is very costly. That's why in general we advise users to do everything vectorized. \n\nThis is also an area in which Stan shines where this is efficiently possible. With `aeppl`, however, we now have support for `aesara`-style for loops which come in the form of the `scan` operator. With this, we can create RVs *inside* `aesara.scan()`. So the above for-loop can be equivalently rewritten as:"
},
{
"metadata": {
"trusted": true
},
"id": "9b41024a",
"cell_type": "code",
"source": "with pm4.Model():\n random_walk, _ = aesara.scan(\n fn=lambda Y_tm1: Y_tm1 + pm4.Normal.dist(sigma=.1),\n outputs_info=[{\"initial\": at.as_tensor([0.0])}],\n n_steps=100,\n name=\"Y\"\n )\n\nplt.plot(random_walk.eval())",
"execution_count": 137,
"outputs": [
{
"output_type": "error",
"ename": "FileNotFoundError",
"evalue": "[Errno 2] No such file or directory: '/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/tensor/c_code/dimshuffle.c'",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m/var/folders/mn/0x4pxw0n61lf479ndp07r0gr0000gn/T/ipykernel_42277/786484366.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mpm4\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m random_walk, _ = aesara.scan(\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mY_tm1\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mY_tm1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mpm4\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mNormal\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m.1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0moutputs_info\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m\"initial\"\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_tensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0.0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mn_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/scan/basic.py\u001b[0m in \u001b[0;36mscan\u001b[0;34m(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict, return_list)\u001b[0m\n",
"\u001b[0;32m~/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/tensor/shape.py\u001b[0m in \u001b[0;36mshape_padleft\u001b[0;34m(t, n_ones)\u001b[0m\n",
"\u001b[0;32m~/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/tensor/var.py\u001b[0m in \u001b[0;36mdimshuffle\u001b[0;34m(self, *pattern)\u001b[0m\n",
"\u001b[0;32m~/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/tensor/elemwise.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, input_broadcastable, new_order, inplace)\u001b[0m\n",
"\u001b[0;32m~/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/graph/op.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, func_files, func_name)\u001b[0m\n",
"\u001b[0;32m~/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/graph/op.py\u001b[0m in \u001b[0;36mload_c_code\u001b[0;34m(self, func_files)\u001b[0m\n",
"\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '/Users/twiecki/miniforge3/envs/pymc4b1/lib/python3.8/site-packages/aesara/tensor/c_code/dimshuffle.c'"
]
}
]
},
{
"metadata": {},
"id": "63b6c9d2",
"cell_type": "markdown",
"source": "The syntax is famously cryptic. In brief, we define a `lambda` function we pass to `fn` that is executed inside the for-loop retrieving the previous value. The key here is that when we define the RV inside the `scan`, that this is now efficient in PyMC 4.0.\n\nThis feature makes writing various other timeseries distributions much simpler and more flexible."
},
{
"metadata": {},
"id": "374a426b",
"cell_type": "markdown",
"source": "## Better (and Dynamic) Shape Support\n\nAnother big improvement in `PyMC` 4.0 is in how shapes are handled internally. Before this was also a bunch of complicated Python code because we did not make proper use of `aesara`. Now, everything is completely offloaded to `aesara`. As a sid-effect, this better shape support also allows dynamic RV shapes, where the shape depends on another RV:"
},
{
"metadata": {
"trusted": false
},
"id": "0daf28f1",
"cell_type": "code",
"source": "with pm4.Model() as m:\n x = pm4.Poisson('x', 2)\n z = pm4.Normal('z', shape=x)\n \nfor _ in range(5):\n print(\"Value of z = {}\".format(z.eval()))",
"execution_count": 66,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": "Value of z = [-0.36 -0.87]\nValue of z = [-0.43]\nValue of z = [-1.25 -0.07 0.6 ]\nValue of z = [-0.4 1.94]\nValue of z = [1.58]\n"
}
]
},
{
"metadata": {},
"id": "5b1cec0b",
"cell_type": "markdown",
"source": "As you can see, the shape of `z` changes with each draw according to the integer sampled by `x`.\n\nNote, however, that this does not yet work for posterior inference (i.e. sampling). The reason is that the trace backend (`arviz.InferenceData`) as well as samplers in this case also must support changing dimensionality (like reversible-jump MCMC). There are plans to add this."
},
{
"metadata": {},
"id": "33c1900e",
"cell_type": "markdown",
"source": "## Conditioning & Causal Inference"
},
{
"metadata": {
"trusted": false
},
"id": "ee948605",
"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": 67,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": "10 samples from P(x | mu=0, sigma=0.1) = [ 0.03 -0.15 0.07 -0.02 0.12 -0.01 0.07 0.08 -0.05 0.01]\n"
}
]
},
{
"metadata": {},
"id": "b80cb4cf",
"cell_type": "markdown",
"source": "## Better NUTS initialization\n\nWe have also fixed an issue with the default NUTS warm-up which sometimes lead to the sampler getting stuck for a while. While fixing this issue, [Adrian Seyboldt](https://twitter.com/aseyboldt) also came up with a new initialization method that uses the gradients to estimate a better mass-matrix. You can use this (still experimental) feature by calling `pm.sample(init=\"jitter+adapt_diag_grad\")`.\n\nLet's try this on the hierarchical regression model from above:"
},
{
"metadata": {
"trusted": false
},
"id": "f7c53fe9",
"cell_type": "code",
"source": "model_pymc4_grad, idata_pymc4_grad = build_model(pm4, init=\"jitter+adapt_diag_grad\")",
"execution_count": 95,
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": "Auto-assigning NUTS sampler...\nInitializing NUTS using jitter+adapt_diag_grad...\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]\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"
},
{
"data": {
"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, 0 divergences]\n </div>\n ",
"text/plain": "<IPython.core.display.HTML object>"
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"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 15 seconds.\nThe number of effective samples is smaller than 25% for some parameters.\n"
}
]
},
{
"metadata": {},
"id": "9789798d",
"cell_type": "markdown",
"source": "The first thing to observe as that we did not get any divergences this time. Comparing the effective sample size of the default and grad-based initialization, we can also see that it leads to much better sampling for certain parameters:"
},
{
"metadata": {
"trusted": false
},
"id": "f65f04c3",
"cell_type": "code",
"source": "pd.DataFrame({\"Default init\": az.summary(idata_pymc4, var_names=[\"~a\", \"~b\"])[\"ess_bulk\"],\n \"Grad-based init\": az.summary(idata_pymc4_grad, var_names=[\"~a\", \"~b\"])[\"ess_bulk\"]}).plot.barh()\nplt.xlabel(\"effective sample size (higher is better)\");",
"execution_count": 126,
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEGCAYAAACZ0MnKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgjElEQVR4nO3de5hV1X3/8ffHAUHkoiJYQeuQxBgRxhEHJUENKsU7jRGDVSo2sQSNWjFiaazpmDxJUUnTRiNeqIHUWxoJKuqvsY9i1QSV2zAMWk3UMSIpIFECCAbl+/tj7xkPw1z2DGc4Z4bP63nOc/bZZ+21v+vsmfnOXuvstRURmJmZZbFXoQMwM7OOw0nDzMwyc9IwM7PMnDTMzCwzJw0zM8usS6EDaG8HHnhglJaWFjoMM7MOZcmSJe9GRL+G6zt90igtLWXx4sWFDsPMrEOR9FZj6909ZWZmmTlpmJlZZk4aZmaWmZOGmZll5qRhZmaZOWmYmVlmThpmZpaZk4aZmWXmpGFmZpl1+ivCWb0MKvsUOgozs92rckO7VOszDTMzy8xJw8zMMnPSMDOzzJw0zMwsMycNMzPLzEnDzMwyc9IwM7PMnDTMzCwzJw0zM8usKJKGpAmSXpJUJelOSSWSNkn6gaSlkp6S1C8te5WklyVVS3qw0LGbme1JCj6NiKQjgfHAyIjYJul24CJgX2BpRHxT0reBfwKuAKYBgyLiQ0n7NVHnJGASQEnvfpRu/cluaInZ7lc7/axCh2B7mGI40zgVOBZYJKkqff0pYDvws7TMvcAJ6XI1cJ+kCcBHjVUYEXdFREVEVJT08LxTZmb5UgxJQ8CciChPH0dERGUj5SJ9Pgv4MUmiWSKp4GdLZmZ7imJIGk8B4yT1B5B0gKTDSGIbl5a5EHhe0l7AoRGxALgO2A/ouftDNjPbMxX8v/SIeFnSPwJPpklhG/ANYDNwlKQlwAaScY8S4F5JfUjOUH4YEe8XJnIzsz1PwZMGQET8jE/GLwCQRETcANzQoPgJmJlZQRRD95SZmXUQRZs0IsJjFWZmRaZok4aZmRUfJw0zM8vMScPMzDIrim9PtaehA/uw2FMtmJnlhc80zMwsMycNMzPLzEnDzMwyc9IwM7PMnDTMzCwzJw0zM8vMScPMzDJz0jAzs8ycNMzMLDMnDTMzy8xJw8zMMnPSMDOzzJw0zMwsMycNMzPLzEnDzMwyc9IwM7PMnDTMzCwzJw0zM8vMScPMzDJz0jAzs8ycNMzMLLMuhQ6g3a1eBpV9Ch1F+6ncUOgIzGwP4jMNMzPLzEnDzMwyc9IwM7PMnDTMzCwzJw0zM8vMScPMzDJz0jAzs8zalDQkzZI0ON/BZNhvraQDd/d+zcws0aaL+yLi0nwHYmZmxa/FMw1J+0p6XNJySTWSxkt6RlJF+v7XJL2Wrrtb0m3p+tmSZkpaIOkNSV+UdI+kVyTNzql/pqTFklZKujFDzFMlvZQ+PtPWhpuZWetlOdM4HVgdEWcBSOoDXJYuDwBuAIYBG4GngeU52+4PnAKMBeYDI4FLgUWSyiOiCrg+Iv4gqQR4SlJZRFQ3E88fI+I4SRcD/wqc3bCApEnAJICS3v0o3fqTDM3soKY9Xr9YO/2sAgZiZnuCLGMaK4DRkm6SdGJE5E52dBzwPxHxh4jYBvy8wbbzIyLSOtZExIqI2A6sBErTMl+RtBRYBhwFtDRW8kDO8+cbKxARd0VERURUlPToxPNOmZntZi2eaUTEa5KOBc4E/lnSkzlvq4XNP0yft+cs173uImkQcC0wPCLeS7uturcUUhPLZmbWzrKMaQwAPoiIe4EZJF1RdV4Cvihpf0ldgPNauf/ewGZgg6SDgDMybDM+53lhK/dnZma7IMuYxlDgFknbgW0k4xkzACLiHUnfB14EVgMvA5nn6o6I5ZKWkXRXvQH8KsNm3SS9SJLw/irrvszMbNcpGXLYhQqknhGxKT3TmAfcExHz8hJdHnQ7+PA4eOK/FjqM3cID4WaWL5KWRERFw/X5uCK8UlIVUAO8CTychzrNzKwI7fKd+yLi2nwEkkvSPGBQg9V/HxG/zPe+zMwsu6K83WtEnFvoGMzMbGeesNDMzDIryjONfBo6sA+LPUBsZpYXPtMwM7PMnDTMzCwzJw0zM8vMScPMzDJz0jAzs8ycNMzMLDMnDTMzy8xJw8zMMnPSMDOzzJw0zMwsMycNMzPLzEnDzMwyc9IwM7PMnDTMzCwzJw0zM8vMScPMzDJz0jAzs8ycNMzMLDMnDTMzy8xJw8zMMnPSMDOzzLoUOoB2t3oZVPYpdBRmnVPlhkJHYLuZzzTMzCwzJw0zM8vMScPMzDJz0jAzs8ycNMzMLDMnDTMzy8xJw8zMMnPSMDOzzDpM0pBUKqmm0HGYme3JOkzSMDOzwmu3aUQklQL/BTwPjACWAz8BbgT6AxcBZwKbImJGuk0NcHZE1DYVr6Q5wDHAa8DFEfFBI/ueBEwCKOndj9KtP8lfw8yaUDv9rEKHYNbu2vtM4zPAvwFlwOeAC4ETgGuBb7WhviOAuyKiDPgjcHljhSLiroioiIiKkh6ed8rMLF/aO2m8GRErImI7sBJ4KiICWAGUtqG+tyPiV+nyvSQJyMzMdpP2Thof5ixvz3m9naRr7KMGMXRvob5o4bWZmbWjQg+E1wLDACQNAwa1UP7PJX0+Xf4rkvESMzPbTQqdNOYCB0iqAi4jGdxuzivAREnVwAHAzPYNz8zMcrXbt6fSb0ANyXl9SRPvjWlFfYPzFZ+ZmbVeoc80zMysAym6271K6gs81chbp0bE+t0dj5mZfaLokkaaGMoLHYeZme3M3VNmZpZZ0Z1p5NvQgX1Y7OkdzMzywmcaZmaWmZOGmZll5qRhZmaZOWmYmVlmThpmZpaZk4aZmWXmpGFmZpk5aZiZWWZOGmZmlpmThpmZZeakYWZmmTlpmJlZZk4aZmaWmZOGmZll5qRhZmaZOWmYmVlmThpmZpaZk4aZmWXmpGFmZpk5aZiZWWZOGmZmllmXQgfQ7lYvg8o+hY7CclVuKHQEZtZGPtMwM7PMnDTMzCwzJw0zM8vMScPMzDJz0jAzs8ycNMzMLDMnDTMzy6xNSUPSLEmD8x2MmZkVtzZd3BcRl+Y7EDMzK34tnmlI2lfS45KWS6qRNF7SM5Iq0ve/Jum1dN3dkm5L18+WNFPSAklvSPqipHskvSJpdk79MyUtlrRS0o0txPJtSYvSOO6SpF1sv5mZtUKWM43TgdURcRaApD7AZenyAOAGYBiwEXgaWJ6z7f7AKcBYYD4wErgUWCSpPCKqgOsj4g+SSoCnJJVFRHUTsdwWEd9J9/0fwNlpvTuQNAmYBFDSux+lW3+SoZm220x7vNm3a6eftZsCMbPWyjKmsQIYLekmSSdGRO7EQccB/xMRf4iIbcDPG2w7PyIirWNNRKyIiO3ASqA0LfMVSUuBZcBRQHNjJSdLelHSCpJkdFRjhSLiroioiIiKkh6ed8rMLF9aPNOIiNckHQucCfyzpCdz3m6pe+jD9Hl7znLd6y6SBgHXAsMj4r2026p7YxVJ6g7cDlRExNuSKpsqa2Zm7SPLmMYA4IOIuBeYQdIVVecl4IuS9pfUBTivlfvvDWwGNkg6CDijmbJ1CeJdST2Bca3cl5mZ7aIsYxpDgVskbQe2kYxnzACIiHckfR94EVgNvAxknvc6IpZLWkbSXfUG8Ktmyr4v6W6Srq5aYFHW/ZiZWX4oGXLYhQqknhGxKT3TmAfcExHz8hJdHnQ7+PA4eOK/FjoMawUPhJsVnqQlEVHRcH0+rgivlFQF1ABvAg/noU4zMytCu3znvoi4Nh+B5JI0DxjUYPXfR8Qv870vMzPLrihv9xoR5xY6BjMz25knLDQzs8yK8kwjn4YO7MNiD6yameWFzzTMzCwzJw0zM8vMScPMzDJz0jAzs8ycNMzMLDMnDTMzy8xJw8zMMnPSMDOzzJw0zMwsMycNMzPLzEnDzMwyc9IwM7PMnDTMzCwzJw0zM8vMScPMzDJz0jAzs8w6/U2YzKz9bNu2jVWrVrF169ZCh2Jt1L17dw455BC6du2aqbyThpm12apVq+jVqxelpaVIKnQ41koRwfr161m1ahWDBg3KtI27p8yszbZu3Urfvn2dMDooSfTt27dVZ4pOGma2S5wwOrbWHj8nDTMzy6zzj2msXgaVfQodhRVK5YZCR7BHKZ32eF7rq51+VotlSkpKGDp0KNu2baNLly5MnDiRq6++mr32av5/4qlTp/LEE09w5plncsstt7Q6tp49e7Jp0yZqa2v59a9/zYUXXrhTmdWrV3PVVVfx0EMPNVvXmWeeyf333w/A/fffz+WXX97qeHaXzp80zKxT22effaiqqgJg7dq1XHjhhWzYsIEbb7yx2e3uvPNO1q1bR7du3XZp/7W1tdx///2NJo0BAwa0mDAAnnjiifq6br/99qJOGu6eMrNOo3///tx1113cdtttRAQff/wxU6dOZfjw4ZSVlXHnnXcCMHbsWDZv3szxxx/Pz372M+bPn8/xxx/PMcccw+jRo1mzZg0AlZWVzJgxo77+IUOGUFtbu8M+p02bxnPPPUd5eTk//OEPd3ivtraWIUOGADB79my+/OUvc/rpp3P44Ydz3XXX1ZcrLS3l3XffZdq0abz++uuUl5czderU9viIdpnPNMysU/nUpz7F9u3bWbt2LY888gh9+vRh0aJFfPjhh4wcOZIxY8bw6KOP0rNnz/ozlPfee48XXngBScyaNYubb76ZH/zgB5n2N336dGbMmMFjjz3WYtmqqiqWLVtGt27dOOKII7jyyis59NBDd6irpqamPq5i5KRhZp1ORADw5JNPUl1dXd9FtGHDBn7zm9/sdE3CqlWrGD9+PL///e/505/+lPmahdY69dRT6dMnGWMdPHgwb7311g5JoyNw0jCzTuWNN96gpKSE/v37ExHceuutnHbaac1uc+WVV3LNNdcwduxYnnnmGSorKwHo0qUL27dvry+3q1e+546flJSU8NFHH+1SfYXgMQ0z6zTWrVvH5MmTueKKK5DEaaedxsyZM9m2bRsAr732Gps3b95puw0bNjBw4EAA5syZU7++tLSUpUuXArB06VLefPPNnbbt1asXGzduzEv8+ayrvfhMw8zyJstXZPNty5YtlJeX13/l9q//+q+55pprALj00kupra1l2LBhRAT9+vXj4Ycf3qmOyspKzj//fAYOHMiIESPqk8N5553HT3/6U8rLyxk+fDif/exnd9q2rKyMLl26cPTRR3PJJZcwZcqUNrelb9++jBw5kiFDhnDGGWe06avA7U11fX+dVcWAklg8qWehw7BC8XUa7eqVV17hyCOPLHQYtosaO46SlkRERcOy7p4yM7PM2i1pSCqV9L+SZkmqkXSfpNGSfiXpN5KOk1Qp6dqcbWoklTZT58OSlkhaKWlSe8VuZmaNa+8xjc8A5wOTgEXAhcAJwFjgW0BVK+v7akT8QdI+wCJJcyNifcNCaUKZBFDSux+lW3/S9hbYblOI/nAza5327p56MyJWRMR2YCXwVCSDKCuA0jbUd5Wk5cALwKHA4Y0Vioi7IqIiIipKenjeKTOzfGnvM40Pc5a357zenu77I3ZMXN2bqkjSKGA08PmI+EDSM82VNzOz/Cv0QHgtMAxA0jCgucsw+wDvpQnjc8CI9g/PzMxyFfo6jbnAxZKqSMY8Xmum7H8BkyVVA6+SdFGZWTHJ920IMnxles2aNUyZMoUXXniB/fffn7333pvrrruOc889t+27raykZ8+eXHvttTusr62t5eyzz6ampqbNdbfGqFGjmDFjBhUVO37z9dJLL+Waa65h8ODBTW57xx130KNHDy6++GJmz57NmDFjGDBgwC7H1G5JIyJqgSE5ry9p4r0xGev7EDgjbwGaWYcXEXzpS19i4sSJ9fejeOutt3j00Ud3KvvRRx/RpUuh/0/Oj1mzZrVYZvLkyfXLs2fPZsiQIXlJGoXunjIza7Onn36avffee4c/kIcddhhXXnklkPyxPP/88znnnHMYM2YMmzZt4tRTT2XYsGEMHTqURx55pH67733vexxxxBGMHj2aV199tcl9fvTRR0ycOJGysjLGjRvHBx98AMB3vvMdhg8fzpAhQ5g0aVL9pIk/+tGPGDx4MGVlZVxwwQUAbN68ma9+9asMHz6cY445pj6OLVu2cMEFF1BWVsb48ePZsmVLozGMGjWKxYsXA8nNoK6//nqOPvpoRowYsdO07g899BCLFy/moosuory8vMk6syq6pCGpr6SqRh59Cx2bmRWXlStXMmzYsGbLLFy4kDlz5vD000/TvXt35s2bx9KlS1mwYAHf/OY3iQiWLFnCgw8+yLJly/jFL37BokWLmqzv1VdfZdKkSVRXV9O7d29uv/12AK644goWLVpETU0NW7ZsqZ8qffr06Sxbtozq6mruuOMOIElQp5xyCosWLWLBggVMnTqVzZs3M3PmTHr06EF1dTXXX389S5YsafEz2Lx5MyNGjGD58uWcdNJJ3H333Tu8P27cOCoqKrjvvvuoqqpin332abHO5hRd0oiI9RFR3shjp+sxzMxyfeMb3+Doo49m+PDh9ev+4i/+ggMOOABIurO+9a1vUVZWxujRo3nnnXdYs2YNzz33HOeeey49evSgd+/ejB07tsl9HHrooYwcORKACRMm8PzzzwOwYMECjj/+eIYOHcrTTz/NypUrgWRuqosuuoh77723vnvsySefZPr06ZSXlzNq1Ci2bt3K7373O5599lkmTJhQv11ZWVmLbd577705++yzATj22GN3uklUvnWODj4z2yMdddRRzJ07t/71j3/8Y959990dBo733Xff+uX77ruPdevWsWTJErp27UppaWn9dOeSdqr/7bff5pxzzgGSMYLTTz99p3KS2Lp1K5dffjmLFy/m0EMPpbKysr7exx9/nGeffZZHH32U7373u6xcuZKIYO7cuRxxxBE77bOxOJrTtWvX+m12x3TrRXemYWaW1SmnnMLWrVuZOXNm/bq6MYbGbNiwgf79+9O1a1cWLFjAW2+9BcBJJ53EvHnz2LJlCxs3bmT+/PlAclZRVVVFVVVV/bjJ7373OxYuXAjAAw88wAknnFCfIA488EA2bdpUf9On7du38/bbb3PyySdz88038/7777Np0yZOO+00br311vpxj2XLltXHcd999wFQU1NDdXV1Xj6nfE653unPNIYO7MNiT09htnvs5lmFJfHwww8zZcoUbr75Zvr168e+++7LTTfd1Gj5iy66iHPOOYeKigrKy8v53Oc+B8CwYcMYP3485eXlHHbYYZx44olN7vPII49kzpw5fP3rX+fwww/nsssuo0ePHvzt3/4tQ4cOpbS0tL577OOPP2bChAls2LCBiGDKlCnst99+3HDDDVx99dWUlZUREZSWlvLYY49x2WWX8Td/8zeUlZVRXl7Occcdl5fP6ZJLLmHy5Mnss88+LFy4cJfGNTr/1OgVFVH3LQMzyy9Pjd45eGp0MzNrF04aZmaWmZOGme2Szt7F3dm19vg5aZhZm3Xv3p3169c7cXRQEcH69evp3j37hOGd/ttTZtZ+DjnkEFatWsW6desKHYq1Uffu3TnkkEMyl3fSMLM269q1K4MGNXdHA+ts3D1lZmaZOWmYmVlmThpmZpZZp78iXNJGkjv9dQYHAu8WOog8cVuKk9tSnArRlsMiol/DlXvCQPirjV0K3xFJWuy2FB+3pTi5Le3D3VNmZpaZk4aZmWW2JySNuwodQB65LcXJbSlObks76PQD4WZmlj97wpmGmZnliZOGmZll1mmThqTTJb0q6beSphU6niwk1UpaIalK0uJ03QGS/lvSb9Ln/XPK/0PavlclnVa4yEHSPZLWSqrJWdfq2CUdm34Gv5X0I0kqkrZUSnonPTZVks7sIG05VNICSa9IWinp79L1He7YNNOWDndsJHWX9JKk5WlbbkzXF/9xiYhO9wBKgNeBTwF7A8uBwYWOK0PctcCBDdbdDExLl6cBN6XLg9N2dQMGpe0tKWDsJwHDgJpdiR14Cfg8IOD/AWcUSVsqgWsbKVvsbTkYGJYu9wJeS2PucMemmbZ0uGOT7rdnutwVeBEY0RGOS2c90zgO+G1EvBERfwIeBP6ywDG11V8Cc9LlOcCXctY/GBEfRsSbwG9J2l0QEfEs8IcGq1sVu6SDgd4RsTCS34af5myz2zTRlqYUe1t+HxFL0+WNwCvAQDrgsWmmLU0p5rZERGxKX3ZNH0EHOC6dNWkMBN7Oeb2K5n+4ikUAT0paImlSuu6giPg9JL80QP90fUdoY2tjH5guN1xfLK6QVJ12X9V1G3SYtkgqBY4h+a+2Qx+bBm2BDnhsJJVIqgLWAv8dER3iuHTWpNFYn15H+G7xyIgYBpwBfEPSSc2U7ahthKZjL+Y2zQQ+DZQDvwd+kK7vEG2R1BOYC1wdEX9srmgj64qqPY20pUMem4j4OCLKgUNIzhqGNFO8aNrSWZPGKuDQnNeHAKsLFEtmEbE6fV4LzCPpblqTnoKSPq9Ni3eENrY29lXpcsP1BRcRa9Jf8u3A3XzSFVj0bZHUleSP7H0R8Yt0dYc8No21pSMfG4CIeB94BjidDnBcOmvSWAQcLmmQpL2BC4BHCxxTsyTtK6lX3TIwBqghiXtiWmwi8Ei6/ChwgaRukgYBh5MMiBWTVsWeno5vlDQi/QbIxTnbFFTdL3LqXJJjA0XelnTf/w68EhH/kvNWhzs2TbWlIx4bSf0k7Zcu7wOMBv6XjnBcduc3BnbnAziT5NsVrwPXFzqeDPF+iuTbEcuBlXUxA32Bp4DfpM8H5Gxzfdq+VynAN3MaxP8ASdfANpL/fr7WltiBCpJf+teB20hnLSiCtvwHsAKoJvkFPriDtOUEku6KaqAqfZzZEY9NM23pcMcGKAOWpTHXAN9O1xf9cfE0ImZmllln7Z4yM7N24KRhZmaZOWmYmVlmThpmZpaZk4aZmWXmpGE7kXR+OpPogvT1A+kUDVNaWc9+ki7PeT1A0kP5jre9KZl9+MBdrGOypIvzEMsxkmaly5WSrm2i3K8z1LXL7WrtPtu6b0mjJH0h5/WXJA1uQ4xn180oa23jpGGN+RpweUScLOnPgC9ERFlE/LCV9ewH1CeNiFgdEePyGGeHERF3RMRP81DVt4BbM+zvCy2V2RWSuuzmfY4Ccuv/EsnMr5mlMT8OjJXUI2+R7WGcNPZgkiYomdO/StKd6QRq3ya5iOoOSbcATwL90zInSvq0pP9KJ1V8TtLn0roOkjRPyf0Blqf/FU4HPp1ue4ukUqX3qJD0oqSjcmJ5Rsl9AfZNJ51bJGmZpJ1mJ5Z0sKRn03prJJ2Yrp8pabFy7k+Qrq+V9H1JC9P3h0n6paTXJU1Oy4xK65wn6WVJd0ja6fejsc+skTLT0zqqJc1I11VKujY926rKeXws6bD0CuG5absXSRrZSL29gLKIWJ6zenD62b0h6aqcspvS570k3Z5+Jo9JekJSbuK+UtJSJfdjqDuWjR4DSZdI+rmk+enPRcP46vbZ6PFpxNT0s3xJ0mfSbXf6HJRMTjgZmJLW+UVgLHBL+vrTzfxczpb0L0rOmm+K5MK0Z4Czm4jJWrI7r4L0o3gewJHAfKBr+vp24OJ0+RmgIl0uZcf7SjwFHJ4uHw88nS7/jGQCOUjuZ9KnkW3rXwNTgBvT5YOB19Ll7wMT0uX9SK7q37dB7N/kkyvmS4Be6fIBOeueIfkDC8l9Si5Ll39IchVuL6AfsDZdPwrYSnJlfgnw38C4nO0PbO4zy4ntAJIrdusunN0vfa6kwT0fgG8A/5ku3w+ckC7/OclUGQ2P2cnA3JzXlcCvSe6xcCCwPie2TenzOOAJkn8Q/wx4r0G7rkyXLwdmNXcMgEtIrpA/oGFsDfbZ6PFpULY2p8zFwGPNfQ4NPz9gdl07Wvi5nA08Rs69ZoCLgFsL/TvYUR87nWLaHuNU4FhgkZIbfe3DJ5OjNUrJ7KJfAH6uT24O1i19PoXkl5+I+BjYoJy7jjXiP0n+MP8T8BXg5+n6MSTdB3V99d1J/3jkbLsIuEfJ5HUPR0RVuv4rSqaU70KSiAaTJAj4ZO6xFSQ3v9lIMmfPVqVzAJHM5fNG2tYHSM64csdgsnxmfyRJPrMkPU7yB2sn6ZnEpUDdf+GjSc4a6or0ltQrjbPOwcC6BlU9HhEfAh9KWgscxI5TZZ8A/DySyfz+L/2PO1fdBIZLgC+ny00dA0im8G7pXiNNHZ+GHsh5ruv6bPRzaG5nLfxcQtL+j3NerwUGtNAGa4KTxp5LwJyI+IdWbLMX8H4k0znvkoh4R9J6SWXAeODrOXGdFxGvNrPts0qmjT8L+A8l3WjPAdcCwyPiPUmzSf7Y1fkwfd6es1z3uu73oOGcOg1ft/iZRcRHko4jSTAXAFeQJNRPKkkm2Pt3YGx8ciOevYDPR8SWpuoGtjRoEw3a8jE7/063dOvPuu1zt230GEg6HtjcQn2NHp9ofDwnGllu9HNQ83cwbennsmHM3Uk+S2sDj2nsuZ4CxknqD/X3Jj6suQ0iuXfBm5LOT7eRpKNz6rssXV8iqTewkaQbqCkPAtcBfSJiRbrulyT97ErrOqbhRmmcayPibpI/vsOA3iR/HDZIOojkniStdZySmZH3Iklkzzd4v8XPLP2vt09EPAFcTXKPh9z3u5KcZf19RLyW89aTJAmmrtwO26VeAT7TyjY9D5yXjm0cRNIN15IWj0Fzmjg+jRmf87wwXW7qc2j4s1T/uoWfy8Z8lk9mwrVWctLYQ0XEy8A/ktwpsJqkq+jg5rcCkv7gr0mqm423bqD674CTJa0g6eo4KiLWA79KB0NvaaSuh0j+G//PnHXfJbn1ZbWSQfPvNrLdKKBK0jLgPODfIhkcXpbGdA/wqwxtaWghyeB9DfAmyT1N6mX8zHoBj6Xv/w/J2E2uLwDDgRv1yWD4AOAqoELJ4PnLJAO/O4iI/wX6tNRd08Bcku6qGuBOkjvdbWhhmyzHoDmjaHB8mijXTdKLJD87dZ9TU5/DfODc9PM6keQfjqlKBuo/TdM/l405meRbVNYGnuXWjOTbUyQDrUX9rRol18psjIhZrdimZ0RsktSX5J4rIyPi/9otyCKWnm3dHxGnFjqWjspjGmYdy0zg/FZu81g62L838N09NWGk/pzk213WRj7TMDOzzDymYWZmmTlpmJlZZk4aZmaWmZOGmZll5qRhZmaZ/X/B2LR1+G2SCgAAAABJRU5ErkJggg==\n",
"text/plain": "<Figure size 432x288 with 1 Axes>"
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
]
},
{
"metadata": {},
"id": "b7f4454d",
"cell_type": "markdown",
"source": "## What's not yet working\n\nCurrently we are still missing a few features from PyMC3. These should not be hard to add and will be there for the 4.0 stable release. The list includes:\n\n* Timeseries distributions\n* Mixture Distributions\n\n## Accolades\n\nWhile many people contributed to this effort, we would like to highlight the outstanding contributions of [Brandon Willard](https://brandonwillard.github.io), [Ricardo Vieira](https://github.com/ricardoV94), and [Kaustubh Chaudhari](https://github.com/kc611) who did most of the work."
}
],
"metadata": {
"_draft": {
"nbviewer_url": "https://gist.github.com/badf23b78fd1237d16f4c1909ddbe3e3"
},
"gist": {
"id": "badf23b78fd1237d16f4c1909ddbe3e3",
"data": {
"description": "PyMC3 vs PyMC 4.0.ipynb",
"public": true
}
},
"kernelspec": {
"name": "pymc4b1",
"display_name": "pymc4b1",
"language": "python"
},
"language_info": {
"name": "python",
"version": "3.9.9",
"mimetype": "text/x-python",
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"pygments_lexer": "ipython3",
"nbconvert_exporter": "python",
"file_extension": ".py"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment