Lesson 2: Monte Carlo Approximation

1 Modelling context

Definite integrals appear throughout applied mathematics, but they cannot always be evaluated analytically. Numerical quadrature handles one-dimensional integrals efficiently, but its cost grows exponentially with dimension: doubling the dimension of a grid-based quadrature multiplies the number of function evaluations by a large factor. Monte Carlo does not have this problem. Its convergence rate is 1/√N regardless of dimension, which makes it the tool of choice for high-dimensional integration.

The Monte Carlo approximation of an integral is the direct application of the estimator introduced in Lesson 1. Any definite integral can be written as an expectation with respect to a suitable probability distribution, and the sample mean of the integrand then estimates the integral.

From integral to expectation

Every Monte Carlo approximation is built on the same identity: ∫ f(x) dx = E[f(X)/p(X)] where X has density p. Choosing p is a modelling decision, and different choices lead to different estimators with different variances.

2 Learning objectives

By the end of this lesson you will be able to

Write a definite integral as an expectation with respect to a specified distribution.
Estimate an integral by Monte Carlo sampling.
Estimate a probability from samples of a random variable.
Estimate the area of a geometric region using the hit-or-miss method.
Explain the dimensional independence of the Monte Carlo convergence rate.

3 Prerequisite review

Definite integral: ∫_a^b f(x) dx, the signed area under the curve f over the interval [a, b].
Probability density: A nonnegative function p such that ∫ p(x) dx = 1. The probability that a continuous random variable X falls in a set A is ∫_A p(x) dx.
Law of large numbers: The sample mean of independent identically-distributed random variables converges to the true expectation. This is the theoretical justification for Monte Carlo.
Indicator function: 1{x ∈ A} is 1 if x ∈ A and 0 otherwise. The expectation of an indicator is a probability.

4 Mathematical formulation

Integrals as expectations

Let f be a function on an interval [a, b], and let p be a probability density on [a, b] that is positive wherever f is nonzero. Then

∫_a^b f(x) dx = ∫_a^b [ f(x) / p(x) ] p(x) dx = E[ f(X) / p(X) ]

where X ~ p. So the integral equals the expectation of the ratio f(X) / p(X). The Monte Carlo estimator is

Î_N = (1/N) Σ f(Xᵢ) / p(Xᵢ)

The simplest choice is the uniform density on [a, b]: p(x) = 1 / (b − a). This gives

Î_N = (b − a) · (1/N) Σ f(Xᵢ),    Xᵢ ~ Uniform(a, b)

In words: sample uniformly on [a, b], evaluate f at each sample, average, and multiply by the interval length. This is the "average height times width" interpretation of the integral.

Probabilities as expectations

If X has density p and A is an event, then

P(X ∈ A) = ∫_A p(x) dx = E[ 1{X ∈ A} ]

The Monte Carlo estimator is the empirical frequency:

p̂_N = (1/N) Σ 1{Xᵢ ∈ A}

This is the fraction of samples falling in A. Its standard error is √( p(1 − p) / N ), which is bounded above by 1 / (2 √N) — this bound is the origin of the "1 over 2 root N" rule of thumb for Monte Carlo estimation of probabilities.

Hit-or-miss estimation of areas

A special case of probability estimation is the estimation of areas. Suppose we want the area of a region R contained in a bounding rectangle B. Sample uniformly in B, count the fraction of samples that fall in R, and multiply by the area of B:

Area(R) ≈ Area(B) · ( number of hits / N )

The hit-or-miss method generalises to any dimension and is the basis for using Monte Carlo to estimate volumes and high-dimensional integrals.

Figure 4.3. Hit-or-miss estimation of a circular area. The fraction of uniformly distributed points that fall inside the quarter circle estimates π/4.

5 Worked example

Estimating `∫₀¹ e^(−x²) dx`

The integral ∫₀¹ e^(−x²) dx has no elementary closed form, but its value is approximately 0.7468. We will estimate it by Monte Carlo.

Step 1 — Express as an expectation.

With p(x) = 1 on [0, 1] (i.e. uniform density on the unit interval):

∫₀¹ e^(−x²) dx = E[ e^(−X²) ],    X ~ Uniform(0, 1)

Step 2 — Sample and average.

Draw X₁, …, X_N uniformly on [0, 1], compute e^(−Xᵢ²) for each, and average:

Î_N = (1/N) Σ e^(−Xᵢ²)

Step 3 — Typical results.

Monte Carlo estimate of `∫₀¹ e^(−x²) dx`
N	Estimate	Standard error
100	0.7421	0.0093
1,000	0.7465	0.0029
10,000	0.74673	0.00093
100,000	0.74681	0.00029

The estimate converges to the true value 0.74682 as N grows, with the standard error decreasing as 1/√N.

Step 4 — Compare with the trapezoidal rule.

The trapezoidal rule with 100 equally-spaced points gives the value 0.74682 to five significant figures, using 100 function evaluations. Monte Carlo with 100 samples has a standard error of about 0.009. In one dimension, deterministic quadrature is far more efficient.

Step 5 — The dimensional argument.

Now consider estimating ∫ ... ∫_D e^(−|x|²) dx over the d-dimensional unit cube. A grid with 100 points per dimension requires 100^d evaluations. For d = 10, this is 10²⁰ — infeasible. Monte Carlo with 100,000 samples has a standard error of order 0.001 regardless of d. In high dimensions, Monte Carlo wins by default.

6 Verification check

Checks to perform

Compare against an analytical value. Choose an integral with a known closed form (e.g. ∫₀¹ x² dx = 1/3) and verify the Monte Carlo estimate converges to it.
Check the standard error scaling. Confirm that the standard error decreases as 1/√N across a range of sample sizes.
Repeat with different seeds. The estimates should scatter around the true value with no systematic drift.
Verify the estimator form. For the uniform choice of p, the estimator is the sample mean of f times the interval width. Confirm this against the derivation.
Test on a degenerate integral. If f is constant, the estimator should return the constant times the interval length, with zero variance.

7 Interpretation

The Monte Carlo approximation of an integral has a clear geometric interpretation: the integral is the "average height" of the function over the interval, multiplied by the width. Sampling uniformly and averaging is a numerical realisation of this interpretation.

The dimensional argument is the key interpretive point. In one or two dimensions, grid-based quadrature is more efficient than Monte Carlo. In higher dimensions, the exponential growth of grid-based methods makes them infeasible, and Monte Carlo's dimension-independent 1/√N rate makes it the only practical choice.

Why 1/√N is both blessing and curse

The 1/√N rate is slow: reducing the standard error by a factor of 10 requires 100 times as many samples. But it is independent of dimension, which means Monte Carlo does not slow down as the problem gets harder in the way grid-based methods do. The tradeoff is worth it whenever the dimension is more than about 4.

8 Uncertainty and limitations

The accuracy of a Monte Carlo approximation depends on three factors.

Factors affecting Monte Carlo accuracy
Factor	Effect
Number of samples `N`	Standard error scales as `1/√N`
Variance of the integrand	Larger variance means larger standard error
Choice of sampling distribution `p`	Poor choice inflates variance; good choice (importance sampling) reduces it

The third factor is the subject of Unit 4.2. Importance sampling and other variance-reduction techniques choose p to minimise the variance of the estimator, often achieving order-of-magnitude improvements over naive uniform sampling.

Small samples are misleading

In low dimensions, a Monte Carlo estimate with N = 100 may look accurate by chance. The standard error for N = 100 is 10 times that for N = 10,000. Never report a Monte Carlo estimate without its standard error, and never assume a small sample is "close enough."

9 Common implementation errors

Watch for these

Forgetting the interval-length factor. The estimator (b − a) · (1/N) Σ f(Xᵢ) includes the factor (b − a). Omitting it gives an estimate of the average height, not the integral.
Using the wrong sampling distribution. If the sampling density p is not uniform, the estimator must weight by 1/p(Xᵢ). Forgetting this weight produces a biased estimate.
Confusing the interval of integration with the support of p. The sampling distribution must cover the integration region; otherwise the estimate is biased.
Reusing the same seed across different analyses. This produces identical "random" sequences and can obscure sampling variability. Use different seeds for different runs.
Ignoring the variance of the integrand. A function with high variance requires many more samples to achieve a given standard error. Estimate the variance before choosing N.
Reporting the estimate without the sample size. The estimate is meaningless without the N and the standard error.

10 Independent modelling task

Your turn

Choose three definite integrals: one with a known analytical value, one with no elementary closed form, and one in two dimensions.
For each, write the integral as an expectation with respect to a suitable sampling distribution.
Implement the Monte Carlo estimator in a spreadsheet, running each at N = 1000 and N = 10,000.
Compute the estimate, the standard error, and a 95% confidence interval for each.
Compare against the analytical value where one exists. Discuss any discrepancy.
For the two-dimensional integral, compare the Monte Carlo effort against a grid-based quadrature at the same accuracy. Determine which method is more efficient.

This exercise forms the basis of Laboratory 7 and the integration component of Assignment 4.

11 Reflection prompt

Reflect before continuing

In one dimension, deterministic quadrature is more efficient than Monte Carlo. At what dimension does Monte Carlo become preferable?

Two or three sentences. Consider the exponential cost of grid-based methods and the dimension-independent rate of Monte Carlo.

12 Lesson summary

Key takeaways

Any definite integral can be written as an expectation. The choice of sampling distribution p is a modelling decision.
The Monte Carlo estimator of an integral is the sample mean of the integrand (with weighting). For uniform sampling, this is (b − a) · mean(f(Xᵢ)).
Probabilities are estimated as empirical frequencies. The standard error is √(p(1−p)/N), bounded above by 1/(2√N).
Areas and volumes are estimated by hit-or-miss. The fraction of samples inside the region estimates the fractional area.
The 1/√N rate is dimension-independent. This is why Monte Carlo is the tool of choice for high-dimensional integration.

Coming next

Lesson 3 develops the sampling-error analysis in detail, including the construction of confidence intervals and the diagnostic methods for verifying convergence.

§ Self-quiz, practice, and homework

Practice exercises

Write the Monte Carlo estimator for ∫₀² x³ dx using uniform sampling. Compute the analytical value and compare to the estimator for N = 1000.
Write the Monte Carlo estimator for P(Z > 1.96) where Z ~ Normal(0,1). Compare the estimate to 0.025.
Use the hit-or-miss method to estimate the area of the region {(x, y) : x² + y² < 1, x > 0, y > 0} inside the unit square. State the expected value.
For the standard error of a probability estimate with p = 0.5 and N = 10,000, compute the numerical value.

Homework problems

Integral estimation. Estimate ∫₀^π sin(x) dx by Monte Carlo at N = 100, 1000, 10000, 100000. Compare against the exact value 2. Tabulate the standard errors and the confidence intervals. Write a 500-word report.
Dimensional comparison. Estimate the volume of a unit ball in 2, 3, 5, and 10 dimensions by Monte Carlo. Use a fixed sample size of 10,000 and report the standard error in each case. Discuss the practical difficulty of the high-dimensional case.
Probability estimation. Estimate P(X > 2) where X ~ Exp(1), at several sample sizes. The true value is e^{−2} ≈ 0.1353. Confirm the 1/√N scaling of the standard error.

Mini project

Choose a definite integral or a probability that has an analytical value. Estimate it by Monte Carlo at five different sample sizes. Tabulate the estimate, standard error, and 95% confidence interval at each. Produce a log-log plot of the standard error against N and confirm the slope is −1/2. Submit the spreadsheet and a one-page report. This is the basis for Laboratory 7.