{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "f5199f38", "metadata": { "tags": [ "hide-output" ] }, "outputs": [], "source": [ "import arviz as az\n", "import pymc as pm\n", "from pymc.math import log, sqr\n", "\n", "%load_ext lab_black" ] }, { "cell_type": "markdown", "id": "a3dab390", "metadata": {}, "source": [ "# The Zero Trick and Custom Likelihoods*\n", "\n", "## Zero-trick Jeremy\n", "\n", "This introduces the \"zero trick\", which is a method for specifying custom likelihoods in BUGS. For a more detailed treatment of these methods, see {cite:t}```ntzoufras2009bayesian```, page 276, which is where I got this explanation. These tricks are unnecessary in PyMC since we can just define custom distributions directly, but they do seem to work.\n", "\n", "Adapted from [Unit 6: zerotrickjeremy.odc](https://raw.githubusercontent.com/areding/6420-pymc/main/original_examples/Codes4Unit6/zerotrickjeremy.odc)." ] }, { "cell_type": "markdown", "id": "0fbf22cd-65e7-400e-9f60-0bc3ed9c3799", "metadata": {}, "source": [ "Here's the model we're using:\n", "\n", "$$\n", "\\begin{align*}\n", "X_i \\mid \\theta &\\sim N(\\theta, 80) \\\\\n", "\\theta &\\sim N(110, 120)\n", "\\end{align*}\n", "$$\n" ] }, { "cell_type": "markdown", "id": "2268b05c-3242-4925-bf4c-0e9fcb95a756", "metadata": { "jp-MarkdownHeadingCollapsed": true }, "source": [ "Of course, BUGS can handle this model just fine. But, let's pretend for a moment that there is no built-in normal distribution. The zero trick takes advantages of some properties of the Poisson distribution to recreate an arbitrary likelihood.\n", "\n", "Given a log-likelihood of the form $l_i = \\log f(y; \\theta)$,\n", "\n", "$$\n", "\\begin{align*}\n", "f(y | \\theta) &= \\prod_{i=1}^{n} e^{l_i} && \\text{Log-likelihood to likelihood.}\\\\\n", "&= \\prod_{i=1}^{n} \\frac{e^{-(-l_i)}(-l_i)^0}{0!} && \\text{Matching the form of the Poisson PMF.} \\\\\n", "&= \\prod_{i=1}^{n} f_P(0; -l_i) && \\text{Poisson evaluated at zero with mean $-l_i$.}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "id": "74d36843-3623-4f7b-bde4-5aec9975ed4f", "metadata": { "jp-MarkdownHeadingCollapsed": true }, "source": [ "–{cite:t}```ntzoufras2009bayesian``` page 276.\n", "\n", "But the rate, $\\lambda$, can't be negative. So we need to add a constant, C, to keep that from happening.\n", " \n", "The normal log-likelihood is:\n", "\n", "$$l_i = -0.5 \\log(2\\pi) - \\log(\\sigma^2) - \\frac{(y_i - \\mu_i)^2}{2\\sigma^2}$$\n", "\n", "But the constant terms won't affect the posterior.\n", "\n", "Here's the model in PyMC:" ] }, { "cell_type": "code", "execution_count": 8, "id": "b99d7cd4", "metadata": {}, "outputs": [], "source": [ "y = 98\n", "μ = 110\n", "σ = 80**0.5\n", "τ = 120**0.5\n", "C = 10000\n", "\n", "inits = {\"θ\": 100}" ] }, { "cell_type": "code", "execution_count": 3, "id": "258031fa", "metadata": { "tags": [ "hide-output" ] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [θ]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [24000/24000 00:01<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 1 seconds.\n" ] } ], "source": [ "with pm.Model() as m:\n", " θ = pm.Flat(\"θ\")\n", "\n", " λ1 = pm.Deterministic(\"λ1\", log(σ) + 0.5 * sqr(((y - θ) / σ)) + C)\n", " λ2 = pm.Deterministic(\"λ2\", log(τ) + 0.5 * sqr(((θ - μ) / τ)) + C)\n", "\n", " z1 = pm.Poisson(\"z1\", λ1, observed=0)\n", " z2 = pm.Poisson(\"z2\", λ2, observed=0)\n", "\n", " trace = pm.sample(5000, tune=1000, initvals=inits, target_accept=0.88)" ] }, { "cell_type": "code", "execution_count": 4, "id": "7dba6172", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_2.5%hdi_97.5%
θ102.8526.97289.6116.956
\n", "
" ], "text/plain": [ " mean sd hdi_2.5% hdi_97.5%\n", "θ 102.852 6.972 89.6 116.956" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(trace, hdi_prob=0.95, var_names=\"θ\", kind=\"stats\")" ] }, { "cell_type": "markdown", "id": "a28f71b2", "metadata": {}, "source": [ "## Custom likelihoods in PyMC\n", "\n", "PyMC is a lot more flexible. For one, many more built-in [distributions](https://www.pymc.io/projects/docs/en/stable/api/distributions.html) are available, including mixtures and zero-inflated models, which is where I've seen the zero trick used most often ({cite:t}```karlisfootball```, {cite:t}```ntzoufras2009bayesian``` page 288).\n", "\n", "If you need a custom likelihood, take a look at the [pm.Potential](https://www.pymc.io/projects/docs/en/stable/api/generated/pymc.Potential.html) or [pm.CustomDist](https://www.pymc.io/projects/docs/en/stable/api/distributions/generated/pymc.CustomDist.html) classes." ] }, { "cell_type": "code", "execution_count": 5, "id": "bcdc0192", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: Sun Aug 06 2023\n", "\n", "Python implementation: CPython\n", "Python version : 3.11.4\n", "IPython version : 8.14.0\n", "\n", "pytensor: 2.14.2\n", "\n", "pymc : 5.7.1\n", "arviz: 0.16.1\n", "\n" ] } ], "source": [ "%load_ext watermark\n", "%watermark -n -u -v -iv -p pytensor" ] }, { "cell_type": "code", "execution_count": null, "id": "ac927526", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" } }, "nbformat": 4, "nbformat_minor": 5 }