-
-
Save ltiao/25661d9e938ee84328bebdb70e29472f to your computer and use it in GitHub Desktop.
A notebook detailing multivariate Gauss-Hermite quadrature
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": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import numpy as np\n", | |
"import itertools\n", | |
"\n", | |
"import matplotlib.pyplot as plt\n", | |
"import seaborn as sns" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"ASPECT = 0.5 * (1 + np.sqrt(5))\n", | |
"WIDTH = 8\n", | |
"HEIGHT = WIDTH / ASPECT" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"rc = {\n", | |
" \"figure.figsize\": (WIDTH, HEIGHT),\n", | |
" \"font.serif\": [\"Times New Roman\"],\n", | |
" \"text.usetex\": True,\n", | |
"}\n", | |
"sns.set(context=\"talk\", style=\"ticks\", palette=\"deep\", font=\"serif\", rc=rc)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Numerical integration with `numpy` and `scipy`\n", | |
"\n", | |
"- Mark van der Wilk (20 Sep 2016)\n", | |
"- Revised: Louis C. Tiao (Oct 2021)\n", | |
"\n", | |
"## Gauss-Hermite quadrature\n", | |
"Used for computing integrals of the type:\n", | |
"\n", | |
"$\\int_{-\\infty}^\\infty e^{-x^2} f(x) dx \\approx \\sum_{n=1}^N w_i f(x_i)$\n", | |
"\n", | |
"### First test: Normalising constant\n", | |
"Simply summing the weights should give the normalising constant for a Gaussian with variance 0.5." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Gauss-Hermite: 1.772454\n", | |
"Ground truth : 1.772454\n" | |
] | |
} | |
], | |
"source": [ | |
"x, w = np.polynomial.hermite.hermgauss(10)\n", | |
"print(\"Gauss-Hermite: %f\" % np.sum(w))\n", | |
"print(\"Ground truth : %f\" % (2 * np.pi * 0.5) ** 0.5)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"### Second test: Gaussian moments\n", | |
"We want to calculate expectations w.r.t. a Gaussian distribution with mean $\\mu$ and std $\\sigma$:\n", | |
"\n", | |
"$\\int_{-\\infty}^\\infty \\frac{1}{\\sqrt{2\\pi\\sigma^2}} e^{-\\frac{(y-\\mu)^2}{2\\sigma^2}} f(y) dy$\n", | |
"\n", | |
"We can get this into the form above through a change of variables $x = \\frac{y - \\mu}{\\sqrt{2}\\sigma}$, resulting in:\n", | |
"\n", | |
"$\\int_{-\\infty}^\\infty \\frac{1}{\\sqrt{2\\pi\\sigma^2}} e^{-x^2} f\\left(\\sqrt{2}\\sigma x +\\mu\\right) \\sqrt{2}\\sigma dx = \n", | |
"\\int_{-\\infty}^\\infty \\frac{1}{\\sqrt{\\pi}} e^{-x^2} f\\left(\\sqrt{2}\\sigma x +\\mu\\right) dx$" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 5, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Normalising constant: 1.000000\n", | |
"Mean : 1.230000\n", | |
"Stddev : 2.340000\n" | |
] | |
} | |
], | |
"source": [ | |
"mu, sigma = 1.23, 2.34\n", | |
"const = np.pi**-0.5\n", | |
"y = 2.0**0.5*sigma*x + mu\n", | |
"print(\"Normalising constant: %f\" % np.sum(w * const))\n", | |
"print(\"Mean : %f\" % np.sum(w * const * y))\n", | |
"print(\"Stddev : %f\" % np.sum(w * const * (y - mu)**2.0)**0.5)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"### Third test: More complicated functions\n", | |
"$\\int_{-\\infty}^\\infty \\frac{1}{\\sqrt{2\\pi\\sigma^2}} e^{-\\frac{(y-\\mu)^2}{2\\sigma^2}} \\sin(y) dy$" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 6, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Exp of sin(y) : 0.060982\n" | |
] | |
} | |
], | |
"source": [ | |
"print(\"Exp of sin(y) : %f\" % np.sum(w * const * np.sin(y)))" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Multidimensional Gauss-Hermite quadrature\n", | |
"Integrating against a multi-dimensional Gaussian can be done using Gauss-Hermite quadrature as well, using the correct change of variables. Starting from the general integral:\n", | |
"\n", | |
"$$\n", | |
"\\int \\frac{1}{\\sqrt{\\det{(2\\pi\\boldsymbol{\\Sigma})}}} e^{-\\frac{1}{2}(\\mathbf{y}-\\boldsymbol{\\mu})^{\\top}\\boldsymbol{\\Sigma}^{-1}(\\mathbf{y}-\\boldsymbol{\\mu})} f(\\mathbf{y}) d\\mathbf{y}\n", | |
"$$\n", | |
"\n", | |
"We do the change of variables $\\mathbf{x} = \\frac{1}{\\sqrt{2}}\\mathbf{L}^{-1}(\\mathbf{y}-\\boldsymbol{\\mu})$, where $\\boldsymbol{\\Sigma} = \\mathbf{L}\\mathbf{L}^{\\top}$, giving:\n", | |
"\n", | |
"$$\n", | |
"\\int \\frac{1}{\\sqrt{\\det{(2\\pi\\boldsymbol{\\Sigma})}}} e^{\\mathbf{x}^{\\top} \\mathbf{x}} f\\left(\\sqrt{2}\\mathbf{L}\\mathbf{x} + \\boldsymbol{\\mu}\\right) \\det \\left(\\sqrt{2} \\mathbf{L}\\right) d\\mathbf{x} = \n", | |
"\\int \\pi^{-\\frac{N}{2}} e^{-\\mathbf{x}^{\\top} \\mathbf{x}} f\\left(\\sqrt{2}\\mathbf{L}\\mathbf{x} + \\boldsymbol{\\mu}\\right) d\\mathbf{x}\n", | |
"$$\n", | |
"\n", | |
"If $x$ is $N$ dimensional, this integral above, can now be split into $N$ nested Gauss-Hermite integrals, since the inner product in the exponent can be written as $\\exp\\left(\\sum_{n=1}^N x_n^2\\right)$, which can be written as a product." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 7, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Normalising constant: 1.000000\n", | |
"Mean:\n", | |
"[ 1.00000000e+00 -1.98262334e-17]\n", | |
"Covariance:\n", | |
"[[ 1.3 -0.213]\n", | |
" [-0.213 1.2 ]]\n" | |
] | |
} | |
], | |
"source": [ | |
"mu = np.array([1, 0])\n", | |
"Sigma = np.array([[1.3, -0.213], [-0.213, 1.2]])\n", | |
"N = len(mu)\n", | |
"const = np.pi**(-0.5*N)\n", | |
"xn = np.array(list(itertools.product(*(x,)*N)))\n", | |
"wn = np.prod(np.array(list(itertools.product(*(w,)*N))), 1)\n", | |
"yn = 2.0**0.5*np.dot(np.linalg.cholesky(Sigma), xn.T).T + mu[None, :]\n", | |
"print(\"Normalising constant: %f\" % np.sum(wn * const))\n", | |
"print(\"Mean:\")\n", | |
"print(np.sum((wn * const)[:, None] * yn, 0))\n", | |
"print(\"Covariance:\")\n", | |
"covfunc = lambda x: np.outer(x - mu, x - mu)\n", | |
"print(np.sum((wn * const)[:, None, None] * np.array(list(map(covfunc, yn))), 0))" | |
] | |
} | |
], | |
"metadata": { | |
"anaconda-cloud": {}, | |
"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.8.10" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 1 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment