A notebook detailing multivariate Gauss-Hermite quadrature

numerical-integration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical integration with `numpy` and `scipy`\n",
    "Mark van der Wilk (20 Sep 2016)\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": 2,
   "metadata": {
    "collapsed": false
   },
   "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": 3,
   "metadata": {
    "collapsed": false
   },
   "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": 4,
   "metadata": {
    "collapsed": false
   },
   "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",
    "$\\int \\frac{1}{\\sqrt{\\det{2\\pi\\Sigma}}} e^{-\\frac{1}{2}(y-\\mu)^T\\Sigma^{-1}(y-\\mu)} f(y) dy$\n",
    "\n",
    "We do the change of variables $x = \\frac{1}{\\sqrt{2}}L^{-1}(y-\\mu)$, where $\\Sigma = LL^T$, giving:\n",
    "\n",
    "$\\int \\frac{1}{\\sqrt{\\det{2\\pi\\Sigma}}} e^{x^T x} f\\left(\\sqrt{2}Lx + \\mu\\right)  \\det \\left(\\sqrt{2} L\\right) dx = \n",
    "\\int \\pi^{\\frac{N}{2}} e^{x^T x} f\\left(\\sqrt{2}Lx + \\mu\\right) dx$\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": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normalising constant: 1.000000\n",
      "Mean:\n",
      "[  1.00000000e+00  -2.47886632e-18]\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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

Is the last equation for Multidimensional Gauss-Hermite quadrature? I think two "-" are missing. One is for exp(x^T x), and the other is for pi^(N/2).

Yes, you are correct. I may edit this at some point in the future. :)