# Quadratic Approximation
Modern Bayesian data analysis is based on efficient Markov Chain Monte Carlo (MCMC) techniques.
Learning Bayesian statistic can be hard for a Frequentist! The [Statistical Rethinking](https://xcelab.net/rm/) book by Richard McElreath tries to avoid cognitive overload for its readers by not requireing them to learn Bayesian statistics and MCMC at the same time. The posterior distribution is in most interesting cases analytically intractable and hence MCMC is used to numerically determine it.
As a simpler alternative to MCMC one can you Quadratic Approximation[^1]. A lot of people has ported the R code examples from Statistical Rethinking to Python using frameworks like PYMC, Pyro, NumPyro, and TensorFlow Probability. Apparently there is no Quadratic Approximation solution available for PYMC[^2]. Numpyro has [AutoLaplaceApproximation](https://num.pyro.ai/en/latest/autoguide.html#numpyro.infer.autoguide.AutoLaplaceApproximation)[^3] but it is not clear to me that this is the same as Quadratic Approximation.
quad5 leverages [PYMC](https://www.pymc.io/welcome.html) and works by adding a custom step for the sample method on PYMC models. By doing so it benefits from standard PYMC functionality that automatically adds constant, deterministic, and observed data to the
[inferencedata](https://python.arviz.org/en/stable/getting_started/WorkingWithInferenceData.html) output from the sample method. This allows for easy usage of the [Arviz](https://www.arviz.org/en/latest/) visualization library for posterior distributions and more general the ecosystem around PYMC.
## Word of warning
This package is not production-grade code. It is primarily meant for educational purposes. Secondary for the author to learn more about the internals of the PYMC library.
I am quite sure there will be valid PYMC models where this package not is able to produce a quadratic approximation for the posterior distribution. You are more than welcome to submit either a PR or create a issue for such cases.
## Example
``` python
import arviz as az
import numpy as np
import pymc as pm
import quad5 as quad5
y = np.array([2642, 3503, 4358], dtype=np.float64)
with pm.Model() as m:
logsigma = pm.Uniform("logsigma", 1, 100)
mu = pm.Uniform("mu", 0, 10000)
_ = pm.Normal("y", mu=mu, sigma=pm.math.exp(logsigma), observed=y)
custom_step = quad5.QuadraticApproximation([mu, logsigma], m)
trace = pm.sample(draws=1000, chains=4, tune=100, step=custom_step)
az.plot_posterior(trace)
```
See more examples in this [notebook](https://colab.research.google.com/github/carsten-j/Rethinking/blob/main/chapter4.ipynb) with examples from chapter 4 in Statistical Rethinking.
[^1]: [The Bernstein-Von Mises Theorem](https://en.wikipedia.org/wiki/Bernstein%E2%80%93von_Mises_theorem) states that under some regularity conditions the posterior distribution is asymptotically normal. If the distribution is unimodal with most of the probability mass located symmetric around the peak then quite often you will get a faily good approximation using Quadratic Approximation.
[^2]: This work is partly based on the Python package [pymc3-quap](https://github.com/rasmusbergpalm/pymc3-quap) but pymc3-quap is based on PYMC3 and a lot happend bewteen version 3 and 5 of PYMC. Optimizers works better when provided with a good initial guess and hence a (optional) starting point has been added to function arguments. Please see [Github](https://github.com/pymc-devs/pymc/issues/5443#issuecomment-1030609090) for a discussion about the differences between PYMC version 3 and 5 for computing the Hessian and in particular the for loop `for var in vars:` used when computing the Hessian.
[^3]:The NumPyro documentation refers to "Automatic Guide Generation" and as I understand it this is a kind
of variational inference of the posterior distribution.
Raw data
{
"_id": null,
"home_page": "https://github.com/carsten-j/quad5",
"name": "quad5",
"maintainer": null,
"docs_url": null,
"requires_python": null,
"maintainer_email": null,
"keywords": "Bayesian Statistics, Probabalistic Programming Language, Python",
"author": "Carsten J\u00f8rgensen",
"author_email": "carstenj@gmail.com",
"download_url": "https://files.pythonhosted.org/packages/65/f6/511d3d18d80bedabc118225fee2659964043d215aeadbecaf7594fc2a1ab/quad5-0.1.14.tar.gz",
"platform": null,
"description": "# Quadratic Approximation\n\nModern Bayesian data analysis is based on efficient Markov Chain Monte Carlo (MCMC) techniques. \n\nLearning Bayesian statistic can be hard for a Frequentist! The [Statistical Rethinking](https://xcelab.net/rm/) book by Richard McElreath tries to avoid cognitive overload for its readers by not requireing them to learn Bayesian statistics and MCMC at the same time. The posterior distribution is in most interesting cases analytically intractable and hence MCMC is used to numerically determine it.\n\nAs a simpler alternative to MCMC one can you Quadratic Approximation[^1]. A lot of people has ported the R code examples from Statistical Rethinking to Python using frameworks like PYMC, Pyro, NumPyro, and TensorFlow Probability. Apparently there is no Quadratic Approximation solution available for PYMC[^2]. Numpyro has [AutoLaplaceApproximation](https://num.pyro.ai/en/latest/autoguide.html#numpyro.infer.autoguide.AutoLaplaceApproximation)[^3] but it is not clear to me that this is the same as Quadratic Approximation.\n\nquad5 leverages [PYMC](https://www.pymc.io/welcome.html) and works by adding a custom step for the sample method on PYMC models. By doing so it benefits from standard PYMC functionality that automatically adds constant, deterministic, and observed data to the \n[inferencedata](https://python.arviz.org/en/stable/getting_started/WorkingWithInferenceData.html) output from the sample method. This allows for easy usage of the [Arviz](https://www.arviz.org/en/latest/) visualization library for posterior distributions and more general the ecosystem around PYMC.\n\n## Word of warning\nThis package is not production-grade code. It is primarily meant for educational purposes. Secondary for the author to learn more about the internals of the PYMC library.\n\nI am quite sure there will be valid PYMC models where this package not is able to produce a quadratic approximation for the posterior distribution. You are more than welcome to submit either a PR or create a issue for such cases.\n\n## Example\n``` python\n\n import arviz as az\n import numpy as np\n import pymc as pm\n import quad5 as quad5\n\n y = np.array([2642, 3503, 4358], dtype=np.float64)\n\n with pm.Model() as m:\n logsigma = pm.Uniform(\"logsigma\", 1, 100)\n mu = pm.Uniform(\"mu\", 0, 10000)\n _ = pm.Normal(\"y\", mu=mu, sigma=pm.math.exp(logsigma), observed=y)\n custom_step = quad5.QuadraticApproximation([mu, logsigma], m)\n trace = pm.sample(draws=1000, chains=4, tune=100, step=custom_step)\n\n az.plot_posterior(trace) \n```\n\n\n\nSee more examples in this [notebook](https://colab.research.google.com/github/carsten-j/Rethinking/blob/main/chapter4.ipynb) with examples from chapter 4 in Statistical Rethinking.\n\n[^1]: [The Bernstein-Von Mises Theorem](https://en.wikipedia.org/wiki/Bernstein%E2%80%93von_Mises_theorem) states that under some regularity conditions the posterior distribution is asymptotically normal. If the distribution is unimodal with most of the probability mass located symmetric around the peak then quite often you will get a faily good approximation using Quadratic Approximation.\n\n[^2]: This work is partly based on the Python package [pymc3-quap](https://github.com/rasmusbergpalm/pymc3-quap) but pymc3-quap is based on PYMC3 and a lot happend bewteen version 3 and 5 of PYMC. Optimizers works better when provided with a good initial guess and hence a (optional) starting point has been added to function arguments. Please see [Github](https://github.com/pymc-devs/pymc/issues/5443#issuecomment-1030609090) for a discussion about the differences between PYMC version 3 and 5 for computing the Hessian and in particular the for loop `for var in vars:` used when computing the Hessian.\n\n[^3]:The NumPyro documentation refers to \"Automatic Guide Generation\" and as I understand it this is a kind\nof variational inference of the posterior distribution.\n",
"bugtrack_url": null,
"license": "MIT",
"summary": "Quadratic Approximation custom step sampler for PYMC.",
"version": "0.1.14",
"project_urls": {
"Homepage": "https://github.com/carsten-j/quad5"
},
"split_keywords": [
"bayesian statistics",
" probabalistic programming language",
" python"
],
"urls": [
{
"comment_text": "",
"digests": {
"blake2b_256": "23c034aeacb9719cf7d71b46ca25dc5e605cb0f3b94e95afa264c54aed4cd521",
"md5": "1eed41c4249ebccbdd944d2150c4144d",
"sha256": "2f864873183e522f602db5ff7d6063eecdd03f3a5b69a53c31b3cf651f00bc93"
},
"downloads": -1,
"filename": "quad5-0.1.14-py3-none-any.whl",
"has_sig": false,
"md5_digest": "1eed41c4249ebccbdd944d2150c4144d",
"packagetype": "bdist_wheel",
"python_version": "py3",
"requires_python": null,
"size": 4962,
"upload_time": "2024-05-25T15:21:37",
"upload_time_iso_8601": "2024-05-25T15:21:37.150138Z",
"url": "https://files.pythonhosted.org/packages/23/c0/34aeacb9719cf7d71b46ca25dc5e605cb0f3b94e95afa264c54aed4cd521/quad5-0.1.14-py3-none-any.whl",
"yanked": false,
"yanked_reason": null
},
{
"comment_text": "",
"digests": {
"blake2b_256": "65f6511d3d18d80bedabc118225fee2659964043d215aeadbecaf7594fc2a1ab",
"md5": "b2323dca9b2f3bea6686fc436dec0e45",
"sha256": "50ca37dd6e5e3d936005e59fd8838587970fb47e36e9f68cc7bf6d48284534f7"
},
"downloads": -1,
"filename": "quad5-0.1.14.tar.gz",
"has_sig": false,
"md5_digest": "b2323dca9b2f3bea6686fc436dec0e45",
"packagetype": "sdist",
"python_version": "source",
"requires_python": null,
"size": 5148,
"upload_time": "2024-05-25T15:21:38",
"upload_time_iso_8601": "2024-05-25T15:21:38.138634Z",
"url": "https://files.pythonhosted.org/packages/65/f6/511d3d18d80bedabc118225fee2659964043d215aeadbecaf7594fc2a1ab/quad5-0.1.14.tar.gz",
"yanked": false,
"yanked_reason": null
}
],
"upload_time": "2024-05-25 15:21:38",
"github": true,
"gitlab": false,
"bitbucket": false,
"codeberg": false,
"github_user": "carsten-j",
"github_project": "quad5",
"travis_ci": false,
"coveralls": false,
"github_actions": true,
"requirements": [],
"lcname": "quad5"
}
JSON
