{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "execute:\n",
    "  enabled: false\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-1 fw=all -->\n",
    "\n",
    "# Integral Calculus\n",
    "\n",
    "Differentiation answered a local question: how does a function change when we\n",
    "nudge its input? Integration answers a global one: how much of something is\n",
    "there in total, whether the area under a curve, the volume under a surface, or\n",
    "the probability under a density. The two look unrelated, yet the *fundamental\n",
    "theorem of calculus* joins them into a single subject: integration is\n",
    "differentiation run backwards. That one fact is what makes integrals\n",
    "computable at all, and it is the reason a deep-learning reader needs them, since\n",
    "every continuous probability is an integral and every expectation is an integral\n",
    "average (that section).\n",
    "\n",
    "We will not need the full machinery of a calculus course. This section builds\n",
    "what an integral *is* (a limit of sums, extended to unbounded domains), the\n",
    "fundamental theorem that lets us *compute* integrals together with the\n",
    "integration-by-parts and change-of-variables rules it yields, and the\n",
    "multiple-integral toolkit (Fubini's theorem and the $n$-dimensional change of\n",
    "variables) that powers the Gaussian normalizer and, later, normalizing flows\n",
    "(that section). The last part turns to\n",
    "probability: densities, expectations, Monte Carlo estimation, and\n",
    "differentiation under the integral sign."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-imports"
   },
   "source": [
    "%matplotlib inline\n",
    "from d2l import jax as d2l\n",
    "import numpy as np\n",
    "import jax\n",
    "from jax import numpy as jnp"
   ],
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-2 fw=all -->\n",
    "\n",
    "## The Definite Integral\n",
    "\n",
    "Take a function $f$ and ask for the area trapped between its graph and the\n",
    "$x$-axis over an interval $[a, b]$. The area under an entire curve is usually\n",
    "infinite or undefined, so we always pin it between two endpoints. The recipe is\n",
    "the same one Archimedes used: chop $[a,b]$ into $N$ thin vertical slices of width\n",
    "$\\epsilon = (b-a)/N$, approximate each slice by a rectangle of height $f(x_i)$,\n",
    "and add them up. As the slices shrink, the staircase of rectangles squeezes onto\n",
    "the curve and the sum approaches a definite number, shown in\n",
    "the figure.\n",
    "\n",
    "![Left to right: the area under $f$ from $a$ to $b$ approximated by rectangles of shrinking width $\\epsilon$. As $\\epsilon\\to 0$ the rectangle sum $\\sum_i \\epsilon\\,f(x_i)$ converges to the definite integral, the exact signed area.](../img/mdl-cal-riemann.svg)\n",
    "\n",
    "We *define* the **definite integral** as this limit and write it with the\n",
    "elongated-$S$ symbol,\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x)\\,dx = \\lim_{\\epsilon\\to 0}\\ \\sum_i \\epsilon\\, f(x_i),\n",
    "$$\n",
    "\n",
    "a continuous analogue of a sum: $\\sum$ becomes $\\int$, the spacing $\\epsilon$\n",
    "becomes $dx$, and the term $f(x_i)$ becomes $f(x)$. The inner variable is a dummy,\n",
    "exactly like a summation index, so $\\int_a^b f(x)\\,dx = \\int_a^b f(z)\\,dz$. We\n",
    "grant the following facts about this limit once and for all\n",
    "[@Folland.1999]. For bounded, piecewise-continuous $f$ (continuous except\n",
    "at finitely many jump points), which covers everything we meet in machine\n",
    "learning, the limit exists and does not depend on how the slices are chosen.\n",
    "The integral it defines is *monotone*, so $f \\le g$ implies\n",
    "$\\int_a^b f\\,dx \\le \\int_a^b g\\,dx$, and *additive over subintervals*, so\n",
    "$\\int_a^c f\\,dx = \\int_a^b f\\,dx + \\int_b^c f\\,dx$ for $a\\le b\\le c$; both are\n",
    "plausible from the rectangle picture, and we use both below. (Some hypothesis\n",
    "on $f$ is genuinely needed: for the function that is $1$ on the rationals and\n",
    "$0$ elsewhere, the rectangle sums never settle on any value at all. Lebesgue's\n",
    "integral, built on measure theory rather than on slicing the axis, handles even\n",
    "that function, assigning it integral $0$ because the rationals have measure\n",
    "zero, and it is the foundation probability theory ultimately stands on; the\n",
    "Riemann picture suffices for every computation in this book.)\n",
    "\n",
    "The definition is not yet a *computation*: only the simplest\n",
    "integrands (a line, $\\int_a^b x\\,dx$) succumb to summing the rectangles by hand.\n",
    "The next section gives the tool that computes the rest. First let us check that\n",
    "the limit the equation is real by watching it converge. Take\n",
    "$\\int_0^2 \\tfrac{x}{1+x^2}\\,dx$, whose exact value is $\\tfrac12\\log 5$ (we will\n",
    "see why shortly), and refine the partition: the rectangle sum should march toward\n",
    "the truth, with the error shrinking roughly in proportion to $\\epsilon$ for this\n",
    "left-endpoint rule."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-riemann-converge"
   },
   "source": [
    "import numpy as onp\n",
    "def riemann(eps, a=0., b=2.):\n",
    "    x = onp.arange(a, b, eps)\n",
    "    return float(onp.sum(eps * x / (1 + x**2)))\n",
    "\n",
    "truth = float(onp.log(onp.array(5.)) / 2)\n",
    "for eps in [0.5, 0.1, 0.05, 0.01, 0.001]:\n",
    "    approx = riemann(eps)\n",
    "    print(f'eps={eps:<6} sum={approx:.6f}  error={abs(approx - truth):.6f}')\n",
    "print(f'truth = (1/2) log 5 = {truth:.6f}')"
   ],
   "execution_count": 2,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "eps=0.5    sum=0.680769  error=0.123950\neps=0.1    sum=0.783785  error=0.020934\neps=0.05   sum=0.794486  error=0.010233\neps=0.01   sum=0.802710  error=0.002009\neps=0.001  sum=0.804519  error=0.000200\ntruth = (1/2) log 5 = 0.804719\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-3 fw=all -->\n",
    "\n",
    "The error falls by about a factor of ten each time we cut $\\epsilon$ by ten,\n",
    "confirming the first-order convergence of the left-rule. Refining the partition\n",
    "*works*, but it is slow and gives no closed form. We need a better idea.\n",
    "\n",
    "## The Fundamental Theorem of Calculus\n",
    "\n",
    "The breakthrough is to stop treating the integral as a fixed number and instead\n",
    "let its upper limit *move*. Define the **area-so-far function**\n",
    "\n",
    "$$\n",
    "F(x) = \\int_a^x f(y)\\,dy,\n",
    "$$\n",
    "\n",
    "the signed area accumulated from the left endpoint $a$ out to $x$. This single\n",
    "function captures *every* definite integral at once: the area over any\n",
    "sub-interval is a difference of accumulated areas,\n",
    "\n",
    "$$\n",
    "\\int_c^b f(x)\\,dx = F(b) - F(c),\n",
    "$$\n",
    "\n",
    "since the area out to $b$ minus the area out to $c$ leaves exactly the area\n",
    "between them. the figure draws this subtraction,\n",
    "accumulating from a base point $c$ at the left edge of the plot.\n",
    "\n",
    "![The area under $f$ from $a$ to $b$ as the difference of two areas accumulated from a base point $c$ at the plot's left edge: $\\int_a^b f\\,dx = \\int_c^b f\\,dx - \\int_c^a f\\,dx$. Every definite integral is a difference of two values of a single accumulation function.](../img/mdl-cal-sub-area.svg)\n",
    "\n",
    "So knowing $F$ solves the whole problem. The *fundamental theorem of calculus*\n",
    "is the fact that $F$ is determined by a derivative.\n",
    "\n",
    "**Theorem (Fundamental theorem of calculus).** *Let $f$ be continuous and let\n",
    "$F(x)=\\int_a^x f(y)\\,dy$. Then $F$ is differentiable and*\n",
    "\n",
    "$$\n",
    "\\frac{dF}{dx}(x) = f(x).\n",
    "$$\n",
    "\n",
    "**Proof.** Differentiate $F$ from the definition: nudge $x$ by $\\epsilon$ and\n",
    "read off how much new area appears. By the equation and additivity of\n",
    "the integral, the increment is a single thin sliver,\n",
    "\n",
    "$$\n",
    "F(x+\\epsilon) - F(x) = \\int_x^{x+\\epsilon} f(y)\\,dy .\n",
    "$$\n",
    "\n",
    "The sliver sits over an interval of width $\\epsilon$. Write\n",
    "$m_\\epsilon = \\inf_{[x,x+\\epsilon]} f$ and\n",
    "$M_\\epsilon = \\sup_{[x,x+\\epsilon]} f$ for the tightest constant bounds on $f$\n",
    "there. Comparing $f$ with those two constants, monotonicity of the integral\n",
    "pins the sliver between two rectangles of width $\\epsilon$:\n",
    "\n",
    "$$\n",
    "\\epsilon\\, m_\\epsilon\n",
    "\\ \\le\\ F(x+\\epsilon) - F(x)\\ \\le\\\n",
    "\\epsilon\\, M_\\epsilon .\n",
    "$$\n",
    "\n",
    "Divide by $\\epsilon$. Because $f$ is continuous at $x$, for any $\\delta>0$ the\n",
    "whole interval $[x,x+\\epsilon]$ eventually lies where $|f(y)-f(x)|<\\delta$, so\n",
    "$m_\\epsilon$ and $M_\\epsilon$ both converge to the single value\n",
    "$f(x)$ as $\\epsilon \\to 0$, and the squeeze forces\n",
    "$\\lim_{\\epsilon\\to 0}\\frac{F(x+\\epsilon)-F(x)}{\\epsilon}=f(x)$ (for\n",
    "$\\epsilon<0$ the sliver sits to the left of $x$ and the same bounds apply, with\n",
    "both sign flips canceling in the quotient). That limit is\n",
    "exactly $F'(x)$. $\\blacksquare$\n",
    "\n",
    "The sliver argument is the whole story: the rate at which accumulated area grows\n",
    "is just the current height of the curve. This *reverses* the problem. Finding\n",
    "areas, hard on its own, becomes the search for an **antiderivative**, a function\n",
    "whose derivative is $f$, which we can read straight off the derivative table of\n",
    "that section. If $G'=f$ then $G$ and $F$ differ by a\n",
    "constant (they have the same derivative), and since that constant cancels in any\n",
    "difference,\n",
    "\n",
    "$$\n",
    "\\int_a^b f(x)\\,dx = G(b) - G(a)\n",
    "$$\n",
    "\n",
    "for *any* antiderivative $G$. The arbitrary constant of integration is real but\n",
    "irrelevant to definite integrals.\n",
    "\n",
    "This turns hard sums into easy lookups. The derivative of $x^n$ is $nx^{n-1}$, so\n",
    "running it backwards,\n",
    "\n",
    "$$\n",
    "\\int_0^{x} n\\,y^{n-1}\\,dy = x^n - 0^n = x^n ,\n",
    "$$\n",
    "\n",
    "and since $e^x$ is its own derivative,\n",
    "\n",
    "$$\n",
    "\\int_0^{x} e^{t}\\,dt = e^{x} - e^{0} = e^x - 1 .\n",
    "$$\n",
    "\n",
    "Every integration rule is a differentiation rule read in reverse, which is why\n",
    "the running example obeyed $\\int_0^2 x/(1+x^2)\\,dx=\\tfrac12\\log 5$: the\n",
    "antiderivative of $x/(1+x^2)$ is $\\tfrac12\\log(1+x^2)$, evaluated at the\n",
    "endpoints. Let us confirm the theorem numerically by checking that the\n",
    "area-so-far function $F$ has derivative $f$: compute $F$ as a cumulative Riemann\n",
    "sum, finite-difference it, and compare against $f$ itself."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-ftc-check"
   },
   "source": [
    "import numpy as onp\n",
    "f = lambda t: t / (1 + t**2)\n",
    "G = lambda t: onp.log(1 + t**2) / 2          # antiderivative: F(x) = G(x) - G(0)\n",
    "eps = 1e-3\n",
    "x = onp.arange(0., 2., eps)\n",
    "F = onp.cumsum(eps * f(x))                    # area-so-far via Riemann sum\n",
    "dFdx = (F[1:] - F[:-1]) / eps                # telescopes to f(x[1:]) exactly\n",
    "print('max |dF/dx - f|      :', float(onp.abs(dFdx - f(x[1:])).max()))\n",
    "print('Riemann F(2) vs G(2)-G(0):', float(F[-1]), float(G(onp.array(2.))))"
   ],
   "execution_count": 3,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "max |dF/dx - f|      : 5.551115123125783e-14\nRiemann F(2) vs G(2)-G(0): 0.8045188628837103 0.8047189562170501\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-4 fw=all -->\n",
    "\n",
    "The first comparison is exact in exact arithmetic: differencing the cumulative\n",
    "sum *telescopes*, $F(x+\\epsilon)-F(x)$ collapsing back to the single term\n",
    "$\\epsilon f(x+\\epsilon)$ that the sum just added. The printed error, a few\n",
    "parts in $10^{5}$, is therefore not discretization error at all but pure\n",
    "float32 roundoff in the accumulated $F$, amplified by the division by\n",
    "$\\epsilon$. The Riemann total $F(2)$ agrees with the antiderivative value\n",
    "$G(2)$ to about $2\\times 10^{-4}$, the expected first-order bias of the\n",
    "left-endpoint rule.\n",
    "\n",
    "### Improper Integrals\n",
    "\n",
    "Many densities live on an unbounded domain (the Gaussian, the exponential, any\n",
    "heavy tail), so we will need to integrate \"all the way to infinity,\" yet\n",
    "the equation only defined integrals over a finite $[a,b]$. We fill the\n",
    "gap the same way we defined the integral itself: as a limit. The **improper\n",
    "integral** is\n",
    "\n",
    "$$\n",
    "\\int_a^\\infty f(x)\\,dx = \\lim_{b\\to\\infty}\\int_a^b f(x)\\,dx ,\n",
    "$$\n",
    "\n",
    "with $\\int_{-\\infty}^\\infty$ defined by splitting at any point and requiring\n",
    "*both* one-sided limits to exist on their own (a symmetric combination of the\n",
    "two is not enough). The limit may be finite, in which case the integral\n",
    "*converges*;\n",
    "otherwise, whether the partial integrals grow without bound or merely fail to\n",
    "settle, it *diverges*. Whether the tail is \"thin enough\" is a\n",
    "genuine question. The test case is the power law $x^{-p}$, for which the\n",
    "antiderivative gives\n",
    "\n",
    "$$\n",
    "\\int_1^\\infty x^{-p}\\,dx = \\lim_{b\\to\\infty}\\frac{b^{1-p}-1}{1-p}\n",
    "= \\begin{cases} \\dfrac{1}{p-1}, & p>1\\ \\text{(converges)},\\\\[1ex] \\infty, & p\\le 1\\ \\text{(diverges)}.\\end{cases}\n",
    "$$\n",
    "\n",
    "The antiderivative $\\tfrac{b^{1-p}-1}{1-p}$ is the $0/0$ indeterminate form at\n",
    "exactly $p=1$, where the integrand $x^{-1}$ has antiderivative $\\log x$ instead, so\n",
    "$\\int_1^\\infty x^{-1}\\,dx=\\lim_{b\\to\\infty}\\log b=\\infty$; the case split above\n",
    "already records this divergent value. So $\\int_1^\\infty x^{-2}\\,dx = 1$ while\n",
    "$\\int_1^\\infty x^{-1}\\,dx=\\infty$: the boundary between convergence and divergence\n",
    "sits exactly at $p=1$. This single threshold is what decides whether a\n",
    "*heavy-tailed* density, one whose tail decays like a power law rather than\n",
    "exponentially, even has a finite\n",
    "normalizer or mean, a recurring concern once we reach probability. The cell\n",
    "watches a convergent improper integral, $\\int_0^\\infty e^{-x}\\,dx = 1$, through\n",
    "two lenses at once: the *exact* partial integrals\n",
    "$\\int_0^b e^{-x}\\,dx = 1-e^{-b}$, read off the antiderivative, and a\n",
    "left-Riemann approximation of each with $\\epsilon = 10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-improper"
   },
   "source": [
    "import numpy as onp\n",
    "for b in [1., 2., 5., 10., 20.]:\n",
    "    x = onp.arange(0., b, 1e-3)\n",
    "    left = float(onp.sum(1e-3 * onp.exp(-x)))\n",
    "    exact = 1 - float(onp.exp(onp.array(-b)))\n",
    "    print(f'b={b:<4}  left-Riemann sum={left:.6f}  exact 1-e^-b={exact:.6f}')"
   ],
   "execution_count": 4,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "b=1.0   left-Riemann sum=0.632437  exact 1-e^-b=0.632121\nb=2.0   left-Riemann sum=0.865097  exact 1-e^-b=0.864665\nb=5.0   left-Riemann sum=0.993759  exact 1-e^-b=0.993262\nb=10.0  left-Riemann sum=1.000455  exact 1-e^-b=0.999955\nb=20.0  left-Riemann sum=1.000500  exact 1-e^-b=1.000000\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-5 fw=all -->\n",
    "\n",
    "Read the two columns against each other. The exact partials approach $1$\n",
    "quickly: the\n",
    "*truncation error* of stopping at $b$ is $e^{-b}$, already invisible at six\n",
    "decimals by $b=20$. The left-Riemann column tells a subtler story: it\n",
    "overshoots and settles near $1.0005$, not $1$. That residue is *discretization\n",
    "error*. The left rule samples each slice at its left edge, where the decreasing\n",
    "integrand $e^{-x}$ is largest, so every rectangle overcounts slightly; the\n",
    "overcounts add up to a bias of\n",
    "$\\approx \\tfrac{\\epsilon}{2}\\bigl(f(0)-f(b)\\bigr) \\to \\tfrac{\\epsilon}{2}\n",
    "= 0.0005$. The two error sources are independent knobs: growing $b$ removes\n",
    "the first, shrinking $\\epsilon$ the second, and no amount of turning one\n",
    "compensates for neglecting the other.\n",
    "\n",
    "### Integration by Parts\n",
    "\n",
    "The fundamental theorem turns every differentiation rule into an integration\n",
    "rule. Converted, the product rule becomes integration by parts.\n",
    "\n",
    "**Proposition (Integration by parts).** *For continuously differentiable $u$\n",
    "and $v$ on $[a,b]$,*\n",
    "\n",
    "$$\n",
    "\\int_a^b u(x)\\,v'(x)\\,dx\n",
    "= \\Bigl[\\,u(x)\\,v(x)\\,\\Bigr]_a^b - \\int_a^b u'(x)\\,v(x)\\,dx ,\n",
    "$$\n",
    "\n",
    "*where $[\\,h(x)\\,]_a^b$ abbreviates $h(b)-h(a)$.*\n",
    "\n",
    "**Proof.** By the product rule, $uv$ is an antiderivative of $u'v + uv'$, so\n",
    "the equation gives\n",
    "$\\int_a^b (u'v + uv')\\,dx = [\\,u(x)v(x)\\,]_a^b$; subtract $\\int_a^b u'v\\,dx$\n",
    "from both sides. $\\blacksquare$\n",
    "\n",
    "The formula is a *trade*, not an evaluation: it moves the derivative from one\n",
    "factor onto the other, which helps whenever the new integrand\n",
    "$u'v$ is simpler than the old one $uv'$. The model example takes $u(x)=x$ and\n",
    "$v'(x)=e^{-x}$, so that $u'=1$ and $v=-e^{-x}$: applying\n",
    "the equation on $[0,b]$ and letting $b\\to\\infty$,\n",
    "\n",
    "$$\n",
    "\\int_0^\\infty x\\,e^{-x}\\,dx\n",
    "= \\Bigl[\\,-x\\,e^{-x}\\,\\Bigr]_0^\\infty + \\int_0^\\infty e^{-x}\\,dx\n",
    "= 0 + 1 = 1 .\n",
    "$$\n",
    "\n",
    "The boundary term vanishes because $x e^{-x}\\to 0$ at both ends, and what\n",
    "survives is exactly the improper integral computed above. This particular\n",
    "value is the mean of the exponential distribution\n",
    "(that section), and it shows the form in which the trick\n",
    "returns throughout probability: expectations $\\int g(x)\\,p(x)\\,dx$ are\n",
    "integrals of products, and shifting a derivative from a factor we cannot handle\n",
    "onto one we can is often the only move available. The identity returns in\n",
    "Hyvärinen's score matching [@Hyvarinen.2005]\n",
    "(that section), where a single integration\n",
    "by parts converts an objective involving the score of the data\n",
    "distribution, the gradient $\\nabla_{\\mathbf{x}} \\log p(\\mathbf{x})$ of the\n",
    "log-density with respect to the *data point*, which no one can evaluate,\n",
    "into one that can be estimated from samples.\n",
    "\n",
    "### A Note on Signed Area\n",
    "\n",
    "The evaluation rule the equation produces negative numbers,\n",
    "which can be jarring if \"area\" is supposed to be positive. The resolution is that\n",
    "integrals measure *signed* area, governed by two independent sign rules. First,\n",
    "where $f<0$ the area counts as negative: $\\int_0^1 (-1)\\,dx = -1$. Second,\n",
    "integrating right-to-left negates the result: $\\int_b^a f = -\\int_a^b f$. Each\n",
    "flip, reflecting across the $x$-axis or reversing the limits, introduces one\n",
    "minus sign; the figure shows both rules on one curve. Apply\n",
    "them in turn to a single example: $\\int_0^{-1} 1\\,dx = -1$ (one\n",
    "flip, from reversed limits), while $\\int_0^{-1}(-1)\\,dx = 1$ (two flips, which\n",
    "cancel). This is the same signed-area convention the\n",
    "determinant followed for transformed regions in\n",
    "that section, and it is exactly what makes the\n",
    "change-of-variables formula below come out right.\n",
    "\n",
    "![Signed area. Lobes above the axis count positive and the lobe below counts negative, so $\\int_a^b f$ totals the shaded regions with their signs; reversing the limits of integration traverses the same area right-to-left and negates the total.](../img/mdl-cal-signed-area.svg)\n",
    "\n",
    "## Change of Variables\n",
    "\n",
    "Every differentiation rule has an integration counterpart obtained by running the\n",
    "fundamental theorem backwards: the sum rule gives linearity, the product rule\n",
    "gave us integration by parts (that section),\n",
    "and the chain rule, most useful of all for us, gives the\n",
    "*change-of-variables* formula: the substitution rule that evaluates many\n",
    "integrals the derivative table alone cannot. The next section carries it into\n",
    "higher dimensions.\n",
    "\n",
    "Suppose we reparametrize the variable through a function $u$. The chain rule says\n",
    "$\\frac{d}{dx}F(u(x)) = F'(u(x))\\,u'(x)$; feeding this through the fundamental\n",
    "theorem with $F'=f$ and integrating both sides yields the\n",
    "**change-of-variables formula**\n",
    "\n",
    "$$\n",
    "\\int_{u(a)}^{u(b)} f(y)\\,dy = \\int_a^b f(u(x))\\,\\frac{du}{dx}\\,dx .\n",
    "$$\n",
    "\n",
    "The picture explains the mysterious factor $\\frac{du}{dx}$. Look at one thin\n",
    "rectangle. On the left side of the equation, the sliver from $x$ to\n",
    "$x+\\epsilon$ has area $\\approx \\epsilon\\,f(u(x))$. On the right, the corresponding\n",
    "sliver runs from $u(x)$ to $u(x+\\epsilon)\\approx u(x)+\\epsilon\\,u'(x)$, so its\n",
    "width is *stretched* by the factor $u'(x)$ and its area is\n",
    "$\\approx \\epsilon\\,u'(x)\\,f(u(x))$. To make the two slivers agree we must multiply\n",
    "by exactly the local stretch $\\frac{du}{dx}$, as the figure\n",
    "shows.\n",
    "\n",
    "![A thin rectangle under the substitution $y=u(x)$. The sliver of width $\\epsilon$ at $x$ maps to a sliver of width $\\epsilon\\,u'(x)$; matching the two areas forces the factor $du/dx$ in the change-of-variables formula.](../img/mdl-cal-rect-trans.svg)\n",
    "\n",
    "With the right $u$ this collapses hard integrals to trivial ones. Taking $f=1$\n",
    "and $u(x)=e^{-x^2}$ (so $u'(x)=-2x\\,e^{-x^2}$),\n",
    "\n",
    "$$\n",
    "e^{-1}-1 = \\int_{e^{0}}^{e^{-1}} 1\\,dy = -2\\int_0^{1} x\\,e^{-x^2}\\,dx\n",
    "\\quad\\Longrightarrow\\quad\n",
    "\\int_0^{1} x\\,e^{-x^2}\\,dx = \\frac{1-e^{-1}}{2},\n",
    "$$\n",
    "\n",
    "a closed form that no amount of rectangle-summing would have revealed. (Keep the\n",
    "two variables apart: $y$ is the transformed variable on the left of\n",
    "the equation, $x$ the original one on the right; the result is a\n",
    "statement about the $x$-integral.)\n",
    "\n",
    "## Multiple Integrals\n",
    "\n",
    "In one dimension the integral was a limit of rectangle sums, and nothing in that\n",
    "recipe was attached to the line. This section runs the same recipe over regions\n",
    "of $\\mathbb{R}^n$, establishes Fubini's theorem, which reduces such\n",
    "integrals to iterated one-dimensional ones, and lifts the change-of-variables\n",
    "rule to $n$ dimensions.\n",
    "\n",
    "### Double Integrals\n",
    "\n",
    "In higher dimensions we integrate over regions. For $f(x,y)$ on a rectangle\n",
    "$U=[a,b]\\times[c,d]$, the integral is the (signed) volume trapped between the\n",
    "surface $z=f(x,y)$ and the base plane, written $\\int_U f(x,y)\\,dx\\,dy$, and\n",
    "the rectangle recipe survives intact: tile $U$ with $\\epsilon\\times\\epsilon$\n",
    "squares, stand a box of height $f$ on each, and total the volumes.\n",
    "the figure shows our running example, the two-dimensional\n",
    "bell $e^{-x^2-y^2}$ over the box $[-2,2]^2$; the cell totals its boxes.\n",
    "\n",
    "![The bell surface $z = e^{-x^2-y^2}$ over the box $U = \\lbrack -2,2\\rbrack^2$. The double integral $\\int_U f\\,dx\\,dy$ is the volume between the surface and the base plane, approximated by tiling the base with $\\epsilon\\times\\epsilon$ squares and adding the volumes of the boxes standing on them.](../img/mdl-cal-bell-surface.svg)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-box-volume"
   },
   "source": [
    "import math\n",
    "import numpy as onp\n",
    "eps = 1e-2\n",
    "x, y = onp.meshgrid(onp.arange(-2., 2., eps), onp.arange(-2., 2., eps),\n",
    "                   indexing='ij')\n",
    "vol = float(onp.sum(eps ** 2 * onp.exp(-x**2 - y**2)))\n",
    "exact = math.pi * math.erf(2.) ** 2\n",
    "print(f'Riemann volume over the box [-2, 2]^2 : {vol:.5f}')\n",
    "print(f'exact volume  pi * erf(2)^2           : {exact:.5f}')"
   ],
   "execution_count": 5,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Riemann volume over the box [-2, 2]^2 : 3.11227\nexact volume  pi * erf(2)^2           : 3.11227\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-6 fw=all -->\n",
    "\n",
    "The exact value involves the **error function**\n",
    "$\\operatorname{erf}(z)=\\tfrac{2}{\\sqrt{\\pi}}\\int_0^z e^{-t^2}\\,dt$, the\n",
    "Gaussian's antiderivative, suitably normalized. It has no elementary formula,\n",
    "but the library provides it as a built-in, which is how the cell computes the\n",
    "exact volume $\\pi\\operatorname{erf}(2)^2$ (an expression derived at the end of\n",
    "this section, where the Gaussian integral is evaluated).\n",
    "\n",
    "### Fubini's Theorem\n",
    "\n",
    "How do we get from grid sums back to the fundamental theorem? Discretize the\n",
    "region into $\\epsilon\\times\n",
    "\\epsilon$ squares indexed by integers $i,j$; the integral is approximately\n",
    "$\\sum_{i,j}\\epsilon^2 f(\\epsilon i, \\epsilon j)$. This is just a sum of numbers\n",
    "on a grid, and **a finite sum can be totalled in any order**. Summing\n",
    "columns-first (add up each column $i$, then add the column totals) we get\n",
    "\n",
    "$$\n",
    "\\sum_j \\epsilon\\!\\left(\\sum_i \\epsilon\\, f(\\epsilon i, \\epsilon j)\\right),\n",
    "$$\n",
    "\n",
    "illustrated in the figure. The inner sum is the Riemann\n",
    "discretization of $\\int_a^b f(x,\\epsilon j)\\,dx$, and the outer sum then\n",
    "integrates that over $y$.\n",
    "\n",
    "![Summing a grid of $\\epsilon\\times\\epsilon$ cells columns-first (1), then adding the column totals together (2). Reordering a finite sum changes nothing; this is the discrete fact behind iterated integration.](../img/mdl-cal-sum-order.svg)\n",
    "\n",
    "For the finite grid sums the reordering is trivial, and that triviality carries\n",
    "none of the theorem's content: the integral is the *limit* of the grid sums,\n",
    "and the substance is that the interchange survives the two limiting processes.\n",
    "That step we take on faith. **Fubini's theorem** [@Folland.1999] requires\n",
    "$f$ to be **absolutely integrable**, meaning\n",
    "\n",
    "$$\n",
    "\\int_U |f(x,y)|\\,dx\\,dy < \\infty ,\n",
    "$$\n",
    "\n",
    "and concludes that the double integral equals either iterated integral,\n",
    "\n",
    "$$\n",
    "\\int_U f(x,y)\\,dx\\,dy\n",
    "= \\int_c^{d}\\!\\left(\\int_a^{b} f(x,y)\\,dx\\right)dy\n",
    "= \\int_a^{b}\\!\\left(\\int_c^{d} f(x,y)\\,dy\\right)dx .\n",
    "$$\n",
    "\n",
    "The hypothesis matters. The standard\n",
    "counterexample is $f(x,y)=(x^2-y^2)/(x^2+y^2)^2$ on $[0,1]^2$: integrating $x$\n",
    "first gives $-\\tfrac{\\pi}{4}$, integrating $y$ first gives $+\\tfrac{\\pi}{4}$, and\n",
    "the orders disagree precisely because the singularity at the origin makes\n",
    "$\\int_U |f|\\,dx\\,dy$ infinite. For the continuous, absolutely integrable\n",
    "functions of machine learning\n",
    "there is no such trouble, and we freely swap orders. When the region $U$ is more\n",
    "complicated than a rectangle we write the integral compactly as $\\int_U\n",
    "f(\\mathbf{x})\\,d\\mathbf{x}$; it is defined by *zero-extension*: extend $f$ by\n",
    "$0$ outside $U$ and integrate the extended function over any rectangle\n",
    "containing $U$.\n",
    "\n",
    "### Change of Variables in Many Dimensions\n",
    "\n",
    "The one-dimensional stretch factor $\\frac{du}{dx}$ has a multidimensional\n",
    "successor, and it is the key that unlocks the Gaussian integral below. Let $U\\subseteq\\mathbb{R}^n$\n",
    "be open and let $\\boldsymbol{\\phi}:U\\to\\mathbb{R}^n$ be a **$C^1$-diffeomorphism**\n",
    "onto its image; that is, $\\boldsymbol{\\phi}$ is continuously differentiable,\n",
    "injective, and has $\\det D\\boldsymbol{\\phi}(\\mathbf{x})\\neq 0$ everywhere on $U$.\n",
    "(Injectivity stops the map folding space onto itself; the nonvanishing Jacobian\n",
    "keeps it from collapsing a region to lower dimension, and together they make\n",
    "$\\boldsymbol{\\phi}$ invertible on $\\boldsymbol{\\phi}(U)$ with a $C^1$ inverse.)\n",
    "Then, for absolutely integrable $f$ (again,\n",
    "$\\int |f|\\,d\\mathbf{x} < \\infty$), a theorem we also take on faith\n",
    "[@Folland.1999] states that\n",
    "\n",
    "$$\n",
    "\\int_{\\boldsymbol{\\phi}(U)} f(\\mathbf{x})\\,d\\mathbf{x}\n",
    "= \\int_{U} f(\\boldsymbol{\\phi}(\\mathbf{x}))\\,\n",
    "  \\bigl|\\det D\\boldsymbol{\\phi}(\\mathbf{x})\\bigr|\\,d\\mathbf{x},\n",
    "$$\n",
    "\n",
    "where $D\\boldsymbol{\\phi}$ is the **Jacobian** of $\\boldsymbol{\\phi}$, the\n",
    "matrix of partial derivatives introduced in that section.\n",
    "\n",
    "Compare the equation with the one-dimensional rule\n",
    "the equation: the scalar stretch $\\frac{du}{dx}$ has been replaced\n",
    "by $|\\det D\\boldsymbol{\\phi}|$, and the determinant of a linear\n",
    "map is exactly the factor by which it scales area or volume, the signed-area\n",
    "reading of the determinant from\n",
    "that section. So a tiny cube of volume\n",
    "$d\\mathbf{x}$ is sent to a tiny parallelepiped of volume\n",
    "$|\\det D\\boldsymbol{\\phi}|\\,d\\mathbf{x}$, and the factor in the integrand\n",
    "compensates for exactly that distortion. the figure puts\n",
    "the one- and two-dimensional pictures side by side.\n",
    "\n",
    "![Left: the 1-D substitution scales a segment by $du/dx$. Right: in 2-D a linear $\\boldsymbol{\\phi}$ sends the unit square to a parallelogram whose area is the local volume-scaling factor $|\\det D\\boldsymbol{\\phi}|$, the multidimensional successor of the stretch factor.](../img/mdl-cal-cov-jacobian.svg)\n",
    "\n",
    "Equation the equation is a statement about *area and volume*: it\n",
    "says how the integral of a fixed function transports under a reparametrization of\n",
    "the domain. Applied instead to a probability *density*, a function that must keep\n",
    "integrating to $1$, the very same Jacobian factor becomes the\n",
    "change-of-variables-for-densities rule of that section, whose\n",
    "$-\\log|\\det D\\boldsymbol{\\phi}|$ correction is the exact mechanism behind\n",
    "**normalizing flows** (that section). We state the\n",
    "area theorem here and defer that density specialization to the probability chapter.\n",
    "\n",
    "### The Gaussian Integral\n",
    "\n",
    "The machinery just assembled evaluates the **Gaussian\n",
    "integral**, which we will meet again as the normalizer of the normal\n",
    "distribution. It is the *one-dimensional* integral\n",
    "\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} e^{-x^2}\\,dx\n",
    "$$\n",
    "\n",
    "that resists direct attack: $e^{-x^2}$ has no elementary antiderivative (its\n",
    "antiderivative is the error function from the bell-volume cell, *defined* as\n",
    "that integral rather than expressed in older functions), so the fundamental\n",
    "theorem does not apply directly. The trick is to go *up* a dimension.\n",
    "\n",
    "One preliminary. The plan involves integrals over all of $\\mathbb{R}^2$, and we\n",
    "have defined improper integrals only on the line the equation. For\n",
    "a *nonnegative* integrand the same limiting idea works in any dimension, by\n",
    "exhaustion: define $\\iint_{\\mathbb{R}^2} f\\,dx\\,dy$ as the limit of the\n",
    "integrals over the growing boxes $[-b,b]^2$ as $b\\to\\infty$. Because $f\\ge 0$,\n",
    "these integrals can only increase as the box grows, so the limit exists (it may\n",
    "be $+\\infty$), and exhausting the plane by discs or any other expanding family\n",
    "of regions gives the same value, since every disc is contained in some box and\n",
    "every box in some disc. Everything we integrate below is nonnegative, so this\n",
    "definition is in force throughout.\n",
    "\n",
    "Now square the integral and rename the dummy variable in the second copy:\n",
    "\n",
    "$$\n",
    "\\left(\\int_{-\\infty}^{\\infty} e^{-x^2}\\,dx\\right)^{\\!2}\n",
    "= \\int_{-\\infty}^{\\infty} e^{-x^2}\\,dx \\int_{-\\infty}^{\\infty} e^{-y^2}\\,dy\n",
    "= \\int_{-\\infty}^{\\infty}\\!\\!\\int_{-\\infty}^{\\infty} e^{-x^2-y^2}\\,dx\\,dy .\n",
    "$$\n",
    "\n",
    "The second equality says that the iterated integral equals the double integral\n",
    "over the whole plane: this is Fubini's theorem again, now on $\\mathbb{R}^2$\n",
    "rather than on a rectangle, and we grant this whole-plane form on the same\n",
    "authority [@Folland.1999].\n",
    "\n",
    "The two-dimensional integral has something its one-dimensional parent\n",
    "lacked, rotational symmetry, and polar coordinates\n",
    "$\\boldsymbol{\\phi}(r,\\theta)=(r\\cos\\theta, r\\sin\\theta)$ are built to exploit\n",
    "it. Take the open set $U=(0,\\infty)\\times(0,2\\pi)$. On $U$ the map\n",
    "$\\boldsymbol{\\phi}$ is injective, and its Jacobian determinant\n",
    "\n",
    "$$\n",
    "\\bigl|\\det D\\boldsymbol{\\phi}\\bigr|\n",
    "= \\left|\\det\\begin{bmatrix}\\cos\\theta & -r\\sin\\theta\\\\ \\sin\\theta & r\\cos\\theta\\end{bmatrix}\\right|\n",
    "= r(\\cos^2\\theta + \\sin^2\\theta) = r\n",
    "$$\n",
    "\n",
    "is nonzero there, so $\\boldsymbol{\\phi}$ is a $C^1$-diffeomorphism of $U$ onto\n",
    "its image. That image is not quite the whole plane: $\\boldsymbol{\\phi}(U)$ is\n",
    "$\\mathbb{R}^2$ minus the ray $\\{(x, 0) : x \\ge 0\\}$, the points the excluded\n",
    "angle $\\theta = 0$ would have produced. A ray, however, has *measure zero*: it\n",
    "can be covered by rectangles of arbitrarily small total area (a notion\n",
    "that section returns to), and we grant, once more from\n",
    "the same source [@Folland.1999], that a set of measure zero contributes\n",
    "nothing to any integral. Integrating over $\\boldsymbol{\\phi}(U)$ is therefore\n",
    "the same as integrating over the plane, and the change of variables\n",
    "the equation on $U$ hands the integrand the factor $r$ that\n",
    "makes it elementary:\n",
    "\n",
    "$$\n",
    "\\iint_{\\mathbb{R}^2} e^{-x^2-y^2}\\,dx\\,dy\n",
    "= \\int_{\\boldsymbol{\\phi}(U)} e^{-x^2-y^2}\\,dx\\,dy\n",
    "= \\int_0^\\infty\\!\\!\\int_0^{2\\pi} r\\,e^{-r^2}\\,d\\theta\\,dr\n",
    "= 2\\pi\\int_0^\\infty r\\,e^{-r^2}\\,dr = \\pi ,\n",
    "$$\n",
    "\n",
    "the last step because $-\\tfrac12 e^{-r^2}$ is an antiderivative of\n",
    "$r\\,e^{-r^2}$, so the improper integral is $\\tfrac12$. (Strictly, one applies\n",
    "the equation on the bounded sets\n",
    "$U_b=(0,b)\\times(0,2\\pi)$, where the continuous bounded integrand is certainly\n",
    "absolutely integrable, and lets $b\\to\\infty$; both sides increase to the values\n",
    "shown.) Taking the square root of the equation (the positive\n",
    "one, since the integrand is positive) gives\n",
    "\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} e^{-x^2}\\,dx = \\sqrt{\\pi} .\n",
    "$$\n",
    "\n",
    "The same factorization also explains the bell-volume cell's exact value:\n",
    "over the box, $\\bigl(\\int_{-2}^{2} e^{-x^2}\\,dx\\bigr)^2\n",
    "= \\bigl(\\sqrt{\\pi}\\,\\operatorname{erf}(2)\\bigr)^2\n",
    "= \\pi\\operatorname{erf}(2)^2 \\approx 3.11227$, shy of the full-plane $\\pi$ by\n",
    "exactly the truncated tails. The cell verifies the two-dimensional integral\n",
    "numerically by summing the integrand over a fine grid on $[-6,6]^2$, a box\n",
    "wide enough that the Gaussian's tails beyond it contribute nothing to six\n",
    "decimals."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-gaussian"
   },
   "source": [
    "import numpy as onp\n",
    "g = onp.arange(-6., 6., 0.01)\n",
    "x, y = onp.meshgrid(g, g, indexing='ij')\n",
    "val = float(onp.sum(0.01 ** 2 * onp.exp(-x**2 - y**2)))\n",
    "print('grid integral of e^(-x^2-y^2):', round(val, 6))\n",
    "print('exact value pi               :', round(float(onp.pi), 6))"
   ],
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "grid integral of e^(-x^2-y^2): 3.141593\nexact value pi               : 3.141593\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-7 fw=all -->\n",
    "\n",
    "## Integration Meets Probability\n",
    "\n",
    "This is why a deep-learning reader needs integration at all. A continuous\n",
    "probability **density** is nothing more than a non-negative function that is\n",
    "*normalized*, meaning its total integral is $1$:\n",
    "\n",
    "$$\n",
    "p(x)\\ge 0, \\qquad \\int_{\\mathcal X} p(x)\\,dx = 1 .\n",
    "$$\n",
    "\n",
    "The Gaussian integral above is precisely what supplies the normalizer: we\n",
    "showed $\\int_{-\\infty}^\\infty e^{-x^2}\\,dx=\\sqrt\\pi$, so the function\n",
    "$p(x)=\\tfrac{1}{\\sqrt\\pi}e^{-x^2}$ integrates to $1$ and is a bona fide density.\n",
    "An **expectation** is then an integral average, weighting each value by its\n",
    "density,\n",
    "\n",
    "$$\n",
    "\\mathbb{E}[X] = \\int_{\\mathcal X} x\\,p(x)\\,dx,\n",
    "\\qquad\n",
    "\\mathbb{E}[g(X)] = \\int_{\\mathcal X} g(x)\\,p(x)\\,dx .\n",
    "$$\n",
    "\n",
    "In several variables a joint density satisfies $\\iint p(x,y)\\,dx\\,dy = 1$, and\n",
    "*marginals* are obtained by integrating a variable out, $p(x)=\\int p(x,y)\\,dy$;\n",
    "Fubini the equation guarantees the order of integration does not matter\n",
    "(densities are absolutely integrable, so the caveat above never bites). The cell\n",
    "verifies that $\\tfrac{1}{\\sqrt\\pi}e^{-x^2}$ is normalized and computes its mean,\n",
    "which symmetry sends to $0$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-density"
   },
   "source": [
    "import numpy as onp\n",
    "x = onp.arange(-8., 8., 1e-3)\n",
    "p = onp.exp(-x**2) / onp.sqrt(onp.array(onp.pi))\n",
    "print('total mass  integral p dx :', round(float(onp.sum(1e-3 * p)), 6))\n",
    "print('mean        integral x p dx:', round(float(onp.sum(1e-3 * x * p)), 6))"
   ],
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "total mass  integral p dx : 1.0\nmean        integral x p dx: 0.0\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-8 fw=all -->\n",
    "\n",
    "When an expectation has no closed form, the rule rather than the exception in\n",
    "deep learning, we estimate it by **Monte Carlo**\n",
    "[@Metropolis.Ulam.1949]: draw samples $x_i\\sim p$ and\n",
    "average,\n",
    "\n",
    "$$\n",
    "\\mathbb{E}[g(X)] \\approx \\frac{1}{n}\\sum_{i=1}^n g(x_i),\n",
    "$$\n",
    "\n",
    "a stochastic counterpart of the Riemann sum that ignores the geometry of the\n",
    "domain. The **law of large numbers** (that section) is the\n",
    "guarantee that it works: the sample\n",
    "average converges to the true expectation. The error scale comes from a\n",
    "variance computation: the estimator's standard deviation is\n",
    "$\\mathrm{sd}(g(X))/\\sqrt{n}$, an order-$1/\\sqrt{n}$ error whose limiting\n",
    "distribution the **central limit theorem** (that section)\n",
    "supplies. The rate depends only on the sample size $n$, not\n",
    "on the dimension. (The fine print: the *constant* in front of\n",
    "$1/\\sqrt{n}$ is the standard deviation of $g(X)$, and *that* can grow with the\n",
    "dimension; what is dimension-free is the exponent, and it is the exponent that\n",
    "grids lose.) That dimension-free rate is decisive. A grid laid down to\n",
    "resolution $\\epsilon$ in $d$ dimensions costs $N=\\epsilon^{-d}$ evaluations.\n",
    "The left rule we have used throughout has error $\\propto\\epsilon$, so it\n",
    "achieves only $N^{-1/d}$; even a *second-order* rule (the midpoint rule, which\n",
    "samples each slice at its center, or the trapezoid rule, which averages the two\n",
    "edge heights) has error $\\propto\\epsilon^2$ and manages just $N^{-2/d}$.\n",
    "Either way the\n",
    "exponent is divided by the dimension, the *curse of\n",
    "dimensionality*, while Monte Carlo keeps its $N^{-1/2}$ in every dimension.\n",
    "This is why sampling, rather than\n",
    "grid quadrature, is how expectations are computed at scale. The sample-mean\n",
    "estimator and its $\\sigma^2/n$ variance return with random variables\n",
    "(that section), and the law of large numbers behind the\n",
    "guarantee is stated in that section.\n",
    "The cell makes this concrete: it estimates $\\int_0^1 e^{-x^2}\\,dx$ by averaging\n",
    "$e^{-x_i^2}$ over uniform samples, prints the running estimate against a\n",
    "midpoint quadrature fine enough to be exact at the printed precision, and then\n",
    "plots the measured error of the running estimate against the sample size on\n",
    "log-log axes, next to the predicted guide line $0.2\\,N^{-1/2}$ (the constant\n",
    "$0.2$ is $\\mathrm{sd}(e^{-X^2})$ for uniform $X$)."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "integral-monte-carlo"
   },
   "source": [
    "np.random.seed(1)\n",
    "# E[e^{-X^2}], X ~ U[0, 1], by Monte Carlo vs a fine midpoint quadrature\n",
    "# (midpoint is second-order: at eps=1e-4 its bias ~ 3e-10, below the print).\n",
    "quad = float(np.sum(1e-4 * np.exp(-(np.arange(0., 1., 1e-4) + 5e-5) ** 2)))\n",
    "draws = np.exp(-np.random.rand(100000) ** 2)\n",
    "running = np.cumsum(draws) / np.arange(1., 100001.)   # running MC estimate\n",
    "for n in [10, 100, 1000, 10000, 100000]:\n",
    "    mc = float(running[n - 1])\n",
    "    print(f'n={n:<7} Monte Carlo={mc:.5f}  error={abs(mc - quad):.1e}')\n",
    "print(f'midpoint quadrature  ={quad:.5f}')\n",
    "\n",
    "# Measured error of the running estimate vs the predicted sd(g)/sqrt(N)\n",
    "# line, sampled at log-spaced N so the log-log plot stays light.\n",
    "Ns = np.logspace(1., 5., 201)\n",
    "err = np.abs(running[Ns.astype('int32') - 1] - quad)\n",
    "d2l.plot(Ns, [err, 0.2 * Ns ** -0.5], 'samples N', 'integration error',\n",
    "         xscale='log', yscale='log',\n",
    "         legend=['measured MC error', '$0.2\\\\,N^{-1/2}$ guide'])"
   ],
   "execution_count": 8,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "n=10      Monte Carlo=0.87741  error=1.3e-01\nn=100     Monte Carlo=0.75454  error=7.7e-03\nn=1000    Monte Carlo=0.74619  error=6.3e-04\nn=10000   Monte Carlo=0.74833  error=1.5e-03\nn=100000  Monte Carlo=0.74719  error=3.7e-04\nmidpoint quadrature  =0.74682\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"253.878125pt\" height=\"183.35625pt\" viewBox=\"0 0 253.878125 183.35625\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n <metadata>\n  <rdf:RDF xmlns:dc=\"http://purl.org/dc/elements/1.1/\" xmlns:cc=\"http://creativecommons.org/ns#\" xmlns:rdf=\"http://www.w3.org/1999/02/22-rdf-syntax-ns#\">\n   <cc:Work>\n    <dc:type rdf:resource=\"http://purl.org/dc/dcmitype/StillImage\"/>\n    <dc:date>2026-07-19T13:53:18.365664</dc:date>\n    <dc:format>image/svg+xml</dc:format>\n    <dc:creator>\n     <cc:Agent>\n      <dc:title>Matplotlib v3.10.8, https://matplotlib.org/</dc:title>\n     </cc:Agent>\n    </dc:creator>\n   </cc:Work>\n  </rdf:RDF>\n </metadata>\n <defs>\n  <style type=\"text/css\">*{stroke-linejoin: round; stroke-linecap: butt}</style>\n </defs>\n <g id=\"figure_1\">\n  <g id=\"patch_1\">\n   <path d=\"M 0 183.35625 \nL 253.878125 183.35625 \nL 253.878125 0 \nL 0 0 \nz\n\" style=\"fill: #ffffff\"/>\n  </g>\n  <g id=\"axes_1\">\n   <g id=\"patch_2\">\n    <path d=\"M 51.378125 145.8 \nL 246.678125 145.8 \nL 246.678125 7.2 \nL 51.378125 7.2 \nz\n\" style=\"fill: #ffffff\"/>\n   </g>\n   <g id=\"matplotlib.axis_1\">\n    <g id=\"xtick_1\">\n     <g id=\"line2d_1\">\n      <path d=\"M 60.255398 145.8 \nL 60.255398 7.2 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_2\">\n      <defs>\n       <path id=\"maafd03ff68\" d=\"M 0 0 \nL 0 3.5 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#maafd03ff68\" x=\"60.255398\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_1\">\n      <!-- $\\mathdefault{10^{1}}$ -->\n      <g transform=\"translate(51.455398 160.398438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-31\" d=\"M 794 531 \nL 1825 531 \nL 1825 4091 \nL 703 3866 \nL 703 4441 \nL 1819 4666 \nL 2450 4666 \nL 2450 531 \nL 3481 531 \nL 3481 0 \nL 794 0 \nL 794 531 \nz\n\" transform=\"scale(0.015625)\"/>\n        <path id=\"DejaVuSans-30\" d=\"M 2034 4250 \nQ 1547 4250 1301 3770 \nQ 1056 3291 1056 2328 \nQ 1056 1369 1301 889 \nQ 1547 409 2034 409 \nQ 2525 409 2770 889 \nQ 3016 1369 3016 2328 \nQ 3016 3291 2770 3770 \nQ 2525 4250 2034 4250 \nz\nM 2034 4750 \nQ 2819 4750 3233 4129 \nQ 3647 3509 3647 2328 \nQ 3647 1150 3233 529 \nQ 2819 -91 2034 -91 \nQ 1250 -91 836 529 \nQ 422 1150 422 2328 \nQ 422 3509 836 4129 \nQ 1250 4750 2034 4750 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_2\">\n     <g id=\"line2d_3\">\n      <path d=\"M 104.641761 145.8 \nL 104.641761 7.2 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_4\">\n      <g>\n       <use xlink:href=\"#maafd03ff68\" x=\"104.641761\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_2\">\n      <!-- $\\mathdefault{10^{2}}$ -->\n      <g transform=\"translate(95.841761 160.398438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-32\" d=\"M 1228 531 \nL 3431 531 \nL 3431 0 \nL 469 0 \nL 469 531 \nQ 828 903 1448 1529 \nQ 2069 2156 2228 2338 \nQ 2531 2678 2651 2914 \nQ 2772 3150 2772 3378 \nQ 2772 3750 2511 3984 \nQ 2250 4219 1831 4219 \nQ 1534 4219 1204 4116 \nQ 875 4013 500 3803 \nL 500 4441 \nQ 881 4594 1212 4672 \nQ 1544 4750 1819 4750 \nQ 2544 4750 2975 4387 \nQ 3406 4025 3406 3419 \nQ 3406 3131 3298 2873 \nQ 3191 2616 2906 2266 \nQ 2828 2175 2409 1742 \nQ 1991 1309 1228 531 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_3\">\n     <g id=\"line2d_5\">\n      <path d=\"M 149.028125 145.8 \nL 149.028125 7.2 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_6\">\n      <g>\n       <use xlink:href=\"#maafd03ff68\" x=\"149.028125\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_3\">\n      <!-- $\\mathdefault{10^{3}}$ -->\n      <g transform=\"translate(140.228125 160.398438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-33\" d=\"M 2597 2516 \nQ 3050 2419 3304 2112 \nQ 3559 1806 3559 1356 \nQ 3559 666 3084 287 \nQ 2609 -91 1734 -91 \nQ 1441 -91 1130 -33 \nQ 819 25 488 141 \nL 488 750 \nQ 750 597 1062 519 \nQ 1375 441 1716 441 \nQ 2309 441 2620 675 \nQ 2931 909 2931 1356 \nQ 2931 1769 2642 2001 \nQ 2353 2234 1838 2234 \nL 1294 2234 \nL 1294 2753 \nL 1863 2753 \nQ 2328 2753 2575 2939 \nQ 2822 3125 2822 3475 \nQ 2822 3834 2567 4026 \nQ 2313 4219 1838 4219 \nQ 1578 4219 1281 4162 \nQ 984 4106 628 3988 \nL 628 4550 \nQ 988 4650 1302 4700 \nQ 1616 4750 1894 4750 \nQ 2613 4750 3031 4423 \nQ 3450 4097 3450 3541 \nQ 3450 3153 3228 2886 \nQ 3006 2619 2597 2516 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-33\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_4\">\n     <g id=\"line2d_7\">\n      <path d=\"M 193.414489 145.8 \nL 193.414489 7.2 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_8\">\n      <g>\n       <use xlink:href=\"#maafd03ff68\" x=\"193.414489\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_4\">\n      <!-- $\\mathdefault{10^{4}}$ -->\n      <g transform=\"translate(184.614489 160.398438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-34\" d=\"M 2419 4116 \nL 825 1625 \nL 2419 1625 \nL 2419 4116 \nz\nM 2253 4666 \nL 3047 4666 \nL 3047 1625 \nL 3713 1625 \nL 3713 1100 \nL 3047 1100 \nL 3047 0 \nL 2419 0 \nL 2419 1100 \nL 313 1100 \nL 313 1709 \nL 2253 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-34\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_5\">\n     <g id=\"line2d_9\">\n      <path d=\"M 237.800852 145.8 \nL 237.800852 7.2 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_10\">\n      <g>\n       <use xlink:href=\"#maafd03ff68\" x=\"237.800852\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_5\">\n      <!-- $\\mathdefault{10^{5}}$ -->\n      <g transform=\"translate(229.000852 160.398438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-35\" d=\"M 691 4666 \nL 3169 4666 \nL 3169 4134 \nL 1269 4134 \nL 1269 2991 \nQ 1406 3038 1543 3061 \nQ 1681 3084 1819 3084 \nQ 2600 3084 3056 2656 \nQ 3513 2228 3513 1497 \nQ 3513 744 3044 326 \nQ 2575 -91 1722 -91 \nQ 1428 -91 1123 -41 \nQ 819 9 494 109 \nL 494 744 \nQ 775 591 1075 516 \nQ 1375 441 1709 441 \nQ 2250 441 2565 725 \nQ 2881 1009 2881 1497 \nQ 2881 1984 2565 2268 \nQ 2250 2553 1709 2553 \nQ 1456 2553 1204 2497 \nQ 953 2441 691 2322 \nL 691 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-35\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_6\">\n     <g id=\"line2d_11\">\n      <defs>\n       <path id=\"m111e93ffd2\" d=\"M 0 0 \nL 0 2 \n\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"53.379863\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_7\">\n     <g id=\"line2d_12\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"55.953915\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_8\">\n     <g id=\"line2d_13\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"58.224389\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_9\">\n     <g id=\"line2d_14\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"73.617025\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_10\">\n     <g id=\"line2d_15\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"81.433075\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_11\">\n     <g id=\"line2d_16\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"86.978651\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_12\">\n     <g id=\"line2d_17\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"91.280135\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_13\">\n     <g id=\"line2d_18\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"94.794702\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_14\">\n     <g id=\"line2d_19\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"97.766227\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_15\">\n     <g id=\"line2d_20\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"100.340278\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_16\">\n     <g id=\"line2d_21\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"102.610753\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_17\">\n     <g id=\"line2d_22\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"118.003388\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_18\">\n     <g id=\"line2d_23\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"125.819439\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_19\">\n     <g id=\"line2d_24\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"131.365015\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_20\">\n     <g id=\"line2d_25\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"135.666498\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_21\">\n     <g id=\"line2d_26\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"139.181066\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_22\">\n     <g id=\"line2d_27\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"142.15259\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_23\">\n     <g id=\"line2d_28\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"144.726642\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_24\">\n     <g id=\"line2d_29\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"146.997116\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_25\">\n     <g id=\"line2d_30\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"162.389752\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_26\">\n     <g id=\"line2d_31\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"170.205803\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_27\">\n     <g id=\"line2d_32\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"175.751379\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_28\">\n     <g id=\"line2d_33\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"180.052862\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_29\">\n     <g id=\"line2d_34\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"183.567429\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_30\">\n     <g id=\"line2d_35\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"186.538954\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_31\">\n     <g id=\"line2d_36\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"189.113006\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_32\">\n     <g id=\"line2d_37\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"191.38348\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_33\">\n     <g id=\"line2d_38\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"206.776115\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_34\">\n     <g id=\"line2d_39\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"214.592166\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_35\">\n     <g id=\"line2d_40\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"220.137742\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_36\">\n     <g id=\"line2d_41\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"224.439225\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_37\">\n     <g id=\"line2d_42\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"227.953793\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_38\">\n     <g id=\"line2d_43\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"230.925318\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_39\">\n     <g id=\"line2d_44\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"233.499369\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_40\">\n     <g id=\"line2d_45\">\n      <g>\n       <use xlink:href=\"#m111e93ffd2\" x=\"235.769844\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_6\">\n     <!-- samples N -->\n     <g transform=\"translate(122.914844 174.076563) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-73\" d=\"M 2834 3397 \nL 2834 2853 \nQ 2591 2978 2328 3040 \nQ 2066 3103 1784 3103 \nQ 1356 3103 1142 2972 \nQ 928 2841 928 2578 \nQ 928 2378 1081 2264 \nQ 1234 2150 1697 2047 \nL 1894 2003 \nQ 2506 1872 2764 1633 \nQ 3022 1394 3022 966 \nQ 3022 478 2636 193 \nQ 2250 -91 1575 -91 \nQ 1294 -91 989 -36 \nQ 684 19 347 128 \nL 347 722 \nQ 666 556 975 473 \nQ 1284 391 1588 391 \nQ 1994 391 2212 530 \nQ 2431 669 2431 922 \nQ 2431 1156 2273 1281 \nQ 2116 1406 1581 1522 \nL 1381 1569 \nQ 847 1681 609 1914 \nQ 372 2147 372 2553 \nQ 372 3047 722 3315 \nQ 1072 3584 1716 3584 \nQ 2034 3584 2315 3537 \nQ 2597 3491 2834 3397 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-61\" d=\"M 2194 1759 \nQ 1497 1759 1228 1600 \nQ 959 1441 959 1056 \nQ 959 750 1161 570 \nQ 1363 391 1709 391 \nQ 2188 391 2477 730 \nQ 2766 1069 2766 1631 \nL 2766 1759 \nL 2194 1759 \nz\nM 3341 1997 \nL 3341 0 \nL 2766 0 \nL 2766 531 \nQ 2569 213 2275 61 \nQ 1981 -91 1556 -91 \nQ 1019 -91 701 211 \nQ 384 513 384 1019 \nQ 384 1609 779 1909 \nQ 1175 2209 1959 2209 \nL 2766 2209 \nL 2766 2266 \nQ 2766 2663 2505 2880 \nQ 2244 3097 1772 3097 \nQ 1472 3097 1187 3025 \nQ 903 2953 641 2809 \nL 641 3341 \nQ 956 3463 1253 3523 \nQ 1550 3584 1831 3584 \nQ 2591 3584 2966 3190 \nQ 3341 2797 3341 1997 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6d\" d=\"M 3328 2828 \nQ 3544 3216 3844 3400 \nQ 4144 3584 4550 3584 \nQ 5097 3584 5394 3201 \nQ 5691 2819 5691 2113 \nL 5691 0 \nL 5113 0 \nL 5113 2094 \nQ 5113 2597 4934 2840 \nQ 4756 3084 4391 3084 \nQ 3944 3084 3684 2787 \nQ 3425 2491 3425 1978 \nL 3425 0 \nL 2847 0 \nL 2847 2094 \nQ 2847 2600 2669 2842 \nQ 2491 3084 2119 3084 \nQ 1678 3084 1418 2786 \nQ 1159 2488 1159 1978 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1356 3278 1631 3431 \nQ 1906 3584 2284 3584 \nQ 2666 3584 2933 3390 \nQ 3200 3197 3328 2828 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-70\" d=\"M 1159 525 \nL 1159 -1331 \nL 581 -1331 \nL 581 3500 \nL 1159 3500 \nL 1159 2969 \nQ 1341 3281 1617 3432 \nQ 1894 3584 2278 3584 \nQ 2916 3584 3314 3078 \nQ 3713 2572 3713 1747 \nQ 3713 922 3314 415 \nQ 2916 -91 2278 -91 \nQ 1894 -91 1617 61 \nQ 1341 213 1159 525 \nz\nM 3116 1747 \nQ 3116 2381 2855 2742 \nQ 2594 3103 2138 3103 \nQ 1681 3103 1420 2742 \nQ 1159 2381 1159 1747 \nQ 1159 1113 1420 752 \nQ 1681 391 2138 391 \nQ 2594 391 2855 752 \nQ 3116 1113 3116 1747 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6c\" d=\"M 603 4863 \nL 1178 4863 \nL 1178 0 \nL 603 0 \nL 603 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-65\" d=\"M 3597 1894 \nL 3597 1613 \nL 953 1613 \nQ 991 1019 1311 708 \nQ 1631 397 2203 397 \nQ 2534 397 2845 478 \nQ 3156 559 3463 722 \nL 3463 178 \nQ 3153 47 2828 -22 \nQ 2503 -91 2169 -91 \nQ 1331 -91 842 396 \nQ 353 884 353 1716 \nQ 353 2575 817 3079 \nQ 1281 3584 2069 3584 \nQ 2775 3584 3186 3129 \nQ 3597 2675 3597 1894 \nz\nM 3022 2063 \nQ 3016 2534 2758 2815 \nQ 2500 3097 2075 3097 \nQ 1594 3097 1305 2825 \nQ 1016 2553 972 2059 \nL 3022 2063 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-20\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-4e\" d=\"M 628 4666 \nL 1478 4666 \nL 3547 763 \nL 3547 4666 \nL 4159 4666 \nL 4159 0 \nL 3309 0 \nL 1241 3903 \nL 1241 0 \nL 628 0 \nL 628 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-73\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(52.099609 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6d\" transform=\"translate(113.378906 0)\"/>\n      <use xlink:href=\"#DejaVuSans-70\" transform=\"translate(210.791016 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6c\" transform=\"translate(274.267578 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(302.050781 0)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(363.574219 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(415.673828 0)\"/>\n      <use xlink:href=\"#DejaVuSans-4e\" transform=\"translate(447.460938 0)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"matplotlib.axis_2\">\n    <g id=\"ytick_1\">\n     <g id=\"line2d_46\">\n      <path d=\"M 51.378125 129.843924 \nL 246.678125 129.843924 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_47\">\n      <defs>\n       <path id=\"m868d4947f2\" d=\"M 0 0 \nL -3.5 0 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m868d4947f2\" x=\"51.378125\" y=\"129.843924\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_7\">\n      <!-- $\\mathdefault{10^{-4}}$ -->\n      <g transform=\"translate(20.878125 133.643142) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-2212\" d=\"M 678 2272 \nL 4684 2272 \nL 4684 1741 \nL 678 1741 \nL 678 2272 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-34\" transform=\"translate(186.855469 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_2\">\n     <g id=\"line2d_48\">\n      <path d=\"M 51.378125 92.505198 \nL 246.678125 92.505198 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_49\">\n      <g>\n       <use xlink:href=\"#m868d4947f2\" x=\"51.378125\" y=\"92.505198\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_8\">\n      <!-- $\\mathdefault{10^{-3}}$ -->\n      <g transform=\"translate(20.878125 96.304417) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-33\" transform=\"translate(186.855469 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_3\">\n     <g id=\"line2d_50\">\n      <path d=\"M 51.378125 55.166473 \nL 246.678125 55.166473 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_51\">\n      <g>\n       <use xlink:href=\"#m868d4947f2\" x=\"51.378125\" y=\"55.166473\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_9\">\n      <!-- $\\mathdefault{10^{-2}}$ -->\n      <g transform=\"translate(20.878125 58.965692) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.765625)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 39.046875) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(186.855469 39.046875) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_4\">\n     <g id=\"line2d_52\">\n      <path d=\"M 51.378125 17.827748 \nL 246.678125 17.827748 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #b0b0b0; stroke-width: 0.8; stroke-linecap: square\"/>\n     </g>\n     <g id=\"line2d_53\">\n      <g>\n       <use xlink:href=\"#m868d4947f2\" x=\"51.378125\" y=\"17.827748\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_10\">\n      <!-- $\\mathdefault{10^{-1}}$ -->\n      <g transform=\"translate(20.878125 21.626966) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(0 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0.684375)\"/>\n       <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(128.203125 38.965625) scale(0.7)\"/>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(186.855469 38.965625) scale(0.7)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_5\">\n     <g id=\"line2d_54\">\n      <defs>\n       <path id=\"m6b91654190\" d=\"M 0 0 \nL -2 0 \n\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"144.702496\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_6\">\n     <g id=\"line2d_55\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"141.084\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_7\">\n     <g id=\"line2d_56\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"138.127473\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_8\">\n     <g id=\"line2d_57\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"135.627765\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_9\">\n     <g id=\"line2d_58\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"133.46242\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_10\">\n     <g id=\"line2d_59\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"131.55245\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_11\">\n     <g id=\"line2d_60\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"118.603847\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_12\">\n     <g id=\"line2d_61\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"112.028824\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_13\">\n     <g id=\"line2d_62\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"107.363771\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_14\">\n     <g id=\"line2d_63\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"103.745275\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_15\">\n     <g id=\"line2d_64\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"100.788748\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_16\">\n     <g id=\"line2d_65\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"98.28904\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_17\">\n     <g id=\"line2d_66\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"96.123695\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_18\">\n     <g id=\"line2d_67\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"94.213725\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_19\">\n     <g id=\"line2d_68\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"81.265122\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_20\">\n     <g id=\"line2d_69\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"74.690099\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_21\">\n     <g id=\"line2d_70\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"70.025046\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_22\">\n     <g id=\"line2d_71\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"66.406549\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_23\">\n     <g id=\"line2d_72\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"63.450023\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_24\">\n     <g id=\"line2d_73\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"60.950315\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_25\">\n     <g id=\"line2d_74\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"58.784969\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_26\">\n     <g id=\"line2d_75\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"56.874999\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_27\">\n     <g id=\"line2d_76\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"43.926397\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_28\">\n     <g id=\"line2d_77\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"37.351374\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_29\">\n     <g id=\"line2d_78\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"32.68632\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_30\">\n     <g id=\"line2d_79\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"29.067824\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_31\">\n     <g id=\"line2d_80\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"26.111297\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_32\">\n     <g id=\"line2d_81\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"23.611589\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_33\">\n     <g id=\"line2d_82\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"21.446244\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_34\">\n     <g id=\"line2d_83\">\n      <g>\n       <use xlink:href=\"#m6b91654190\" x=\"51.378125\" y=\"19.536274\" style=\"stroke: #000000; stroke-width: 0.6\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_11\">\n     <!-- integration error -->\n     <g transform=\"translate(14.798438 117.657031) rotate(-90) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-69\" d=\"M 603 3500 \nL 1178 3500 \nL 1178 0 \nL 603 0 \nL 603 3500 \nz\nM 603 4863 \nL 1178 4863 \nL 1178 4134 \nL 603 4134 \nL 603 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6e\" d=\"M 3513 2113 \nL 3513 0 \nL 2938 0 \nL 2938 2094 \nQ 2938 2591 2744 2837 \nQ 2550 3084 2163 3084 \nQ 1697 3084 1428 2787 \nQ 1159 2491 1159 1978 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1366 3272 1645 3428 \nQ 1925 3584 2291 3584 \nQ 2894 3584 3203 3211 \nQ 3513 2838 3513 2113 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-74\" d=\"M 1172 4494 \nL 1172 3500 \nL 2356 3500 \nL 2356 3053 \nL 1172 3053 \nL 1172 1153 \nQ 1172 725 1289 603 \nQ 1406 481 1766 481 \nL 2356 481 \nL 2356 0 \nL 1766 0 \nQ 1100 0 847 248 \nQ 594 497 594 1153 \nL 594 3053 \nL 172 3053 \nL 172 3500 \nL 594 3500 \nL 594 4494 \nL 1172 4494 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-67\" d=\"M 2906 1791 \nQ 2906 2416 2648 2759 \nQ 2391 3103 1925 3103 \nQ 1463 3103 1205 2759 \nQ 947 2416 947 1791 \nQ 947 1169 1205 825 \nQ 1463 481 1925 481 \nQ 2391 481 2648 825 \nQ 2906 1169 2906 1791 \nz\nM 3481 434 \nQ 3481 -459 3084 -895 \nQ 2688 -1331 1869 -1331 \nQ 1566 -1331 1297 -1286 \nQ 1028 -1241 775 -1147 \nL 775 -588 \nQ 1028 -725 1275 -790 \nQ 1522 -856 1778 -856 \nQ 2344 -856 2625 -561 \nQ 2906 -266 2906 331 \nL 2906 616 \nQ 2728 306 2450 153 \nQ 2172 0 1784 0 \nQ 1141 0 747 490 \nQ 353 981 353 1791 \nQ 353 2603 747 3093 \nQ 1141 3584 1784 3584 \nQ 2172 3584 2450 3431 \nQ 2728 3278 2906 2969 \nL 2906 3500 \nL 3481 3500 \nL 3481 434 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-72\" d=\"M 2631 2963 \nQ 2534 3019 2420 3045 \nQ 2306 3072 2169 3072 \nQ 1681 3072 1420 2755 \nQ 1159 2438 1159 1844 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1341 3275 1631 3429 \nQ 1922 3584 2338 3584 \nQ 2397 3584 2469 3576 \nQ 2541 3569 2628 3553 \nL 2631 2963 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6f\" d=\"M 1959 3097 \nQ 1497 3097 1228 2736 \nQ 959 2375 959 1747 \nQ 959 1119 1226 758 \nQ 1494 397 1959 397 \nQ 2419 397 2687 759 \nQ 2956 1122 2956 1747 \nQ 2956 2369 2687 2733 \nQ 2419 3097 1959 3097 \nz\nM 1959 3584 \nQ 2709 3584 3137 3096 \nQ 3566 2609 3566 1747 \nQ 3566 888 3137 398 \nQ 2709 -91 1959 -91 \nQ 1206 -91 779 398 \nQ 353 888 353 1747 \nQ 353 2609 779 3096 \nQ 1206 3584 1959 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-69\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(27.783203 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(91.162109 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(130.371094 0)\"/>\n      <use xlink:href=\"#DejaVuSans-67\" transform=\"translate(191.894531 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(255.371094 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(296.484375 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(357.763672 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(396.972656 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(424.755859 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(485.9375 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(549.316406 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(581.103516 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(642.626953 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(681.990234 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(720.853516 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(782.035156 0)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"line2d_84\">\n    <path d=\"M 60.255398 13.5 \nL 62.030852 13.5 \nL 62.91858 13.941281 \nL 63.806307 16.827049 \nL 64.694034 16.827049 \nL 65.581761 15.630696 \nL 66.469489 15.630696 \nL 67.357216 20.269365 \nL 68.244943 18.302376 \nL 69.13267 18.302376 \nL 70.020398 20.608816 \nL 70.908125 20.505959 \nL 71.795852 21.60183 \nL 72.68358 20.018374 \nL 73.571307 20.018374 \nL 74.459034 18.877747 \nL 75.346761 21.694834 \nL 76.234489 26.36871 \nL 77.122216 25.210548 \nL 78.89767 37.973128 \nL 79.785398 33.989082 \nL 80.673125 30.927885 \nL 81.560852 32.839881 \nL 82.44858 30.349131 \nL 83.336307 34.565135 \nL 85.111761 35.474522 \nL 85.999489 42.277403 \nL 86.887216 38.437652 \nL 87.774943 50.511619 \nL 88.66267 51.222717 \nL 89.550398 51.037302 \nL 90.438125 60.27094 \nL 91.325852 45.212327 \nL 92.21358 43.609756 \nL 93.101307 39.093588 \nL 93.989034 37.292711 \nL 94.876761 37.588537 \nL 95.764489 36.573651 \nL 96.652216 34.693615 \nL 97.539943 38.6335 \nL 98.42767 40.267626 \nL 99.315398 39.926844 \nL 100.203125 41.961457 \nL 101.090852 50.723645 \nL 101.97858 49.949891 \nL 103.754034 59.776006 \nL 104.641761 59.380134 \nL 105.529489 59.905191 \nL 106.417216 82.65244 \nL 107.304943 85.314188 \nL 108.19267 68.810343 \nL 109.080398 87.991072 \nL 109.968125 80.172502 \nL 110.855852 118.219317 \nL 111.74358 81.086054 \nL 112.631307 57.255384 \nL 113.519034 57.540068 \nL 114.406761 64.944616 \nL 115.294489 62.169112 \nL 116.182216 58.9018 \nL 117.069943 58.935774 \nL 117.95767 61.643248 \nL 119.733125 106.91093 \nL 120.620852 88.780311 \nL 121.50858 88.915192 \nL 122.396307 78.94282 \nL 123.284034 59.072606 \nL 124.171761 58.606559 \nL 125.059489 55.002154 \nL 125.947216 56.005245 \nL 126.834943 59.052329 \nL 127.72267 66.324512 \nL 128.610398 57.873652 \nL 129.498125 55.174398 \nL 130.385852 55.647379 \nL 131.27358 58.655277 \nL 132.161307 58.222144 \nL 133.049034 64.068197 \nL 133.936761 57.534056 \nL 134.824489 59.336935 \nL 135.712216 59.60053 \nL 136.599943 60.599519 \nL 137.48767 62.753729 \nL 138.375398 63.088879 \nL 139.263125 68.566276 \nL 140.150852 70.291255 \nL 141.03858 87.160558 \nL 141.926307 72.324643 \nL 142.814034 75.768799 \nL 143.701761 66.23964 \nL 144.589489 68.570844 \nL 145.477216 70.168585 \nL 146.364943 68.139562 \nL 147.25267 78.185084 \nL 148.140398 100.880446 \nL 149.028125 99.873313 \nL 149.915852 109.053785 \nL 150.80358 89.95342 \nL 151.691307 86.878385 \nL 152.579034 81.253669 \nL 154.354489 87.159238 \nL 155.242216 78.113882 \nL 156.129943 72.434195 \nL 157.01767 78.429082 \nL 157.905398 81.888922 \nL 158.793125 73.008846 \nL 159.680852 66.416604 \nL 160.56858 62.686879 \nL 161.456307 60.994532 \nL 162.344034 60.769973 \nL 163.231761 65.387815 \nL 164.119489 69.263263 \nL 165.007216 68.725129 \nL 165.894943 71.961913 \nL 166.78267 82.605441 \nL 167.670398 98.669465 \nL 168.558125 89.003922 \nL 169.445852 98.755284 \nL 170.33358 112.680103 \nL 171.221307 91.519158 \nL 172.109034 100.782367 \nL 172.996761 105.362858 \nL 173.884489 116.853272 \nL 174.772216 102.807446 \nL 175.659943 94.291076 \nL 176.54767 103.11504 \nL 177.435398 114.097713 \nL 178.323125 82.171068 \nL 179.210852 118.038027 \nL 180.09858 129.276468 \nL 180.986307 107.995242 \nL 181.874034 103.887786 \nL 182.761761 98.73603 \nL 183.649489 85.513509 \nL 184.537216 93.180217 \nL 185.424943 115.755476 \nL 186.31267 108.629879 \nL 187.200398 96.428458 \nL 188.088125 90.746774 \nL 188.975852 89.021973 \nL 189.86358 91.620598 \nL 190.751307 79.606649 \nL 191.639034 80.403802 \nL 192.526761 78.809699 \nL 193.414489 85.851585 \nL 194.302216 90.793089 \nL 195.189943 93.951043 \nL 196.07767 94.423365 \nL 196.965398 104.705957 \nL 197.853125 103.652556 \nL 198.740852 125.27564 \nL 199.62858 134.499636 \nL 200.516307 108.234247 \nL 201.404034 105.121194 \nL 202.291761 99.31393 \nL 203.179489 123.23679 \nL 204.067216 105.501078 \nL 204.954943 104.981831 \nL 205.84267 93.912359 \nL 206.730398 92.121381 \nL 207.618125 95.36851 \nL 208.505852 94.67826 \nL 209.39358 96.513916 \nL 210.281307 94.378557 \nL 211.169034 98.247522 \nL 212.056761 105.393663 \nL 212.944489 114.650548 \nL 213.832216 104.359641 \nL 214.719943 115.927938 \nL 215.60767 139.5 \nL 216.495398 121.168926 \nL 217.383125 133.46339 \nL 218.270852 117.655663 \nL 219.15858 97.43419 \nL 220.046307 102.446364 \nL 220.934034 104.300172 \nL 221.821761 109.437258 \nL 222.709489 112.684939 \nL 223.597216 125.505659 \nL 224.484943 114.219932 \nL 225.37267 130.591453 \nL 226.260398 136.5329 \nL 227.148125 122.728868 \nL 228.035852 111.526621 \nL 228.92358 106.939333 \nL 229.811307 115.65111 \nL 230.699034 125.499447 \nL 231.586761 130.781258 \nL 232.474489 126.046785 \nL 233.362216 110.36758 \nL 234.249943 102.760391 \nL 235.13767 99.238033 \nL 236.025398 102.322645 \nL 236.913125 103.344969 \nL 237.800852 108.663586 \nL 237.800852 108.663586 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square\"/>\n   </g>\n   <g id=\"line2d_85\">\n    <path d=\"M 60.255398 25.257034 \nL 237.800852 99.934485 \nL 237.800852 99.934485 \n\" clip-path=\"url(#p95e8260ecd)\" style=\"fill: none; stroke-dasharray: 5.55,2.4; stroke-dashoffset: 0; stroke: #bf00bf; stroke-width: 1.5\"/>\n   </g>\n   <g id=\"patch_3\">\n    <path d=\"M 51.378125 145.8 \nL 51.378125 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_4\">\n    <path d=\"M 246.678125 145.8 \nL 246.678125 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_5\">\n    <path d=\"M 51.378125 145.8 \nL 246.678125 145.8 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_6\">\n    <path d=\"M 51.378125 7.2 \nL 246.678125 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"legend_1\">\n    <g id=\"patch_7\">\n     <path d=\"M 111.55 46.078125 \nL 239.678125 46.078125 \nQ 241.678125 46.078125 241.678125 44.078125 \nL 241.678125 14.2 \nQ 241.678125 12.2 239.678125 12.2 \nL 111.55 12.2 \nQ 109.55 12.2 109.55 14.2 \nL 109.55 44.078125 \nQ 109.55 46.078125 111.55 46.078125 \nz\n\" style=\"fill: #ffffff; opacity: 0.8; stroke: #cccccc; stroke-linejoin: miter\"/>\n    </g>\n    <g id=\"line2d_86\">\n     <path d=\"M 113.55 20.298438 \nL 123.55 20.298438 \nL 133.55 20.298438 \n\" style=\"fill: none; stroke: #1f77b4; stroke-width: 1.5; stroke-linecap: square\"/>\n    </g>\n    <g id=\"text_12\">\n     <!-- measured MC error -->\n     <g transform=\"translate(141.55 23.798438) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-75\" d=\"M 544 1381 \nL 544 3500 \nL 1119 3500 \nL 1119 1403 \nQ 1119 906 1312 657 \nQ 1506 409 1894 409 \nQ 2359 409 2629 706 \nQ 2900 1003 2900 1516 \nL 2900 3500 \nL 3475 3500 \nL 3475 0 \nL 2900 0 \nL 2900 538 \nQ 2691 219 2414 64 \nQ 2138 -91 1772 -91 \nQ 1169 -91 856 284 \nQ 544 659 544 1381 \nz\nM 1991 3584 \nL 1991 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-64\" d=\"M 2906 2969 \nL 2906 4863 \nL 3481 4863 \nL 3481 0 \nL 2906 0 \nL 2906 525 \nQ 2725 213 2448 61 \nQ 2172 -91 1784 -91 \nQ 1150 -91 751 415 \nQ 353 922 353 1747 \nQ 353 2572 751 3078 \nQ 1150 3584 1784 3584 \nQ 2172 3584 2448 3432 \nQ 2725 3281 2906 2969 \nz\nM 947 1747 \nQ 947 1113 1208 752 \nQ 1469 391 1925 391 \nQ 2381 391 2643 752 \nQ 2906 1113 2906 1747 \nQ 2906 2381 2643 2742 \nQ 2381 3103 1925 3103 \nQ 1469 3103 1208 2742 \nQ 947 2381 947 1747 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-4d\" d=\"M 628 4666 \nL 1569 4666 \nL 2759 1491 \nL 3956 4666 \nL 4897 4666 \nL 4897 0 \nL 4281 0 \nL 4281 4097 \nL 3078 897 \nL 2444 897 \nL 1241 4097 \nL 1241 0 \nL 628 0 \nL 628 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-43\" d=\"M 4122 4306 \nL 4122 3641 \nQ 3803 3938 3442 4084 \nQ 3081 4231 2675 4231 \nQ 1875 4231 1450 3742 \nQ 1025 3253 1025 2328 \nQ 1025 1406 1450 917 \nQ 1875 428 2675 428 \nQ 3081 428 3442 575 \nQ 3803 722 4122 1019 \nL 4122 359 \nQ 3791 134 3420 21 \nQ 3050 -91 2638 -91 \nQ 1578 -91 968 557 \nQ 359 1206 359 2328 \nQ 359 3453 968 4101 \nQ 1578 4750 2638 4750 \nQ 3056 4750 3426 4639 \nQ 3797 4528 4122 4306 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-6d\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(97.412109 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(158.935547 0)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(220.214844 0)\"/>\n      <use xlink:href=\"#DejaVuSans-75\" transform=\"translate(272.314453 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(335.693359 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(374.556641 0)\"/>\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(436.080078 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(499.556641 0)\"/>\n      <use xlink:href=\"#DejaVuSans-4d\" transform=\"translate(531.34375 0)\"/>\n      <use xlink:href=\"#DejaVuSans-43\" transform=\"translate(617.623047 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(687.447266 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(719.234375 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(780.757812 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(820.121094 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(858.984375 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(920.166016 0)\"/>\n     </g>\n    </g>\n    <g id=\"line2d_87\">\n     <path d=\"M 113.55 36.478125 \nL 123.55 36.478125 \nL 133.55 36.478125 \n\" style=\"fill: none; stroke-dasharray: 5.55,2.4; stroke-dashoffset: 0; stroke: #bf00bf; stroke-width: 1.5\"/>\n    </g>\n    <g id=\"text_13\">\n     <!-- $0.2\\,N^{-1/2}$ guide -->\n     <g transform=\"translate(141.55 39.978125) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-2e\" d=\"M 684 794 \nL 1344 794 \nL 1344 0 \nL 684 0 \nL 684 794 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-Oblique-4e\" d=\"M 1081 4666 \nL 1931 4666 \nL 3219 666 \nL 4000 4666 \nL 4616 4666 \nL 3706 0 \nL 2853 0 \nL 1569 4025 \nL 788 0 \nL 172 0 \nL 1081 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-2f\" d=\"M 1625 4666 \nL 2156 4666 \nL 531 -594 \nL 0 -594 \nL 1625 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(0 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(89.910156 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-Oblique-4e\" transform=\"translate(169.768879 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-2212\" transform=\"translate(251.590338 39.046875) scale(0.7)\"/>\n      <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(310.242681 39.046875) scale(0.7)\"/>\n      <use xlink:href=\"#DejaVuSans-2f\" transform=\"translate(354.778814 39.046875) scale(0.7)\"/>\n      <use xlink:href=\"#DejaVuSans-32\" transform=\"translate(375.825299 39.046875) scale(0.7)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(423.095806 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-67\" transform=\"translate(454.882916 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-75\" transform=\"translate(518.359478 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(581.738385 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(609.521588 0.765625)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(672.99815 0.765625)\"/>\n     </g>\n    </g>\n   </g>\n  </g>\n </g>\n <defs>\n  <clipPath id=\"p95e8260ecd\">\n   <rect x=\"51.378125\" y=\"7.2\" width=\"195.3\" height=\"138.6\"/>\n  </clipPath>\n </defs>\n</svg>\n"
     },
     "metadata": {}
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-9 fw=all -->\n",
    "\n",
    "The running estimate closes in on the quadrature value as $n$\n",
    "grows, unevenly, as randomness must: the $1/\\sqrt{n}$ rate holds *in\n",
    "expectation* (about $6\\times 10^{-4}$ at $n=10^5$ for this integrand), and an\n",
    "individual prefix can land lucky, as the $n=1000$ row here does, its error a\n",
    "tenth of the predicted scale. The log-log plot shows the measured error\n",
    "scattering around the $0.2\\,N^{-1/2}$ guide line while tracking its slope of\n",
    "$-\\tfrac12$ across four decades. A second-order grid would beat that slope at\n",
    "$d=2$ (error $N^{-1}$), tie it exactly at $d=4$ ($N^{-2/4}=N^{-1/2}$), and by\n",
    "$d=6$ and $d=10$ flatten to slopes $-\\tfrac13$ and $-\\tfrac15$ and lose\n",
    "badly; for the first-order left rule the break-even is already $d=2$. Past a\n",
    "handful of dimensions, sampling is the only practical choice, which is the\n",
    "quantitative face of the curse of dimensionality.\n",
    "\n",
    "### Differentiating under the Integral Sign\n",
    "\n",
    "::: {.callout-note title=\"Assumptions behind common calculus interchanges\"}\n",
    "Several formulas in this appendix are conditional statements. The following\n",
    "sufficient conditions cover the cases used most often in machine learning:\n",
    "\n",
    "* A derivative or Jacobian formula assumes differentiability on a neighborhood\n",
    "  of the point. At a ReLU kink, autograd returns a library convention rather\n",
    "  than a classical derivative.\n",
    "* Equality of mixed partials and symmetry of the Hessian follow when the second\n",
    "  partial derivatives are continuous on a neighborhood.\n",
    "* Swapping two integrals is valid for a non-negative integrand (Tonelli) or an\n",
    "  absolutely integrable one (Fubini). Conditional convergence needs separate\n",
    "  care.\n",
    "* Passing a limit or derivative through an integral is justified when an\n",
    "  integrable dominating function controls the relevant sequence or derivative\n",
    "  throughout a neighborhood of the parameter. Moving boundaries contribute the\n",
    "  Leibniz boundary terms below.\n",
    "* A change of variables requires a one-to-one differentiable map on the region\n",
    "  of interest, or a partition into such regions, and a nonzero Jacobian\n",
    "  determinant away from negligible boundary sets.\n",
    "\n",
    "These are convenient sufficient conditions, not the weakest possible ones.\n",
    "When a model changes support, approaches a singular distribution, or contains a\n",
    "nonsmooth branch, pause and check the relevant theorem rather than applying the\n",
    "formal rule mechanically.\n",
    ":::\n",
    "\n",
    "One more interchange matters in deep learning: moving a gradient through an\n",
    "integral. When an integrand\n",
    "depends on a parameter $\\boldsymbol{\\theta}$, we will routinely write\n",
    "\n",
    "$$\n",
    "\\nabla_{\\boldsymbol{\\theta}} \\int f(x, \\boldsymbol{\\theta})\\,dx\n",
    "= \\int \\nabla_{\\boldsymbol{\\theta}} f(x, \\boldsymbol{\\theta})\\,dx .\n",
    "$$\n",
    "\n",
    "The swap needs a hypothesis. The standard sufficient condition, a consequence\n",
    "of the **dominated convergence theorem** that we also state on faith\n",
    "[@Folland.1999], is that, near the current $\\boldsymbol{\\theta}$, the\n",
    "integrands' derivatives are bounded in magnitude by a single integrable\n",
    "function of $x$. Many smooth, rapidly decaying densities used in deep learning comply, but the\n",
    "condition must be checked rather than inferred from smoothness alone. Moving\n",
    "support and increasingly sharp parameter-dependent spikes are two common ways\n",
    "it can fail; failure is not limited to those cases.\n",
    "\n",
    "When it is the *limits* of integration\n",
    "that move, the correct derivative is known in closed form. The **Leibniz\n",
    "rule** [@Folland.1999] states that\n",
    "\n",
    "$$\n",
    "\\frac{d}{d\\theta} \\int_{a(\\theta)}^{b(\\theta)} f(x, \\theta)\\,dx\n",
    "= \\int_{a(\\theta)}^{b(\\theta)} \\frac{\\partial f}{\\partial\\theta}(x,\\theta)\\,dx\n",
    "+ f\\bigl(b(\\theta), \\theta\\bigr)\\,b'(\\theta)\n",
    "- f\\bigl(a(\\theta), \\theta\\bigr)\\,a'(\\theta) .\n",
    "$$\n",
    "\n",
    "The two boundary terms are precisely what the naive swap forgets: each moving\n",
    "endpoint sweeps area in or out at the rate the limit moves, weighted by the\n",
    "integrand's value there. A density supported on $[0, \\theta]$ is the classic\n",
    "trap: dropping the boundary term produces a confidently wrong gradient.\n",
    "\n",
    "Two results we lean on later are direct applications of this swap.\n",
    "Differentiating the normalization\n",
    "$\\int p_{\\boldsymbol{\\theta}}(x)\\,dx = 1$ under the integral shows the score\n",
    "$\\nabla_{\\boldsymbol{\\theta}}\\log p_{\\boldsymbol{\\theta}}$ has mean zero, the\n",
    "seed of Fisher information (that section). This score\n",
    "is a gradient with respect to the *parameters* $\\boldsymbol{\\theta}$; it is a\n",
    "different object from the data score\n",
    "$\\nabla_{\\mathbf{x}}\\log p(\\mathbf{x})$ that score matching estimates, which is\n",
    "a gradient with respect to the data point. Applying the same swap to an expectation gives, when $g$ itself does not\n",
    "depend on $\\boldsymbol\\theta$,\n",
    "\n",
    "$$\n",
    "\\nabla_{\\boldsymbol{\\theta}}\\,\\mathbb{E}_{x\\sim p_{\\boldsymbol{\\theta}}}[g(x)]\n",
    "= \\mathbb{E}\\bigl[g(x)\\,\\nabla_{\\boldsymbol{\\theta}}\\log\n",
    "p_{\\boldsymbol{\\theta}}(x)\\bigr].\n",
    "$$\n",
    "\n",
    "This is the **score-function**, or REINFORCE, identity\n",
    "[@Williams.1992], at the heart of policy-gradient methods. If\n",
    "$g(x,\\boldsymbol\\theta)$ also depends directly on the parameter, add the\n",
    "ordinary term\n",
    "$\\mathbb E[\\nabla_{\\boldsymbol\\theta}g(x,\\boldsymbol\\theta)]$ to the right.\n",
    "\n",
    "#### Pathwise and Score-Function Gradients\n",
    "\n",
    "The score-function identity is only one way to differentiate an expectation.\n",
    "If a random variable can be generated by transforming parameter-free noise,\n",
    "\n",
    "$$\n",
    "\\boldsymbol\\epsilon\\sim q,\\qquad\n",
    "x=T_{\\boldsymbol\\theta}(\\boldsymbol\\epsilon),\n",
    "$$\n",
    "\n",
    "then\n",
    "\n",
    "$$\n",
    "\\nabla_{\\boldsymbol\\theta}\n",
    "\\mathbb E_{x\\sim p_{\\boldsymbol\\theta}}[g(x,\\boldsymbol\\theta)]\n",
    "=\n",
    "\\mathbb E_{\\boldsymbol\\epsilon\\sim q}\n",
    "\\left[\\frac{d}{d\\boldsymbol\\theta}\n",
    " g(T_{\\boldsymbol\\theta}(\\boldsymbol\\epsilon),\\boldsymbol\\theta)\\right].\n",
    "$$\n",
    "\n",
    "This is the **pathwise**, or **reparameterization**, gradient. For example,\n",
    "$X\\sim\\mathcal N(\\mu,\\sigma^2)$ can be written\n",
    "$X=\\mu+\\sigma\\epsilon$ with $\\epsilon\\sim\\mathcal N(0,1)$; automatic\n",
    "differentiation then follows an ordinary differentiable path from $\\mu$ and\n",
    "$\\sigma$ to the sample. Variational autoencoders use exactly this device to\n",
    "backpropagate through a Gaussian latent variable.\n",
    "\n",
    "The two estimators are both Monte Carlo estimators of the same derivative, but\n",
    "their variances can differ sharply. For\n",
    "$L(\\theta)=\\mathbb E_{X\\sim\\mathcal N(\\theta,1)}[X^2]=\\theta^2+1$, the true\n",
    "derivative is $2\\theta$. With $X=\\theta+\\epsilon$, one sample contributes\n",
    "$2X$ to the pathwise estimator. The score of the normal location family is\n",
    "$\\partial_\\theta\\log p_\\theta(X)=X-\\theta=\\epsilon$, so REINFORCE contributes\n",
    "$X^2\\epsilon$. Subtracting any baseline $b$ that does not depend on $X$\n",
    "leaves the expectation unchanged because $\\mathbb E[\\epsilon]=0$; choosing a\n",
    "good baseline can remove substantial variance."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-integral-stochastic-gradient-estimators"
   },
   "source": [
    "import numpy as onp\n",
    "rng = onp.random.default_rng(7)\n",
    "theta, n = 1.5, 200000\n",
    "eps = rng.standard_normal(n)\n",
    "x = theta + eps\n",
    "pathwise = 2 * x\n",
    "score = x**2 * eps\n",
    "b_opt = theta**2 + 3                         # minimizes Var[(X^2-b) epsilon]\n",
    "score_baseline = (x**2 - b_opt) * eps\n",
    "\n",
    "def report(name, samples):\n",
    "    se = samples.std(ddof=1) / onp.sqrt(len(samples))\n",
    "    print(f'{name:22s} mean={samples.mean():.4f}  '\n",
    "          f'sample sd={samples.std(ddof=1):.3f}  MC se={se:.4f}')\n",
    "\n",
    "print(f'true derivative = {2 * theta:.4f}')\n",
    "report('pathwise', pathwise)\n",
    "report('score function', score)\n",
    "report('score + baseline', score_baseline)"
   ],
   "execution_count": 9,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "true derivative = 3.0000\npathwise               mean=3.0019  sample sd=1.998  MC se=0.0045\nscore function         mean=2.9974  sample sd=7.162  MC se=0.0160\nscore + baseline       mean=2.9924  sample sd=4.878  MC se=0.0109\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-10 fw=all -->\n",
    "\n",
    "All three sample means agree with $2\\theta$ within Monte Carlo error. The\n",
    "pathwise estimate is much less noisy on this smooth problem; the baseline\n",
    "helps REINFORCE but does not generally close the gap. This is a common pattern,\n",
    "not a universal ordering: variance depends on the model, transformation, and\n",
    "integrand.\n",
    "\n",
    "The score estimator remains valuable because it does not differentiate through\n",
    "the sampled value. It therefore works for Bernoulli and categorical choices,\n",
    "where an exact differentiable reparameterization is generally unavailable.\n",
    "Continuous relaxations such as Gumbel--softmax\n",
    "[@Jang.Gu.Poole.2017; @Maddison.Mnih.Teh.2017] create a differentiable\n",
    "surrogate; straight-through estimators use a discrete forward value and a\n",
    "surrogate backward derivative. These can reduce variance, but they are biased\n",
    "for the original discrete objective unless a separate argument establishes\n",
    "otherwise.\n",
    "\n",
    "Both identities inherit the assumptions of differentiation under the integral\n",
    "sign. The score form also assumes the relevant support does not move with\n",
    "$\\boldsymbol\\theta$, or that the missing boundary terms are included. The\n",
    "pathwise form requires a differentiable transformation and an integrable bound\n",
    "for its derivative. In practice, “differentiate through randomness” therefore\n",
    "means choosing an estimator, checking its assumptions, and reporting its\n",
    "variance or Monte Carlo standard error—not merely asking an autodiff library\n",
    "for a gradient.\n",
    "\n",
    "## Summary\n",
    "\n",
    "* The **definite integral** $\\int_a^b f\\,dx$ is the limit of Riemann rectangle\n",
    "  sums the equation: the signed area under $f$, with negative\n",
    "  contributions where $f<0$ or when the limits run backwards.\n",
    "* The **fundamental theorem of calculus** the equation says the area-so-far\n",
    "  function $F(x)=\\int_a^x f$ has derivative $F'=f$. Integration is therefore\n",
    "  differentiation reversed: $\\int_a^b f = G(b)-G(a)$ for any antiderivative $G$.\n",
    "* **Improper integrals** extend the definition to infinite domains as a limit;\n",
    "  $\\int_1^\\infty x^{-p}\\,dx$ converges exactly when $p>1$.\n",
    "* **Integration by parts** the equation converts $\\int u\\,v'$ into a\n",
    "  boundary term minus $\\int u'\\,v$: the product rule run backwards, and the\n",
    "  standard device for moving a derivative onto the factor that can absorb it\n",
    "  (score matching, expectations of products).\n",
    "* **Change of variables** multiplies by the local stretch: $du/dx$ in one\n",
    "  dimension the equation, the Jacobian determinant\n",
    "  $|\\det D\\boldsymbol{\\phi}|$ for a $C^1$-diffeomorphism in many\n",
    "  the equation, the same volume-scaling read off the\n",
    "  determinant.\n",
    "* **Fubini's theorem** evaluates a multiple integral as iterated single integrals\n",
    "  in either order, for absolutely integrable functions\n",
    "  ($\\int_U |f|\\,d\\mathbf{x} < \\infty$).\n",
    "* **Integration is the language of continuous probability:** a density is a\n",
    "  normalized non-negative function $\\int p = 1$, and an expectation is the integral\n",
    "  $\\int x\\,p\\,dx$. When it has no closed form, Monte Carlo estimates it by\n",
    "  sampling, converging at the dimension-free rate $1/\\sqrt{n}$, unlike a grid,\n",
    "  whose cost explodes with dimension.\n",
    "* Gradients of expectations admit two principal Monte Carlo estimators. The\n",
    "  **score-function estimator** differentiates the log-density and works for\n",
    "  discrete samples but can have high variance; the **pathwise estimator**\n",
    "  differentiates a transformation of parameter-free noise and is usually\n",
    "  preferred when such a smooth reparameterization exists. Baselines reduce\n",
    "  score-estimator variance without changing its expectation.\n",
    "\n",
    "## Exercises\n",
    "\n",
    "1. For $X\\sim\\mathcal N(\\theta,1)$ and $g(X)=X^4$, derive the pathwise and\n",
    "   score-function gradient estimators. Verify analytically that both have mean\n",
    "   $4\\theta^3+12\\theta$, then compare their Monte Carlo variances in NumPy.\n",
    "2. Show that replacing $g(X)$ by $g(X)-b$ in\n",
    "   the equation leaves the expected gradient\n",
    "   unchanged when $b$ is independent of $X$. Derive the scalar baseline that\n",
    "   minimizes the estimator's second moment.\n",
    "3. Evaluate $\\int_1^2 \\tfrac{1}{x}\\,dx$ using an antiderivative, then confirm with\n",
    "   a Riemann sum.\n",
    "4. Use the change-of-variables formula to integrate\n",
    "   $\\int_0^{\\sqrt{\\pi}} x\\sin(x^2)\\,dx$.\n",
    "5. Compute $\\int_{[0,1]^2} xy\\,dx\\,dy$ by Fubini's theorem.\n",
    "6. For which $p$ does $\\int_1^\\infty x^{-p}\\,dx$ converge? Verify the boundary\n",
    "   $p=1$ numerically by watching the partial integrals as the upper limit grows.\n",
    "7. Find the constant $c$ that makes $c\\,e^{-x^2}$ a probability density on\n",
    "   $\\mathbb{R}$, then compute its mean and $\\mathbb{E}[X^2]$.\n",
    "8. Let $f(x,y) = (x^2-y^2)/(x^2+y^2)^2$. Compute the iterated integrals\n",
    "   $\\int_0^1\\!\\bigl(\\int_0^1 f\\,dx\\bigr)dy$ and $\\int_0^1\\!\\bigl(\\int_0^1\n",
    "   f\\,dy\\bigr)dx$ over $[0,1]^2$ and confirm they equal $-\\tfrac{\\pi}{4}$ and\n",
    "   $+\\tfrac{\\pi}{4}$. Why does this not contradict Fubini's theorem? (Hint:\n",
    "   $\\frac{\\partial}{\\partial x}\\frac{-x}{x^2+y^2} = f(x,y)$; the singularity at the\n",
    "   origin breaks absolute integrability.)\n",
    "9. Estimate $\\int_0^1 e^{-x^2}\\,dx$ by Monte Carlo (average $e^{-x_i^2}$ over\n",
    "   $x_i$ uniform on $[0,1]$) and compare with a Riemann sum.\n",
    "10. Use integration by parts the equation twice to evaluate\n",
    "   $\\int_0^\\infty u^2 e^{-u}\\,du$, the second moment behind the exponential\n",
    "   distribution's variance in that section. The integral you\n",
    "   have computed is $\\Gamma(3)$, an instance of the **Gamma function**\n",
    "   $\\Gamma(t) = \\int_0^\\infty u^{t-1}e^{-u}\\,du$; the same\n",
    "   integration-by-parts step, run once at general $t$, gives the recursion\n",
    "   $\\Gamma(t+1) = t\\,\\Gamma(t)$ and hence $\\Gamma(n+1) = n!$, the function that\n",
    "   names the Gamma distribution of that section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-integral-calculus-md-11 fw=jax -->\n",
    "\n",
    "[Discussions](https://d2l.discourse.group/t/1093)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "d2l-jax",
   "display_name": "d2l (JAX)",
   "language": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
