4 Estimating Integrals

Monte Carlo integration transforms the problem of computing definite integrals into a sampling problem, making it particularly powerful for high-dimensional integration where traditional methods fail.

Consider the problem of estimating the integral \(\ell = \int_a^b f(x) \, dx\). The key insight is to rewrite this integral as an expectation. If \(X \sim \text{Uniform}(a, b)\), then: \[\int_a^b f(x) \, dx = (b-a) \mathbb{E}[f(X)]\]

This reformulation allows us to estimate the integral using sample averages. Given \(N\) independent samples \(X_1, \ldots, X_N \sim \text{Uniform}(a, b)\), our Monte Carlo estimator is: \[\hat{\ell}_N = \frac{b - a}{N} \sum_{i=1}^{N} f(X_i)\]

Geometric Interpretation

We’re approximating the area under \(f(x)\) by averaging function values at random points and scaling by interval width.

4.1 Statistical Properties

The Monte Carlo integration estimator possesses several important statistical properties. First, it is unbiased: \[\begin{aligned} \mathbb{E}[\hat{\ell}_N] &= \frac{b - a}{N} \sum_{i=1}^{N} \mathbb{E}[f(X_i)] = \frac{b - a}{N} \sum_{i=1}^{N} \int_a^b f(x) \frac{1}{b - a} \, dx = \ell \end{aligned}\]

The variance of our estimator is: \[\text{Var}(\hat{\ell}_N) = \frac{(b - a)^2}{N} \cdot \sigma^2_f\] where \(\sigma^2_f = \text{Var}(f(X))\) is the variance of \(f(X)\) under the uniform distribution.

Key Properties

Standard Error: \(\text{SE} = \frac{(b-a)\sigma_f}{\sqrt{N}}\)
Convergence Rate: \(O(1/\sqrt{N})\) regardless of dimension
CLT: \(\sqrt{N}(\hat{\ell}_N - \ell) \xrightarrow{d} \mathcal{N}(0, (b-a)^2\sigma^2_f)\)

4.2 Function Characteristics and Performance

Monte Carlo integration’s performance depends heavily on the function’s characteristics. While the convergence rate is theoretically independent of function complexity, practical performance varies significantly with function smoothness and oscillatory behavior.

Functions with high-frequency oscillations present particular challenges for Monte Carlo methods because random sampling may miss important features or fail to adequately capture rapid variations. This example compares Monte Carlo integration performance on two trigonometric functions with dramatically different oscillation frequencies, illustrating how function characteristics affect integration accuracy and convergence behavior.

4.2.1 Example: Low vs. High Frequency Functions

The following example integrates two functions over the interval \([0, \pi]\):

Smooth decay: \(f_1(x) = \sin(x)\) — a smooth, slowly varying function
Oscillatory: \(f_2(x) = 25\sin(25x)\) — a rapidly oscillating function with many periods

VM = window.VM
vmTheme = Generators.observe(notify => VM.plotting.onThemeChange(notify))
chartColors = VM.plotting.colors(vmTheme)

intSmoothPlot = VM.plotting.persistentPlot({height: "300px"})
intOscPlot = VM.plotting.persistentPlot({height: "300px"})
intEstimatePlot = VM.plotting.persistentPlot({height: "380px"})
intErrorPlot = VM.plotting.persistentPlot({height: "380px"})

// The two integrands on [0, pi]; both integrate to exactly 2.
intFunctions = [
  {name: "Smooth", title: "Smooth: sin(x)", f: x => Math.sin(x), trueValue: 2, color: "fn", marker: "circle"},
  {name: "Oscillatory", title: "Oscillatory: 25 sin(25x)", f: x => 25 * Math.sin(25 * x), trueValue: 2, color: "alt", marker: "square"}
]

// (b - a) times the average of f at n uniform points in [a, b].
function monteCarloIntegrate(f, a, b, n, rng) {
  let sum = 0
  for (let i = 0; i < n; i++) sum += f(a + (b - a) * rng())
  return (b - a) * sum / n
}

// n = 10^2 ... 10^5, 20 sizes evenly spaced on a log scale.
intSampleSizes = {
  const sizes = []
  for (let i = 0; i < 20; i++) sizes.push(Math.floor(Math.pow(10, 2 + 3 * i / 19)))
  return sizes
}

intResults = {
  const rng = VM.sampling.seededRandom(0)
  const results = []
  for (const fn of intFunctions) {
    const estimates = [], errors = []
    for (const n of intSampleSizes) {
      const estimate = monteCarloIntegrate(fn.f, 0, Math.PI, n, rng)
      estimates.push(estimate)
      errors.push(Math.abs(estimate - fn.trueValue))
    }
    // 25 independent estimates at n = 10,000, for the spread of the estimator.
    const repeats = []
    for (let run = 0; run < 25; run++) repeats.push(monteCarloIntegrate(fn.f, 0, Math.PI, 10000, rng))
    results.push({...fn, estimates, errors, repeat: VM.distributions.sampleStats(repeats, {ddof: 0})})
  }
  return results
}

function intFunctionFigure(fn, plot) {
  const xs = [], ys = []
  for (let i = 0; i <= 2000; i++) {
    const x = Math.PI * i / 2000
    xs.push(x)
    ys.push(fn.f(x))
  }
  const color = chartColors[fn.color]
  return plot([{x: xs, y: ys, type: "scatter", mode: "lines", fill: "tozeroy", line: {color, width: 2},
    fillcolor: VM.plotting.alpha(color, 0.3)}], {title: {text: fn.title}, margin: {t: 40}, showlegend: false, xaxis: {title: "x"}})
}

intFunctionFigure(intFunctions[0], intSmoothPlot)
intFunctionFigure(intFunctions[1], intOscPlot)

intEstimateTraces = {
  const traces = []
  for (const r of intResults) {
    const color = chartColors[r.color]
    traces.push({x: intSampleSizes, y: r.estimates, type: "scatter", mode: "lines+markers", name: r.name,
      line: {color}, marker: {color, symbol: r.marker}})
  }
  traces.push({x: [intSampleSizes[0], intSampleSizes[19]], y: [2, 2], type: "scatter", mode: "lines",
    name: "True value", line: {color: chartColors.ink, dash: "dash"}})
  return traces
}

VM.plotting.plotWithLegend(intEstimatePlot, intEstimateTraces, {title: {text: "Estimate vs sample size"}, margin: {t: 40},
  xaxis: {title: "Sample size", type: "log"}, yaxis: {title: "Integral estimate"}})

intErrorTraces = {
  const traces = []
  for (const r of intResults) {
    const color = chartColors[r.color]
    traces.push({x: intSampleSizes, y: r.errors, type: "scatter", mode: "lines+markers", name: r.name,
      line: {color, width: 2}, marker: {color, symbol: r.marker}})
  }
  return traces
}

VM.plotting.plotWithLegend(intErrorPlot, intErrorTraces, {title: {text: "Error vs sample size"}, margin: {t: 40},
  xaxis: {title: "Sample size", type: "log"}, yaxis: {title: "Absolute error", type: "log"}})

{
  const rows = []
  const smoothSd = intResults[0].repeat.sd
  for (const r of intResults) {
    let ratio = "–"
    if (r.name === "Oscillatory") ratio = (r.repeat.sd / smoothSd).toFixed(4)
    rows.push([r.name, r.trueValue.toFixed(4), r.repeat.mean.toFixed(4),
      Math.abs(r.repeat.mean - r.trueValue).toFixed(4), r.repeat.sd.toFixed(4), ratio])
  }
  return VM.ui.renderTable({html, rows, filename: "smooth-vs-oscillatory.csv",
    headers: ["Function", "True value", "Mean est. (25 runs, n = 10,000)", "Error", "Std dev", "Std ratio"]})
}

The contrast between the smooth \(\sin(x)\) and rapidly oscillating \(25\sin(25x)\) demonstrates why integration difficulty varies. The high-frequency function has many more sign changes and local extrema that random sampling must capture.

The log-log error plot reveals that while both functions follow the expected \(O(1/\sqrt{n})\) convergence rate, the high-frequency function consistently maintains higher absolute errors across all sample sizes. The variance ratio indicates that the oscillatory function’s variance is around 70 times higher than the smooth function’s variance, illustrating how oscillations increase uncertainty in Monte Carlo estimates.

4.2.2 Remarks

Smooth functions are well-suited for Monte Carlo integration, converging quickly with relatively few samples.
Highly oscillatory functions present challenges, requiring larger sample sizes and exhibiting higher variance.
For oscillatory integrals, consider alternative approaches such as:
- Stratified sampling to ensure adequate coverage of oscillation periods
- Importance sampling with distributions that account for function behavior
- Transformation techniques to reduce oscillations.

This example demonstrates why understanding your function’s characteristics is crucial for effective Monte Carlo integration.

4.2.3 Variance Bounds Using Derivatives

For smooth functions, we can bound the variance using derivative information.

Theorem (Variance Bound via Derivatives)

If \(g\) is differentiable on \([a,b]\) with derivative \(g'\), then: \[\text{Var}(g(X)) \leq (b-a)^2 \cdot \max_{x \in [a,b]} |g'(x)|^2\]

Proof. Let \(M = \max_{x \in [a,b]} |g'(x)|\) and \(\mu = \mathbb{E}[g(X)]\).

Since \(g\) is continuous on the compact interval \([a,b]\), it attains its minimum and maximum values. Let \(g_{\min} = \min_{x \in [a,b]} g(x)\) and \(g_{\max} = \max_{x \in [a,b]} g(x)\).

By the mean value theorem, for any two points in \([a,b]\): \[g_{\max} - g_{\min} \leq M(b - a)\]

Since \(\mu = \mathbb{E}[g(X)]\) lies between \(g_{\min}\) and \(g_{\max}\), we have for any \(x \in [a,b]\): \[|g(x) - \mu| \leq g_{\max} - g_{\min} \leq M(b - a)\]

Squaring both sides and taking expectations: \[\text{Var}(g(X)) = \mathbb{E}[(g(X) - \mu)^2] \leq [M(b - a)]^2\]

This gives us the desired bound.

Note that this is a very weak bound and not useful in practice. There are no universal techniques that give us good bounds on \(\text{Var}(g(X))\). In practice, we rely on the heuristic that oscillatory functions are high variance.

4.3 Multidimensional Integration

One of Monte Carlo integration’s greatest advantages becomes apparent in higher dimensions. While traditional numerical integration methods (like the trapezoidal rule or Simpson’s rule) suffer from the “curse of dimensionality”—requiring exponentially more grid points as dimensions increase—Monte Carlo methods maintain their \(O(1/\sqrt{n})\) convergence rate regardless of dimension.

In \(d\) dimensions, a grid-based method with \(m\) points per dimension requires \(m^d\) total evaluations. For even modest accuracy in high dimensions, this becomes computationally prohibitive. Monte Carlo integration, however, uses the same number of random samples regardless of dimension, making it the method of choice for high-dimensional problems in physics, finance, and machine learning.

4.3.1 Example: Gaussian-like Function in Multiple Dimensions

The following example demonstrates Monte Carlo integration of the function \(f(\mathbf{x}) = e^{-\|\mathbf{x}\|^2}\) over the unit hypercube \([0,1]^d\) for dimensions \(d = 1, 2, 3, 4, 5\).

intErrorDimPlot = VM.plotting.persistentPlot({height: "380px"})
intSdDimPlot = VM.plotting.persistentPlot({height: "380px"})

// Integrate f(x) = exp(-|x|²) over [0, 1]^d, d = 1..5. The exact value is
// (∫_0^1 e^{-t²} dt)^d = (sqrt(pi)/2 · erf(1))^d, and erf(1) = 2 Φ(√2) - 1.
intDimensions = {
  const oneDim = Math.sqrt(Math.PI) / 2 * (2 * VM.distributions.normalCdf(Math.SQRT2) - 1)
  const sizes = []
  for (let i = 0; i < 15; i++) sizes.push(Math.floor(Math.pow(10, 2 + 1.5 * i / 14)))
  const runs = 30
  const rng = VM.sampling.seededRandom(0)
  const results = []
  for (let d = 1; d <= 5; d++) {
    const exact = Math.pow(oneDim, d)
    const relErrors = [], sds = []
    for (const n of sizes) {
      const estimates = []
      for (let run = 0; run < runs; run++) {
        let sum = 0
        for (let i = 0; i < n; i++) {
          let normSq = 0
          for (let k = 0; k < d; k++) {
            const x = rng()
            normSq += x * x
          }
          sum += Math.exp(-normSq)
        }
        estimates.push(sum / n)
      }
      const stats = VM.distributions.sampleStats(estimates, {ddof: 0})
      relErrors.push(Math.abs(stats.mean - exact) / exact)
      sds.push(stats.sd)
    }
    results.push({d, relErrors, sds})
  }
  return {sizes, results}
}

intDimColors = ["fn", "alt", "ok", "warn", "accent2"]

intErrorDimTraces = {
  const traces = []
  for (let i = 0; i < intDimensions.results.length; i++) {
    const r = intDimensions.results[i]
    const color = chartColors[intDimColors[i]]
    traces.push({x: intDimensions.sizes, y: r.relErrors, type: "scatter", mode: "lines+markers", name: `Dim ${r.d}`,
      line: {color, width: 2}, marker: {color, size: 5}})
  }
  return traces
}

intSdDimTraces = {
  const sizes = intDimensions.sizes
  const traces = []
  let firstSdTotal = 0
  for (let i = 0; i < intDimensions.results.length; i++) {
    const r = intDimensions.results[i]
    const color = chartColors[intDimColors[i]]
    firstSdTotal += r.sds[0]
    traces.push({x: sizes, y: r.sds, type: "scatter", mode: "lines+markers", name: `Dim ${r.d}`,
      line: {color, width: 2}, marker: {color, size: 5}})
  }
  // A 1/sqrt(n) line through the average first standard deviation.
  const scale = firstSdTotal / intDimensions.results.length * Math.sqrt(sizes[0])
  const reference = []
  for (const n of sizes) reference.push(scale / Math.sqrt(n))
  traces.push({x: sizes, y: reference, type: "scatter", mode: "lines", name: "∝ 1/√n",
    line: {color: chartColors.ink, dash: "dash", width: 2}})
  return traces
}

VM.plotting.plotWithLegend(intErrorDimPlot, intErrorDimTraces, {title: {text: "Relative error vs sample size"}, margin: {t: 40},
  xaxis: {title: "Number of samples", type: "log"}, yaxis: {type: "log"}})

VM.plotting.plotWithLegend(intSdDimPlot, intSdDimTraces, {title: {text: "Standard deviation vs sample size"}, margin: {t: 40},
  xaxis: {title: "Number of samples", type: "log"}, yaxis: {title: "Standard deviation", type: "log"}})

Figure 4.1: Monte Carlo advantage in higher dimensions

The example reveals several important properties of multidimensional integration:

Dimension-Independent Convergence: The relative error plots demonstrate Monte Carlo’s key advantage - convergence behavior remains consistent across all dimensions, avoiding the curse of dimensionality that plagues grid-based methods.
Theoretical Validation: The standard deviation follows the expected \(1/\sqrt{n}\) scaling (shown by the dashed reference line) regardless of dimension, confirming Monte Carlo’s theoretical foundation.
Constant Computational Cost: Unlike deterministic methods requiring \(n^d\) evaluations, Monte Carlo uses the same number of function evaluations per sample regardless of dimension.

This dimensional robustness makes Monte Carlo indispensable for high-dimensional applications in physics, statistics, and finance where traditional methods become computationally intractable.

When to Use Monte Carlo Integration

Advantages:

High dimensions: Performance doesn’t degrade with dimension
Complex domains: Easy to handle irregular integration regions
Rough functions: Works with discontinuous or highly oscillatory functions

Disadvantages:

Slow convergence: \(O(1/\sqrt{N})\) rate
Random error: Introduces stochastic uncertainty
Low dimensions: Often outperformed by deterministic methods

Practical Considerations:

Sample size: Need large \(N\) for high precision
Function smoothness: Smoother functions → lower variance → better estimates
Variance reduction: Techniques like antithetic variables can improve efficiency
Error estimation: Use sample variance to estimate uncertainty

4.4 Conclusion

Monte Carlo integration transforms integration into a sampling problem, making it invaluable for high-dimensional problems where traditional quadrature fails. While it converges slowly, its dimension-independent performance and ability to handle complex functions make it essential for modern computational statistics and machine learning applications.