Monte Carlo integration

Stephan Schmidt

March 21, 2024

1 Monte Carlo Integration

Consider a definite integral in the form:

Z ∞
Ex∼p(x) {g(x)} = g(x) · p(x)dx (1)
−∞

We can approximate this integral using samples from p(x):

N
1 X
Ex∼p(x) {g(x)} ≈ g(xi ), where xi ∼ p(x) (2)
N
i=1

If the domain of p(x) is not defined over the whole real axis (e.g. uniform distribution with
finite a and b, Gamma distribution, Beta distribution) then we can easily change the integration
bounds since these functions are 0 outside of the domain. Let’s consider a uniform distribution
U (x; a, b) with a domain a ≤ x ≤ b. This means that
Z ∞
Ex∼p(x) {g(x)} = g(x) · p(x)dx (3)
−∞

can be written as Z ∞
Ex∼p(x) {g(x)} = g(x) · U (x; a, b)dx (4)
−∞
which can be simplified to
Z b
Ex∼U (x;a,b) {g(x)} = g(x) · U (x; a, b)dx (5)
a

1
Since U (x; a, b) = b−a for all a ≤ x ≤ b,
Z b
1
Ex∼U (x;a,b) {g(x)} = g(x) · dx (6)
a b−a
Z b
1
Ex∼U (x;a,b) {g(x)} = g(x)dx (7)
b−a a
We can estimate this using Monte Carlo integration
N
1 X
Ex∼U (x;a,b) {g(x)} ≈ g(xi ), where xi ∼ U (x; a, b) (8)
N
i=1

1
We can also use an indicator function that is similar to a uniform distribution. An indicator
function I(x; a, b) is a function over x with two parameters a and b. The indicator function has
the following property:
I(x; a, b) = 1, if a ≤ x ≤ b (9)
and
I(x; a, b) = 0, if a < x or x > b (10)
Hence, the indicator function is the uniform distribution without the normalisation constant
1/(b − a).

We can use this indicator function to calculate the following integral:

Z c
P (x ≤ c) = p(x)dx (11)
−∞

The domain of the integral can be changed using the indicator function
Z ∞
P (x ≤ c) = p(x)I(x; −∞, c)dx (12)
−∞

This is because the indicator function I(x; −∞, c) is only non-zero between −∞ and c. As a
result, we can calculate P (x ≤ c) using Monte Carlo integration:

P (x ≤ c) = Ex∼p(x) {I(x; −∞, c)} (13)

N
1 X
≈ I(x; −∞, c), where xi ∼ p(x) (14)
N
i=1
PN 1
i.e. i I(x; −∞, c) calculates the number of samples inside the domain −∞ ≤ x ≤ c and N
forces the quantity to be between 0 and 1. This means that we calculate the fraction of samples
that is inside the domain −∞ ≤ x ≤ c and this is our estimate of the probability P (x ≤ c).

Additional examples:
N
1 X
Ex∼p(x) {x} ≈ xi , where xi ∼ p(x) (15)
N
i
N
1 X 2
Ex∼p(x) x2 ≈

xi , where xi ∼ p(x) (16)
N
i
  N
q(x) 1 X q(xi )
Ex∼p(x) ≈ , where xi ∼ p(x) (17)
g(x) N g(xi )
i
N
1 X
P (x2 ≤ c) = I(x2i ; −∞, c), where xi ∼ p(x) (18)
N
i

Lastly, the αth percentile of g(x) can be calculated with the empirical percentile function of the
samples [g(x1 ), g(x2 ), g(x3 ), . . . , g(xN )].

2
2 Exercises

2.1 Exercise 1
Consider the following model:
y =w+ϵ (19)
where
ϵ ∼ N 0, 22


(20)
and
w ∼ N 10, 12


(21)

Use Monte Carlo integration to answer the following questions:

1. E {w}

2. E {ϵ}

3. E {y}

4. E y 2


5. E {sin(y)}

6. 95th percentile of y

7. 95th percentile of y 2

8. P (w > 11)

9. P (y > 11)

If you know what the analytical answer is, you can compare it against your Monte Carlo esti-
mates.

2.2 Exercise 2
Consider the following integral: Z ∞
p(x)q(x)dx (22)
−∞

where p(x) = N (x|1, 22 ) and q(x) = U(x; 1, 3)1 .

This integral needs to be estimated using Monte Carlo integration.

• Use scipy.integrate.quad over the domain −10 ≤ x ≤ 10.

• Perform MC integration using samples from p(x).

• Perform MC integration using samples from q(x).

1
This is a uniform distribution between 1 and 3.

3
3 Answers

3.1 Exercise 1
Consider the following model:
y =w+ϵ (23)
where
ϵ ∼ N 0, 22


(24)
and
w ∼ N 10, 12


(25)

Since the MC estimates are random, the mean value µ and six standard deviations 6 · σ over
1000 runs are reported. The number of samples used for integration is N = 104 . This means if
you use N = 104 , your answer should be inside the reported range:

1. E {w} = 10.000122543911267 ± 0.061506466688857356

2. E {ϵ} = 0.0008070484213724224 ± 0.11989848950245952

3. E {y} = 10.00092959233264 ± 0.13472165779098402

4. E y 2 = 105.02004394008677 ± 2.7312467026892273


5. E {sin(y)} = −0.04474594507631378 ± 0.0432528124564715

6. 95th percentile of y = 13.678401371972711 ± 0.2807707587566037

7. 95th percentile of y 2 = 187.10085434033346 ± 7.68105005655212

8. P (w > 11) = 0.15880669999999997 ± 0.02233649802363835

9. P (y > 11) = 0.3276443 ± 0.028294696859305635

Note: The analytical distribution of y is given by

y ∼ N (10, 5) (26)

The analytical answers are given by

1. E {y} = 10

2. 95th percentile of y = 13.678004522900572

3. P (y > 11) = 0.32736042300928847

3.2 Exercise 2
• Perform MC integration using samples from p(x). Answer: ≈ 0.1706.

• Perform MC integration using samples from q(x). Answer: ≈ 0.1706.

Both should be equal.