34_integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Foundations of Computational Economics #34\n",
    "\n",
    "by Fedor Iskhakov, ANU\n",
    "\n",
    "<img src=\"_static/img/dag3logo.png\" style=\"width:256px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Numerical integration, quadrature\n",
    "\n",
    "<img src=\"_static/img/lecture.png\" style=\"width:64px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"_static/img/youtube.png\" style=\"width:65px;\">\n",
    "\n",
    "[https://youtu.be/cc14S679x2M](https://youtu.be/cc14S679x2M)\n",
    "\n",
    "Description: Gaussian quadrature. Monte Carlo integration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Integration in economics\n",
    "\n",
    "- Expected (discounted) utility  \n",
    "- Expected (discounted) profits  \n",
    "- Bayesian posterior  \n",
    "- Likelihood function with unobservables  \n",
    "- Stochastic elements in (dynamic) economic models  \n",
    "\n",
    "\n",
    "*Most integrals can not be evaluated analytically*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Two main approaches: Monte Carlo and quadrature\n",
    "\n",
    "1. Based on simulations – **Monte Carlo integration**  \n",
    "1. Based on the fixed points and weights – **quadrature integration**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Newton-Cotes formulas\n",
    "\n",
    "Goal: definite integral $ \\int_a^b f(x) dx $\n",
    "\n",
    "Idea: Approximate the function with low order polynomial, then integrate\n",
    "approximation\n",
    "\n",
    "1. First order >> Step function approximation  \n",
    "  - Constant, level at midpoint of $ [a,b] $  \n",
    "1. Second order >> Linear approximation  \n",
    "  - Trapezoid rule  \n",
    "1. Third order >> Quadratic approximation  \n",
    "  - Simpson rule  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"_static/img/composite_simpsons.png\" style=\"width:800px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Composite Newton-Cotes\n",
    "\n",
    "Preform Newton-Cotes on a grid separately on each interval\n",
    "\n",
    "- Equally spaced points  \n",
    "- Newton-Cotes on each sub-interval  \n",
    "\n",
    "\n",
    "*Note that the points are placed exogenously*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"_static/img/newton-cotes.jpg\" style=\"width:800px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Gaussian quadrature\n",
    "\n",
    "General formula\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x) dx = \\sum_{i=1}^n \\omega_i f(x_i)\n",
    "$$\n",
    "\n",
    "- $ x_i \\in [a,b] $ quadrature nodes  \n",
    "- $ \\omega_i $ quadrature weights  \n",
    "\n",
    "\n",
    "*Note that the points are placed endogenously*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Quadrature accuracy\n",
    "\n",
    "Suppose that $ \\{\\phi_k(x)\\}_{k=1,2,\\dots} $ is family of\n",
    "polynomials of degree $ k $ *orthogonal* with respect to the\n",
    "weighting function $ w(x) $\n",
    "\n",
    "- let $ q_k $ denote the leading coefficients so that $ \\phi_k(x)=q_k x^k + \\dots $  \n",
    "- let $ x_i $, $ i=1,\\dots,n $ be $ n $ roots of $ \\phi_n(x) $  \n",
    "- let $ \\omega_i = - \\frac{q_{n+1}/q_n}{\\phi'_n(x_i)\\phi_{n+1}(x_i)}>0 $  \n",
    "\n",
    "\n",
    "Then\n",
    "\n",
    "- $ a<x_1<x_2<\\dots<x_n<b $  \n",
    "- for $ f(x) \\in C^{(2n)}[a,b] $, for some $ \\xi\\in[a,b] $  \n",
    "\n",
    "\n",
    "$$\n",
    "\\int_a^b w(x) f(x) dx = \\sum_{i=1}^n \\omega_i f(x_i) + \\frac{f^{(2n)}(\\xi)}{q_n^2(2n)!}\n",
    "$$\n",
    "\n",
    "- the right hand side is unique on $ n $ nodes  \n",
    "- exact integral for all polynomial $ f(x) $ of degree $ 2n-1 $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Gauss-Chebyshev Quadrature\n",
    "\n",
    "- Domain $ [-1,1] $  \n",
    "- Weighting function $ (1-x^2)^{(-1/2)} $  \n",
    "\n",
    "\n",
    "$$\n",
    "\\int_{-1}^1 \\frac{f(x)dx}{\\sqrt{1-x^2}} = \\frac{\\pi}{n}\\sum_{i=1}^{n} f(x_i) + \\frac{\\pi}{2^{2n-1}}\\frac{f^{(2n)}(\\xi)}{(2n)!}\n",
    "$$\n",
    "\n",
    "- quadrature nodes $ x_i = \\cos(\\frac{2i-1}{2n}\\pi) $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "Want to integrate $ f(x) $ on $ [a,b] $, no weighting function.\n",
    "1. Change of variable $ y=2(x-a)/(b-a)-1 $\n",
    "2. Multiply and divide by weighting function\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x)dx = \\frac{b-a}{2}\\int_{-1}^{1}f\\big(\\frac{(y+1)(b-a)}{2}+a\\big)\\frac{\\sqrt{1-y^2}}{\\sqrt{1-y^2}}dy\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x)dx = \\frac{\\pi(b-a)}{2n}\\sum_{i=1}^{n}f\\big(\\frac{(y_i+1)(b-a)}{2}+a\\big)\\sqrt{1-y_i^2}\n",
    "$$\n",
    "\n",
    "where $ y_i $ are Gauss-Chebyshev nodes over $ [-1,1] $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Gauss-Legendre Quadrature\n",
    "\n",
    "- Domain $ [-1,1] $  \n",
    "- Weighting function $ 1 $  \n",
    "\n",
    "\n",
    "$$\n",
    "\\int_{-1}^1 f(x)dx = \\sum_{i=1}^{n} \\omega_i f(x_i) + \\frac{2^{2n+1}(n!)^4}{(2n+1)!(2n)!} \\frac{f^{(2n)}(\\xi)}{(2n)!}\n",
    "$$\n",
    "\n",
    "- Nodes and weights come from Legendre polynomials, values tabulated  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Gauss-Hermite Quadrature\n",
    "\n",
    "- Domain $ [-\\infty,\\infty] $  \n",
    "- Weighting function $ \\exp(-x^2) $  \n",
    "\n",
    "\n",
    "$$\n",
    "\\int_{-\\infty}^\\infty f(x) \\exp(-x^2)dx = \\sum_{i=1}^{n} \\omega_i f(x_i) + \\frac{n!\\sqrt{\\pi}}{2^n}\\frac{f^{(2n)}(\\xi)}{(2n)!}\n",
    "$$\n",
    "\n",
    "- Nodes and weights come from Hermite polynomials, values tabulated  \n",
    "- Good for computing expectation with Normal distribution  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Gauss-Laguerre Quadrature\n",
    "\n",
    "- Domain $ [0,\\infty] $  \n",
    "- Weighting function $ \\exp(-x) $  \n",
    "\n",
    "\n",
    "$$\n",
    "\\int_{0}^\\infty f(x) \\exp(-x)dx = \\sum_{i=1}^{n} \\omega_i f(x_i) + (n!)^2\\frac{f^{(2n)}(\\xi)}{(2n)!}\n",
    "$$\n",
    "\n",
    "- Nodes and weights come from Laguerre polynomials, values tabulated  \n",
    "- Good for computing expectation exponential discounting  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Multidimensional quadrature\n",
    "\n",
    "Much more complication, simple methods are subject to curse of\n",
    "dimensionality\n",
    "\n",
    "- Generic product rule  \n",
    "\n",
    "\n",
    "$$\n",
    "\\int_{[a,b]^d}f(x)dx=\\sum_{i_1=1}^n \\dots\\sum_{i_d=1}^n \\omega_{i_1}^1\\dots\\omega_{i_1}^d f(x_{i_1}^1,\\dots,x^d_{i_d})\n",
    "$$\n",
    "\n",
    "- Product Gaussian quadrature based on product orthogonal polynomials  \n",
    "- Sparse methods  \n",
    "- Monte Carlo integration!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Monte Carlo integration\n",
    "\n",
    "- Stochastic algorithm for computing integrals  \n",
    "- Main idea: approximate the expectation of a function with an average computed from a sample of random draws  \n",
    "- (convergence in the number of draws is due to the law of large numbers)  \n",
    "- Then convert the integral in expectation to the integral of interest  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Expectation of a function of random variable\n",
    "\n",
    "- let continuous random variable $ \\tilde{X} $ be distributed with pdf $ p(x) $ over domain $ \\Omega $  \n",
    "- we are interested in the expectation of the function $ f(\\tilde{X}) $ which is in turn a random variable itself  \n",
    "\n",
    "\n",
    "$$\n",
    "E\\big(f(\\tilde{X})\\big) = \\int_{\\Omega} f(x) \\,dF(x) = \\int_{\\Omega} f(x) p(x) \\,dx\n",
    "$$\n",
    "\n",
    "- variance of $ f(\\tilde{X}) $ is  \n",
    "\n",
    "\n",
    "$$\n",
    "\\operatorname{Var} \\big(f(\\tilde{X})\\big) = E\\Big(f(\\tilde{X})- E\\big(f(\\tilde{X})\\big)\\Big)^2 = E\\big(f(\\tilde{X})\\big)^2 - \\Big(E\\big(f(\\tilde{X})\\big)\\Big)^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### From expectation to integration\n",
    "\n",
    "- imagine we want to compute the integral denoted $ I_f $ of function $ f(x) $ over some set $ \\Omega $  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f =\\int_{\\Omega} f(x)\\,dx\n",
    "$$\n",
    "\n",
    "- Step 1: represent the integral as an expectation of a function of some random variable $ \\tilde{X} $ defined over domain $ \\Omega $  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f =\\int_{\\Omega} f(x)\\,dx =\n",
    "\\int_{\\Omega} \\frac{f(x)}{p(x)} p(x) \\,dx =E \\left[ \\frac{f(\\tilde{X})}{p(\\tilde{X})} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### From expectation to integration\n",
    "\n",
    "- Step 2: compute the expectation using $ N $ independent draws $ x_i $ of $ \\tilde{X} $ — from the distribution with pdf $ p(x) $  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f =\\int_{\\Omega} f(x)\\,dx = E \\left[ \\frac{f(\\tilde{X})}{p(\\tilde{X})} \\right] \\approx \\frac{1}{N} \\sum_{i=1}^{N} \\frac{f(x_i)}{p(x_i)}\n",
    "\\underset{N \\rightarrow \\infty}{\\rightarrow} E \\left[ \\frac{f(\\tilde{X})}{p(\\tilde{X})} \\right]\n",
    "$$\n",
    "\n",
    "- convergence due to the law of large numbers  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Special simple case (naive Monte Carlo integration)\n",
    "\n",
    "- not to have to deal with pdf $ p(x) $ we can use **uniform distribution**  \n",
    "- then pdf $ p(x) $ is independent of $ x $ and can be treated as a constant  \n",
    "\n",
    "\n",
    "$$\n",
    "p(x) = \\left( \\int_{\\Omega} dx \\right) ^{-1} = \\frac{1}{V}\n",
    "$$\n",
    "\n",
    "- $ V $ is a measure of the set $ \\Omega $: length in one dimension, volume in two dimensions, etc.  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f =\\int_{\\Omega} f(x)\\,dx = E \\left[ Vf(\\tilde{X}) \\right] \\approx \\frac{V}{N} \\sum_{i=1}^{N} f(x_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Special even simpler case with unit hypercube\n",
    "\n",
    "- let $ \\Omega \\subset \\mathbb{R}^n $ be a unit hypercube denoted $ H_n $ in $ n $-dimensional space  \n",
    "- then $ V = 1 $  \n",
    "- integral is the same as simple average of the function computed on a random set of points uniformly distributed over the hypercube  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f =\\int_{H_n} f(x)\\,dx  \\approx \\frac{1}{N} \\sum_{i=1}^{N} f(x_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### One dimensional example\n",
    "\n",
    "- let $ \\Omega $ be an interval $ [a,b] \\subset \\mathbb{R} $ in one dimensional space  \n",
    "- then $ V = b-a $  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f =\\int_{a}^{b} f(x)\\,dx  \\approx \\frac{b-a}{N} \\sum_{i=1}^{N} f(x_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### One dimensional example visual\n",
    "\n",
    "<img src=\"_static/img/MCIntegration.png\" style=\"width:1000px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### General Monte Carlo integration algorithm\n",
    "\n",
    "1. sample $ N $ points $ x_1,\\cdots,x_N $ from distribution $ p(x) $ of $ \\tilde{X} $ on $ \\Omega $  \n",
    "1. approximate the expectation $ E \\left[ \\frac{f(\\tilde{X})}{p(\\tilde{X})} \\right] $ by the sample average  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f = \\int_{\\Omega} f(x)\\,dx \\approx \\frac{1}{N} \\sum_{i=1}^{N} \\frac{f(x_i)}{p(x_i)} = Q_f(N)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### General Monte Carlo integration algorithm (simple naive approach)\n",
    "\n",
    "1. sample $ N $ points $ x_1,\\cdots,x_N $ uniformly on $ \\Omega $  \n",
    "1. approximate the expectation $ E \\left[ V f(\\tilde{X})\\right] $ by the sample average  \n",
    "\n",
    "\n",
    "$$\n",
    "I_f = \\int_{\\Omega} f(x)\\,dx \\approx \\frac{V}{N} \\sum_{i=1}^{N} f(x_i), \\; V =\\int_{\\Omega} \\,dx\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = [12, 8]\n",
    "\n",
    "def mc_int_cube(f,ndims=1,N=1000):\n",
    "    '''Computes the integral of function f on a hypercube of dimension ndims\n",
    "       using Monte Carlo integration with N uniformly distributed points\n",
    "       Assume function f uses axis=0 for inputs, and can be vectorized in other axis\n",
    "       Return: value and standard error\n",
    "    '''\n",
    "    # generate uniform numbers on the hypercube\n",
    "    x = np.random.random(ndims*N).reshape(ndims,N)  # uniform random numbers in a matrix\n",
    "    y = f(x)                   # function value\n",
    "    Q = y.mean()               # sample average\n",
    "    seQ = y.std()/np.sqrt(N)   # standard error of sample average\n",
    "    return Q,seQ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Monte Carlo samples   :  1000\n",
      "Estimate (pi_hat)               :  3.088\n",
      "Standard error (pi_hat)         :  0.0530684087\n",
      "Approximation error (pi_hat-pi) :  -0.0535926536\n"
     ]
    }
   ],
   "source": [
    "# pi example from video 33 as two-dim integral\n",
    "# Approximate pi using 2-d Monte Carlo integration\n",
    "\n",
    "N=1000 # Number of Monte carlo Samples\n",
    "g = lambda x: (x[0,:]**2 + x[1,:]**2)<1  # indicator function to inegrate\n",
    "\n",
    "q,se = mc_int_cube(g,ndims=2,N=N)\n",
    "pi_hat = 4*q\n",
    "se_pi_hat = 4*se\n",
    "\n",
    "print('Number of Monte Carlo samples   : ', N);\n",
    "print('Estimate (pi_hat)               : ', pi_hat.round(10));\n",
    "print('Standard error (pi_hat)         : ', se_pi_hat.round(10));\n",
    "print('Approximation error (pi_hat-pi) : ', (pi_hat-np.pi).round(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Properties of Monte Carlo integral\n",
    "\n",
    "**Consistency**: Law of large numbers ensures that the sample average converge to the mean\n",
    "\n",
    "$$\n",
    "\\lim _{{N\\to \\infty }}Q_f(N)\n",
    "=\\lim _{{N\\to \\infty }}{\\frac{1}{N}}\\sum _{{i=1}}^{N}\\frac{f(x_i)}{p(x_i)}\n",
    "=E\\left[\\frac{f(\\tilde{x})}{p(\\tilde{x})}\\right] = \\int_{\\Omega} f(x)\\,dx = I_f\n",
    "$$\n",
    "\n",
    "**Assymptotic Normality**: By the central limit theorem we have\n",
    "\n",
    "$$\n",
    "\\sqrt{N}\\left(Q_f(N)-I_f \\right)\\ \\xrightarrow {d} \\ N\\left(0,\\sigma ^{2}\\right), \\; \\sigma^2= \\operatorname{Var}\\left[\\frac{f(\\tilde{x})}{p(\\tilde{x})}\\right]\n",
    "$$\n",
    "\n",
    "The standard error of $ Q_f(N) $ is then given by $ \\sigma_{Q_f(N)}=\\sigma \\big/ \\sqrt{N} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Standard error of Monte Carlo integral\n",
    "\n",
    "Given our estimate $ Q_f(N) $ of $ I_f $, we can obtain an unbiased estimate of $ \\sigma^2= \\operatorname{Var}\\left[\\frac{f(\\tilde{x})}{p(\\tilde{x})}\\right] $\n",
    "\n",
    "$$\n",
    "{\\hat{\\sigma}}^2_N=\\frac{1}{N-1}\\sum _{i=1}^N \\left(\\frac{f(x_i)}{p(x_i)}-Q_f(N)\\right)^2\n",
    "$$\n",
    "\n",
    "and the estimate of the standard error of $ Q_f(N) $\n",
    "\n",
    "$$\n",
    "{\\hat{\\sigma}}_{Q_f(N)}={\\hat{\\sigma}}_N   \\big/ \\sqrt{N}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convergence of Monte Carlo integrals\n",
    "\n",
    "The standard error of $ Q_f(N) $:\n",
    "\n",
    "- is given by $ \\sigma_{Q_f(N)}=\\sigma  \\big/ \\sqrt{N} $  \n",
    "- can be estimated by $ {\\hat{\\sigma}}_{Q_f(N)}={\\hat{\\sigma}}_N   \\big/ \\sqrt{N} $  \n",
    "\n",
    "\n",
    "**Decreases with the standard parametric rate** $ \\sqrt{N} $\n",
    "\n",
    "- doubling of precision requires 4 time as many random points  \n",
    "- but does not depend on the dimensionality of the integral, $ \\Omega $ can be high dimensional  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True value (pi)                  : 3.1415926536\n",
      "Average estimate across all runs : 3.13976\n",
      "Mean bias                        : -0.0018326536\n",
      "Average std err across all runs  : 0.0519299844\n",
      "Std dev of bias                  : 0.0533409261\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAHiCAYAAADh4aRaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7jsdV0v8PcnNmgpCMiGEIWlhSX1iNaOLM1jYabhhec8VpolGR2yzDrHOrm9lNnNbZ3T/aKcLNG8YCaPHLfHJJLMShQKL4SG4U4JFARUUFPB7/nj99s6Ltbaa2avmbUW+/t6Pc88a+Z3/cx3frPmPd/5zm+qtRYAAOjBV2x2AQAAsFGEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCL6xTVb2oqn5hTts6vqpuqaqDxtsXVdWPzWPb4/b+X1WdMa/tzbDfX62qj1XVRzZ63+xbVT2sqq7e7DoWaZ7P0Sn29R1V9f6N2Nc0avBnVXVTVb1jP9bflP8ZsEjCL+xDVe2pqs9U1c1V9fGq+oeqempVffG501p7amvtV6bc1sP3tUxr7UOttbu21m6bQ+2/VFV/vmz7j2qtnbPebc9Yx72S/GySk1prX73C/EOq6rVj+7Sqetiy+VVVL6yqG8bLb1RVTcx/QFVdWlWfHv8+YNp1l+3nYeP+X7ds+snj9IvW1xJJVS2N29q2jm0cMj62V1bVp8Z2+9OqWlpvfQeqaZ+j+2N8PL92Yl9/11r7ukXsaz89JMl3J7lna+2UWVfejP8ZsGjCL6ztMa21Q5OckGRXkmcmecm8d7KeQLTFnZDkhtbadftY5m1JfijJSj3DZyU5PcnJSe6f5NFJfjwZgmCS1yf58yRHJDknyevH6ftcdxXXJ/n2qrr7xLQzkvzrPtbZaK9N8tgkP5jkbhnu26VJTp11Q3eUY+6OUucWdUKSPa21T212IbBltNZcXFxWuSTZk+Thy6adkuQLSb5xvP3SJL86Xj8qyRuSfDzJjUn+LsObzJeP63wmyS1Jfj7JUpKW5MwkH0ry1olp28btXZTkBUnekeQTGYLekeO8hyW5eqV6kzwyyeeSfH7c37smtvdj4/WvSPLcJP+e5LokL0tyt3He3jrOGGv7WJLn7KOd7jauf/24veeO23/4eJ+/MNbx0jXa++okD1s27R+SnDVx+8wkbx+vPyLJfySpifkfSvLItdZdYd8PG/f/oiRPG6cdNE77xSQXTSz77UneOT4m70zy7RPzLkryK0n+PsnNSd6c5KiJ2trYFrck+bZx+o8muSLJTUn+KskJq9S4tz3vtY82fMq4rZuTXJXkx1e4j8/M8Ebj5cuPoyT3G+/Dx5NcnuSx69zXs8fjZ0+SJ03Mf+nY1heM6//t5P0e2+lpSa5M8sF9tXuSI8d9PWa8fdckH0jy5BWeo3vr+vkMx/21Gd4gfW+GNzk3Jnn2suf7P47tcW2SP0hyyDjvrWOdnxofzx+YpT3Huv4wye6xDS5O8jXjvEry22ONn0jy7oz/c1Z4HO6R5Pyx9g8k+W8Tx/t/JrltrO/5K6z7IxmO1d8f9/O+JKcuO55/bNH/a11cNvKi5xdm1Fp7R4YXz+9YYfbPjvO2Jzkmwwt/a639cIbg85g2DGv4jYl1/kuGF8jvWWWXT84Qju6R5NYkvzdFjW9K8utJzh33d/IKi/3IePnOJPfJEBj+YNkyD0nydRl6FX+xqu63yi5/P0MAvs94f56c5Cmttb9O8qgk14x1/Mhata/gG5K8a+L2u8Zpe+e9u7U2+Tvt7142f7V1V/Oysf5keEwuT3LN3plVdWSGsPJ7Se6e5LeS7F7WW/yDGYLh0UkOSfJz4/SHjn8PH9vjH6vq9AzHyX/NcNz8XZJXrVLbw5O8o7X24X3Uf12GHu7Dxhp+u6q+aWL+V2cIiydk6Bn/oqo6OMn/zRDYj07y9CSvqKrVPsafZl9HJTkuwxups5dt60kZ3igcleSyJK9Ytv3Tk3xrkpP21e6ttRszPEf+T1UdnSE0XtZae9kqdX91kjuPdf1ikv+T4ZOHb87wvP7FqrrPuOxtSf7HWOO3ZXgu/GSStNb2Pp4nj4/nuZM7mbI9n5jk+Rk+ufhAkl8bpz8iw/Fy3ySHZwjWN6xyf16V4f/OPZI8PsmvV9WprbWXJHlqkn8c63veKut/a4Y3L0cleV6S143tDQck4Rf2zzUZAsRyn09ybIYerM+3YfxfW2G5Sb/UWvtUa+0zq8x/eWvtvW342PIXknz/3i/ErdOTkvxWa+2q1totSZ6V5AnLPmJ+fmvtM621d2UIjrcL0WMtP5DkWa21m1tre5L87yQ/PIcakyGUf2Li9ieS3HUcu7t83t75h06x7opaa/+Q5MgxoDw5QxiedFqSK1trL2+t3dpae1WG3rLHTCzzZ621fx0f09ckeUBW9+NJXtBau6K1dmuGNy0PqKoTVlj27hl6H1fVWtvdWvu3NvjbDMFr8o3aF5I8r7X22RWOuQdlaLNdrbXPtdb+JsMnGU/cz30lyS+M+/rbDOH1+yfm7W6tvbW19tkkz0nybeMY8b1e0Fq7caxzn+3eWntzkr9IcuG47L6Gt3w+ya+11j6f5NUZQt/vjsfv5Rne8Nx/3O6lrbW3j/vck+TFGd7gTWOa9nxda+0d42P/inzpWPl8huP46zN8snFFa+12j/3YXg9J8szW2n+21i5L8ieZ7fl3XZLfGf9nnZvk/RnaEA5Iwi/sn+MyfMS43G9m6L15c1VdVVU7p9jWvnrxls//9yQHZ3ixXq97jNub3Pa2DD3We02Owf10hhfy5Y7K0Lu5fFvHzaHGZPi49rCJ24cluWV8U7F83t75N0+x7r68PMlPZegVP2/ZvOXtltz+/k7TbnudkOR3xy9U7h0uU1m5/W7I8OZqVVX1qKp6e1XdOG7ve/Plx8v1rbX/XGX1eyT5cGvtCxPTVn0sp9jXTe3Lx5r++7iPvb54bI9vwG5cbX6ma/ezk3xjhjcfq/WSJsMY9L1fKt37BuCjE/M/k/Exq6r7VtUbquojVfXJDG9Opn3+TdOeKx4rY1D+gwzDIj5aVWdX1fJjfe8+bmyt3Twxbdbn338se04sf5zggCL8woyq6lsyvLC8bfm8sefoZ1tr98nQI/WMqtr7RaTVAtdaQWyyJ+z4DD1CH8swzvCrJuo6KMPH5tNu95oMwWty27fmy0PAND421rR8W/8x43ZWc3m+vMf55HHa3nn3X9aTe/9l81dbd19enuGj7Te21j69bN7ydkumv78rPSYfzjBW9vCJy1eOPdDL/XWSU6rqnittvKrulOQvk/yvJMe01g5P8sYMYXpfNex1TZJ7TZ7NJKvctyn3dURV3WXZtq6ZuP3FY7uq7prh05TJ+ZO17rPdx+P/xRl66n9i8gwM6/THGXqYT2ytHZZhiMqqnxwsM3V7rqS19nuttW/OMFTnvkn+5yr7OLKqDp2YNuvz77hlz6HljxMcUIRfmFJVHVZVj87wMemft9bes8Iyj66qrx1fSD6ZYbzg3h6mj2YYEzurH6qqk6rqq5L8cpLXjr1W/5rkzlV12ji28LlJ7jSx3keTLC174Z30qiT/o6ruPQaPvWOEb52luLGW1yT5tao6dPy4/hkZzsAwlaq6U1Xdebx5SFXdeeLF+GUZ3kQcV1X3yDCu+qXjvIsytO9Pj9v4qXH630yx7r7u0wczfLT9nBVmvzHJfavqB6tqW1X9QJKTMnycvZbrMww7mDwOXpTkWVX1DUlSVXerqu9bpa6/zvAFsfOq6pvH/R9aw+n3fjRDD/ydxv3cWlWPyjB2dFoXZ3hT9fNVdXANp517TIZjfrlp9/X88fRs35FhfPBfTMz73qp6yHh2jl9JcvE+xjOv1e7PHv/+aIZA/rI5DQ86NMNz+Zaq+vokP7Fs/r6e17O055epqm+pqm8dn9ufype+uPZlxvb6hyQvGJ8398/wRbfl46f35egMz6GDx2PvfhnaGw5Iwi+s7f9W1c0Zeuiek+GLNk9ZZdkTM/TO3ZLhG+J/1Fq7aJz3giTPHT/e/rlV1l/JyzMEto9k+JLOTydJa+0TGXon/yRDL8+nMnzpZa+9IeOGqvqnFbb7p+O235rkgxleXJ8+Q12Tnj7u/6oMPeKvHLc/rfdn+Kj5uAxnO/hMvtTL9+IMXxp6T5L3Zhg3+uIkaa19LsOXop6c4dv0P5rk9HH6PtddS2vtba212/V+jR+nPzpDkL4hw1kDHt1a+9gU2/x0hi80/f14HDyotXZekhcmefX4sfp7M3xJcDWPzxBMzs0whvm9SXYk+evxo++fzvBm5KYMX7w7f5r7O9b3uQynUXtUhh79P8pwxoT3rbDsNPv6yDjvmgxh7KnLtvXKDF+wujHDl82etI/aVm33qvrmDG+4njy+GXthhl7jaYYdreXnMty3mzN8Me7cZfN/Kck54+M5OZ55pvZcwWHj/m7KMAzhhgyhfiVPzHCGlmsyDNN5Xmvtgin2sdfFGf53fSzD8fn4NYaNwB1arT30DQBmM/Zy/nlrbbUhGi/NcEqw525kXXy5qvqRDKcye8hm1wIbRc8vAADdEH4BAOiGYQ8AAHRDzy8AAN0QfgEA6Ma2tReZn6OOOqotLS1t5C4BAOjMpZde+rHW2vaV5m1o+F1aWsoll1yykbsEAKAzVbX859C/yLAHAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbmzb7AIA2LelnbtnWn7PrtMWVAnAHZ+eXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6MZU4beqDq+q11bV+6rqiqr6tqo6sqouqKorx79HLLpYAABYj2l7fn83yZtaa1+f5OQkVyTZmeTC1tqJSS4cbwMAwJa1ZvitqsOSPDTJS5Kktfa51trHkzwuyTnjYuckOX1RRQIAwDxsm2KZ+yS5PsmfVdXJSS5N8jNJjmmtXZskrbVrq+rolVauqrOSnJUkxx9//FyKBliUpZ27Z1p+z67TFrp9AOZrmmEP25J8U5I/bq09MMmnMsMQh9ba2a21Ha21Hdu3b9/PMgEAYP2mCb9XJ7m6tXbxePu1GcLwR6vq2CQZ/163mBIBAGA+1gy/rbWPJPlwVX3dOOnUJP+S5PwkZ4zTzkjy+oVUCAAAczLNmN8keXqSV1TVIUmuSvKUDMH5NVV1ZpIPJfm+xZQIAADzMVX4ba1dlmTHCrNOnW85AACwOH7hDQCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDoxrbNLgCAA9vSzt0zLb9n12kLqgRAzy8AAB0RfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDoxrbNLgBgkZZ27t7sEgDYQvT8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohlOdAayDU6kB3LHo+QUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A2nOgM4wMx6+rU9u05bUCUAW4+eXwAAujFVz29V7Ulyc5LbktzaWttRVUcmOTfJUpI9Sb6/tXbTYsoEAID1m6Xn9ztbaw9ore0Yb+9McmFr7cQkF463AQBgy1rPsIfHJTlnvH5OktPXXw4AACzOtOG3JXlzVV1aVWeN045prV2bJOPfoxdRIAAAzMu0Z3t4cGvtmqo6OskFVfW+aXcwhuWzkuT444/fjxIB2EpmPZsEwFYyVc9va+2a8e91Sc5LckqSj1bVsUky/r1ulXXPbq3taK3t2L59+3yqBgCA/bBm+K2qu1TVoXuvJ3lEkvcmOT/JGeNiZyR5/aKKBACAeZhm2MMxSc6rqr3Lv7K19qaqemeS11TVmUk+lOT7FlcmAACs35rht7V2VZKTV5h+Q5JTF1EUAAAsgl94AwCgG8IvAADdmPZUZwBbgtNsAbAeen4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0I1tm10AAJtraefuzS4BYMPo+QUAoBvCLwAA3RB+AQDohvALAEA3hF8AALrhbA/ApnKmAQA2kp5fAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuTB1+q+qgqvrnqnrDePveVXVxVV1ZVedW1SGLKxMAANZvlp7fn0lyxcTtFyb57dbaiUluSnLmPAsDAIB5myr8VtU9k5yW5E/G25Xku5K8dlzknCSnL6JAAACYl2l7fn8nyc8n+cJ4++5JPt5au3W8fXWS4+ZcGwAAzNWa4beqHp3kutbapZOTV1i0rbL+WVV1SVVdcv311+9nmQAAsH7T9Pw+OMljq2pPkldnGO7wO0kOr6pt4zL3THLNSiu31s5ure1ore3Yvn37HEoGAID9s2b4ba09q7V2z9baUpInJPmb1tqTkrwlyePHxc5I8vqFVQkAAHOwnvP8PjPJM6rqAxnGAL9kPiUBAMBibFt7kS9prV2U5KLx+lVJTpl/SQAAsBh+4Q0AgG7M1PMLsC9LO3dvdgkAsE96fgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdMOpzgDYUjbilHl7dp228H0AW5OeXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAurFtswsANs7Szt0zLb9n12kLqgQ2l+cC9EvPLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBurBl+q+rOVfWOqnpXVV1eVc8fp9+7qi6uqiur6tyqOmTx5QIAwP6bpuf3s0m+q7V2cpIHJHlkVT0oyQuT/HZr7cQkNyU5c3FlAgDA+q0ZftvglvHmweOlJfmuJK8dp5+T5PSFVAgAAHMy1Zjfqjqoqi5Lcl2SC5L8W5KPt9ZuHRe5OslxiykRAADmY9s0C7XWbkvygKo6PMl5Se630mIrrVtVZyU5K0mOP/74/SwTWMnSzt2bXQIA3KHMdLaH1trHk1yU5EFJDq+qveH5nkmuWWWds1trO1prO7Zv376eWgEAYF2mOdvD9rHHN1X1lUkenuSKJG9J8vhxsTOSvH5RRQIAwDxMM+zh2CTnVNVBGcLya1prb6iqf0ny6qr61ST/nOQlC6wTAADWbc3w21p7d5IHrjD9qiSnLKIoAABYBL/wBgBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOjGNL/wBnRqaefuzS4BAOZKzy8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBubNvsAgDgQLO0c/dMy+/ZddqCKgGW0/MLAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbTnUGAGuY9dRlwNal5xcAgG4IvwAAdGPN8FtV96qqt1TVFVV1eVX9zDj9yKq6oKquHP8esfhyAQBg/03T83trkp9trd0vyYOSPK2qTkqyM8mFrbUTk1w43gYAgC1rzfDbWru2tfZP4/Wbk1yR5Lgkj0tyzrjYOUlOX1SRAAAwDzON+a2qpSQPTHJxkmNaa9cmQ0BOcvQq65xVVZdU1SXXX3/9+qoFAIB1mDr8VtVdk/xlkv/eWvvktOu11s5ure1ore3Yvn37/tQIAABzMVX4raqDMwTfV7TWXjdO/mhVHTvOPzbJdYspEQAA5mOasz1UkpckuaK19lsTs85PcsZ4/Ywkr59/eQAAMD/T/MLbg5P8cJL3VNVl47RnJ9mV5DVVdWaSDyX5vsWUCAAA87Fm+G2tvS1JrTL71PmWAwAAi+MX3gAA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALoh/AIA0A3hFwCAbmzb7AJgsyzt3D3T8nt2nbagSgCAjaLnFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IazPcACOaMEsCj+v8D+0fMLAEA3hF8AALoh/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG5s2+wCgC9Z2rl7s0sANoHnPmwcPb8AAHRD+AUAoBvCLwAA3RB+AQDohvALAEA3hF8AALrhVGcAwO3Mevq1PbtOW1AlMF96fgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdWDP8VtWfVtV1VfXeiWlHVtUFVXXl+PeIxZYJAADrN03P70uTPHLZtJ1JLmytnZjkwvE2AABsaWuG39baW5PcuGzy45KcM14/J8npc64LAADmbn/H/B7TWrs2Sca/R8+vJAAAWIyFf+Gtqs6qqkuq6pLrr79+0bsDAIBV7W/4/WhVHZsk49/rVluwtXZ2a21Ha23H9u3b93N3AACwfvsbfs9PcsZ4/Ywkr59POQAAsDjTnOrsVUn+McnXVdXVVXVmkl1Jvruqrkzy3eNtAADY0rattUBr7YmrzDp1zrUAAMBC+YU3AAC6IfwCANCNNYc9AIOlnbs3uwQAYJ30/AIA0A3hFwCAbgi/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG9s2uwAAYPGWdu7e7BJgS9DzCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdMPZHtiyZv1m8p5dpy2oEgDmzf94NoueXwAAuiH8AgDQDeEXAIBuCL8AAHRD+AUAoBvCLwAA3XCqMw4Ys542BwDoj55fAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBuCL8AAHRj22YXAADc8S3t3L3ltr9n12kLqIQ7Oj2/AAB0Q/gFAKAbwi8AAN0QfgEA6IbwCwBAN4RfAAC6IfwCANAN4RcAgG4IvwAAdEP4BQCgG8IvAADdEH4BAOjGts0uYCta2rl74fvYs+u0hW5/1vuw6HqSjWlXANhrK74WLtL+vM7e0e/z/tDzCwBAN4RfAAC6sa7wW1WPrKr3V9UHqmrnvIoCAIBF2O/wW1UHJfnDJI9KclKSJ1bVSfMqDAAA5m09Pb+nJPlAa+2q1trnkrw6yePmUxYAAMzfesLvcUk+PHH76nEaAABsSes51VmtMK3dbqGqs5KcNd68parev459bqSjknxsURuvFy5qy/tnDvUstL0OQNprNtprNtprdtpsNgdkey3wtXnLttdWyyOjebTXCavNWE/4vTrJvSZu3zPJNcsXaq2dneTsdexnU1TVJa21HZtdxx2F9pqN9pqN9pqN9pqdNpuN9pqN9prNottrPcMe3pnkxKq6d1UdkuQJSc6fT1kAADB/+93z21q7tap+KslfJTkoyZ+21i6fW2UAADBn6/p549baG5O8cU61bDV3uKEam0x7zUZ7zUZ7zUZ7zU6bzUZ7zUZ7zWah7VWt3e47agAAcEDy88YAAHSju/BbVfeqqrdU1RVVdXlV/cwKyxxRVedV1bur6h1V9Y0T8/ZU1Xuq6rKqumRjq994VXXnsQ3eNbbX81dY5k5Vde74M9cXV9XSxLxnjdPfX1Xfs5G1b4b1tFdVLVXVZ8Zj67KqetFG17/Rpmyvh1bVP1XVrVX1+GXzzqiqK8fLGRtX+eaYQ3vdNnF8HfBfUJ6yvZ5RVf8y/r+/sKpOmJjn+Lr9MvtqL8fX7Zd56kRmeNvkL+F6fZy+veb++tha6+qS5Ngk3zRePzTJvyY5adkyv5nkeeP1r09y4cS8PUmO2uz7sYHtVUnuOl4/OMnFSR60bJmfTPKi8foTkpw7Xj8pybuS3CnJvZP8W5KDNvs+beH2Wkry3s2+D1uwvZaS3D/Jy5I8fmL6kUmuGv8eMV4/YrPv01Ztr3HeLZt9H7Zge31nkq8ar//ExPPR8TVDezm+Vm2vwyauPzbJm8brXh9na6+5vj521/PbWru2tfZP4/Wbk1yR2/8y3UlJLhyXeV+Spao6ZkML3SLa4Jbx5sHjZflA8cclOTz33/IAAAQfSURBVGe8/tokp1ZVjdNf3Vr7bGvtg0k+kOFnsQ9Y62yv7kzTXq21Pa21dyf5wrLVvyfJBa21G1trNyW5IMkjF13zZlpne3VnyvZ6S2vt0+PNt2c4Z33i+Jq1vbozZXt9cuLmXSbme32crb3mqrvwO2n8uPmBGd59THpXkv86LnNKhl8J2fsEb0neXFWX1vDrdQe8qjqoqi5Lcl2GF4Pl7fXFn7purd2a5BNJ7p5OfwJ7He2VJPeuqn+uqr+tqu/YsKI30RTttRrH12ztlSR3rqpLqurtVXX6gkrcUmZsrzOT/L/xuuNrtvZKHF8rtldVPa2q/i3JbyT56XGy42u29krm+PrYbfitqrsm+csk/33ZO40k2ZXkiPEBenqSf05y6zjvwa21b0ryqCRPq6qHblTNm6W1dltr7QEZ3gCcUhNjoEer/dT1VD+BfaBZR3tdm+T41toDkzwjySur6rDFVrv5pmiv1Ti+ZmuvZDi+diT5wSS/U1Vfs5Ait5Bp26uqfijJjgzD3hLH16ztlTi+Vmyv1toftta+Jskzkzx3nOz4mq295vr62GX4raqDMwTfV7TWXrd8fmvtk621p4wP0JOTbE/ywXHeNePf65KclwP8Y4pJrbWPJ7kot//o74s/dV1V25LcLcmNmfInsA9Us7bX+PHXDeO6l2YYA3bfDSt4k+2jvVbj+JqtvSb/f101rvvARdS2Fe2rvarq4Umek+SxrbXPjpMdX7O1l+Nr7efjq5Ps7RF3fM3QXvN+fewu/I5jK1+S5IrW2m+tsszhNfxkc5L8WJK3ttY+WVV3qapDx2XukuQRSd67EXVvlqraXlWHj9e/MsnDk7xv2WLnJ9n7TejHJ/mb1lobpz+hhrMb3DvJiUnesTGVb471tNe47kHjuvfJ0F5XbUzlm2PK9lrNXyV5RA1nZzkiw/PxrxZT6dawnvYa2+lO4/Wjkjw4yb8sqtatYJr2qqoHJnlxhiB33cQsx9cM7eX4WrW9Tpy4eVqSK8frXh9naK95vz6u6xfe7qAenOSHk7xnHNaQJM9OcnyStNZelOR+SV5WVbdlePKeOS53TJLzxu8mbUvyytbamzaw9s1wbJJzxoPuK5K8prX2hqr65SSXtNbOz/Bm4uVV9YEMPb5PSJLW2uVV9ZoMbXhrkqe11m7blHuxcfa7vZI8NMkvV9WtSW5L8tTW2o0bfxc21JrtVVXfkuFTliOSPKaqnt9a+4bW2o1V9StJ3jlu65e11+rtleH/2our6gvjurtaawd0OMl0z8ffTHLXJH8x/m//UGvtsY6v2dorjq/V2uunxp7yzye5KWPHh9fH2dorc3599AtvAAB0o7thDwAA9Ev4BQCgG8IvAADdEH4BAOiG8AsAQDeEXwAAuiH8AgDQDeEXAIBu/H951TSwopBiYAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# distribution of Monte Carlo integral\n",
    "\n",
    "N=1000; # number of Monte Carlo samples used to simulate the integral\n",
    "S=1000; # number of runs to generate the distribution of estimates\n",
    "\n",
    "qs = np.empty(S,dtype=float)\n",
    "ses = np.empty(S,dtype=float)\n",
    "\n",
    "for i in range(S):\n",
    "    q,se = mc_int_cube(g,ndims=2,N=N)\n",
    "    qs[i] = 4*q\n",
    "    ses[i] = 4*se\n",
    "\n",
    "plt.hist(qs,bins=50,range=(np.pi-.2, np.pi+.2))\n",
    "plt.title('Distribution of %d Monte Carlo approximations of pi'%S)\n",
    "\n",
    "print('True value (pi)                  :', np.round(np.pi,10))\n",
    "print('Average estimate across all runs :', qs.mean().round(10))\n",
    "print('Mean bias                        :', np.mean(qs-np.pi).round(10))\n",
    "print('Average std err across all runs  :', ses.mean().round(10))\n",
    "print('Std dev of bias                  :', np.std(qs-np.pi).round(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Further learning resources\n",
    "\n",
    "- SciPy docs [https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html](https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html)  \n",
    "- [https://en.wikipedia.org/wiki/Gaussian_quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature)  \n",
    "- Monte Carlo integration [https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration](https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration)  \n",
    "- Useful library for Monte Carlo methods [https://chaospy.readthedocs.io/en/master/index.html](https://chaospy.readthedocs.io/en/master/index.html)  "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "date": 1612589586.2915368,
  "download_nb": false,
  "filename": "34_integration.rst",
  "filename_with_path": "34_integration",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "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.7.6"
  },
  "title": "Foundations of Computational Economics #34"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}