Monte Carlo

Cumulative Distribution Function

Suppose that

X

is a random variable. The cumulative distribution function of

X

is the function

F_{X} : ℝ \to [0,1]

defined by

F_{X} (x) : = P (X \leq x) .

Joint Density Function

Suppose that

X

and

Y

are real-valued random variables. A function

p : ℝ^{2} \to [0, \infty)

is a joint density function for

(X, Y)

if, for every suitable set

A \subseteq ℝ^{2}

P ((X, Y) \in A) = \int_{A} p (x, y) d x d y .

Marginal Density Function

Suppose that

p (x, y)

is a joint density function. The marginal density function of

X

p_{X} (x) : = \int_{- \infty}^{\infty} p (x, y) d y .

Conditional Density Function

Suppose that

p (x, y)

is a joint density function and

p_{X} (x) > 0

. The conditional density function of

Y

given

X = x

p_{Y | X} (y | x) : = \frac{p (x, y)}{p_{X} (x)} .

Joint Density Factors as Marginal Times Conditional Density

Suppose that

p (x, y)

is a joint density. Then

p (x, y) = p_{X} (x) p_{Y | X} (y | x)

whenever

p_{X} (x) > 0

Estimator

Suppose that

Q \in ℝ

is a fixed quantity to be estimated. An estimator for

Q

is a real-valued random variable

\hat{Q}

whose values are used as approximations to

Q

. In practice,

\hat{Q}

is often a function of random samples, so its value changes with the sampled outcomes while

Q

stays fixed.

For example, $Q$ might be the true average height of every adult in a country. Measuring every adult is unrealistic, so we randomly sample some people and compute a sample average $\hat{Q}$ . The true average $Q$ is fixed, but the sample average $\hat{Q}$ is random because a different sample can give a different value.

Unbiased Estimator

Suppose that

\hat{Q}

is an estimator for

Q

. The estimator

\hat{Q}

is unbiased if

E (\hat{Q}) = Q .

Biased Estimator

Suppose that

\hat{Q}

is an estimator for a fixed

Q \in ℝ

. The bias of

\hat{Q}

is the difference between the average value of the estimator and the target value:

E (\hat{Q}) - Q .

The estimator is biased if its bias is not

0

Monte Carlo Estimator

Suppose that

D \subseteq ℝ^{n}

f : D \to ℝ

, and

p

is a density function on

D

such that

p (x) > 0

whenever

f (x) \neq 0

. If

X_{1}, \dots, X_{N}

are independent samples with density

p

, the Monte Carlo estimator for

Q : = \int_{D} f (x) d x

F_{N} : = \frac{1}{N} \sum_{i = 1}^{N} \frac{f (X_{i})}{p (X_{i})} .

Monte Carlo Estimator is Unbiased

Under the hypotheses of the Monte Carlo estimator,

E (F_{N}) = \int_{D} f (x) d x .

By the linearity of expectation for finite sums,

E (F_{N}) = E (\frac{1}{N} \sum_{i = 1}^{N} \frac{f (X_{i})}{p (X_{i})}) = \frac{1}{N} \sum_{i = 1}^{N} E (\frac{f (X_{i})}{p (X_{i})}) .

Since each

X_{i}

has density

p

, the expectation of the weighted sample is

E (\frac{f (X_{i})}{p (X_{i})}) = \int_{D} \frac{f (x)}{p (x)} p (x) d x = \int_{D} f (x) d x .

Therefore

E (F_{N}) = \frac{1}{N} \sum_{i = 1}^{N} \int_{D} f (x) d x = \int_{D} f (x) d x .

By the definition of unbiased estimator,

F_{N}

is unbiased for

Q

Uniform Monte Carlo Estimator is Unbiased

Suppose that

a, b \in ℝ

a < b

, and

f : [a, b] \to ℝ

. If

X_{1}, \dots, X_{N}

are independent samples from the uniform density on

[a, b]

, then

F_{N} : = \frac{b - a}{N} \sum_{i = 1}^{N} f (X_{i})

satisfies

E (F_{N}) = \int_{a}^{b} f (x) d x .

The uniform density on

[a, b]

p (x) = \frac{1}{b - a}

for

x \in [a, b]

. By linearity of expectation for finite sums,

E (F_{N}) = E (\frac{b - a}{N} \sum_{i = 1}^{N} f (X_{i})) = \frac{b - a}{N} \sum_{i = 1}^{N} E (f (X_{i})) .

Since each

X_{i}

has density

p

, expectation of a function of a continuous random variable gives

E (f (X_{i})) = \int_{a}^{b} f (x) p (x) d x = \int_{a}^{b} f (x) \frac{1}{b - a} d x .

Thus

E (F_{N}) = \frac{b - a}{N} \sum_{i = 1}^{N} \int_{a}^{b} f (x) \frac{1}{b - a} d x = \frac{1}{N} \sum_{i = 1}^{N} \int_{a}^{b} f (x) d x = \int_{a}^{b} f (x) d x .

Variance of an Average of Independent Identically Distributed Random Variables

Suppose that

X_{1}, \dots, X_{N}

are independent identically distributed random variables with finite variance. Then

Var (\frac{1}{N} \sum_{i = 1}^{N} X_{i}) = \frac{Var (X_{1})}{N} .

Let

S : = \sum_{i = 1}^{N} X_{i} .

By linearity of expectation,

E (S) = \sum_{i = 1}^{N} E (X_{i}) .

Therefore

S - E (S) = \sum_{i = 1}^{N} (X_{i} - E (X_{i})) .

Now apply the definition of variance:

Var (S) = E ({(\sum_{i = 1}^{N} (X_{i} - E (X_{i})))}^{2}) .

To see the square expansion, first use distributivity on an arbitrary finite sum:

{(\sum_{i = 1}^{N} Y_{i})}^{2} = (\sum_{i = 1}^{N} Y_{i}) (\sum_{j = 1}^{N} Y_{j}) = \sum_{i = 1}^{N} \sum_{j = 1}^{N} Y_{i} Y_{j} .

The terms with

i = j

are the diagonal terms

Y_{i}^{2}

. The terms with

i \neq j

occur in pairs: for each

i < j

, the double sum contains both

Y_{i} Y_{j}

and

Y_{j} Y_{i}

, and these are equal. Hence

{(\sum_{i = 1}^{N} Y_{i})}^{2} = \sum_{i = 1}^{N} Y_{i}^{2} + 2 \sum_{1 \leq i < j \leq N} Y_{i} Y_{j} .

Applying this with

Y_{i} = X_{i} - E (X_{i})

gives

{(\sum_{i = 1}^{N} (X_{i} - E (X_{i})))}^{2} = \sum_{i = 1}^{N} {(X_{i} - E (X_{i}))}^{2} + 2 \sum_{1 \leq i < j \leq N} (X_{i} - E (X_{i})) (X_{j} - E (X_{j})) .

Taking expectations and using linearity of expectation,

Var (\sum_{i = 1}^{N} X_{i}) = \sum_{i = 1}^{N} E ({(X_{i} - E (X_{i}))}^{2}) + 2 \sum_{1 \leq i < j \leq N} E ((X_{i} - E (X_{i})) (X_{j} - E (X_{j}))) .

The first expectation is

Var (X_{i})

by the definition of variance. The second expectation is

Cov (X_{i}, X_{j})

by the definition of covariance. Hence

Var (\sum_{i = 1}^{N} X_{i}) = \sum_{i = 1}^{N} Var (X_{i}) + 2 \sum_{1 \leq i < j \leq N} Cov (X_{i}, X_{j}) .

Since the random variables are independent, they are uncorrelated, so each covariance term is

0

. Since they are identically distributed, each variance is

Var (X_{1})

. Therefore

Var (\sum_{i = 1}^{N} X_{i}) = N Var (X_{1}) .

Finally, scaling by

1 / N

scales variance by

1 / N^{2}

, so

Var (\frac{1}{N} \sum_{i = 1}^{N} X_{i}) = \frac{1}{N^{2}} N Var (X_{1}) = \frac{Var (X_{1})}{N} .

Monte Carlo Estimator Variance

Under the hypotheses of the Monte Carlo estimator, define

Y : = \frac{f (X_{1})}{p (X_{1})} .

Y

has finite variance, then

Var (F_{N}) = \frac{Var (Y)}{N} .

The random variables

f (X_{i}) / p (X_{i})

are independent and identically distributed because the samples

X_{1}, \dots, X_{N}

are independent and identically distributed. Apply variance of an average of independent identically distributed random variables to

F_{N} = \frac{1}{N} \sum_{i = 1}^{N} \frac{f (X_{i})}{p (X_{i})} .

Mean Squared Error

Suppose that

\hat{Q}

is an estimator for

Q

. The mean squared error of

\hat{Q}

MSE (\hat{Q}) : = E ({(\hat{Q} - Q)}^{2}) .

Unbiased Estimator MSE is Variance

\hat{Q}

is an unbiased estimator for

Q

, then

MSE (\hat{Q}) = Var (\hat{Q}) .

Since

\hat{Q}

is unbiased,

E (\hat{Q}) = Q

. Substituting this into the mean squared error gives

MSE (\hat{Q}) = E ({(\hat{Q} - E (\hat{Q}))}^{2}),

which is exactly the variance of

\hat{Q}

Monte Carlo Efficiency

F

is an estimator with finite positive variance and expected runtime

T (F)

, its Monte Carlo efficiency is

ϵ [F] : = \frac{1}{Var (F) T (F)} .

Sample Variance

Suppose that

x_{1}, \dots, x_{N} \in ℝ

, with

N \geq 2

, and let

x : = \frac{1}{N} \sum_{i = 1}^{N} x_{i} .

The sample variance of

x_{1}, \dots, x_{N}

s^{2} : = \frac{1}{N - 1} \sum_{i = 1}^{N} {(x_{i} - x)}^{2} .

Inversion Method

Suppose that

F

is a cumulative distribution function that has an inverse

F^{- 1}

, and suppose that

U \sim unif [0,1]

. The inversion method samples from

F

by setting

X : = F^{- 1} (U) .

The inversion method is one way to turn a uniform random number into a random variable with a desired distribution. Start with $U$ , which is uniformly distributed on $[0,1]$ . If the desired cumulative distribution function is $F$ , then the interval $[0, F (x)]$ has length $F (x)$ , so $P (U \leq F (x)) = F (x)$ . Choosing $X = F^{- 1} (U)$ makes the event $X \leq x$ match the event $U \leq F (x)$ . Thus the transformation $F^{- 1}$ reshapes uniform randomness so that the output $X$ has distribution $F$ .

Discrete Inversion Method

Suppose that

x_{1}, x_{2}, \dots

are possible values with probabilities

p_{i} \geq 0

and

\sum_{i} p_{i} = 1.

Let

P_{i} : = \sum_{j = 1}^{i} p_{j} .

U \sim unif [0,1]

, the discrete inversion method returns

x_{i}

when

P_{i - 1} \leq U < P_{i},

where

P_{0} : = 0

Discrete Inversion Method Has the Target Probabilities

Under the hypotheses of the discrete inversion method,

P (X = x_{i}) = p_{i} .

The event

{X = x_{i}}

is the event

{P_{i - 1} \leq U < P_{i}}

. Since

U

is uniform on

[0,1]

, the probability of this interval is its length:

P (X = x_{i}) = P_{i} - P_{i - 1} = \sum_{j = 1}^{i} p_{j} - \sum_{j = 1}^{i - 1} p_{j} = p_{i} .

Inversion Method Samples from the Target Distribution

Suppose that

F

is an invertible cumulative distribution function,

U \sim unif [0,1]

, and

X = F^{- 1} (U)

. Then

X

has cumulative distribution function

F

For any

x \in ℝ

P (X \leq x) = P (F^{- 1} (U) \leq x) .

Since a cumulative distribution function is non-decreasing, applying

F

to both sides gives the equivalent event

{U \leq F (x)}

. Because

U \sim unif [0,1]

P (F^{- 1} (U) \leq x) = P (U \leq F (x)) = F (x) .

Therefore the cumulative distribution function of

X

F

For example, suppose we want density $p (x) = 2 x$ on $[0,1]$ . Its cumulative distribution function is $F (x) = x^{2}$ . If $U$ is uniform on $[0,1]$ , set $X = \sqrt{U}$ . Then $P (X \leq x) = P (\sqrt{U} \leq x) = P (U \leq x^{2}) = x^{2},$ so $X$ has density $2 x$ . The random number generator only produced a uniform value; the transformation $u \mapsto \sqrt{u}$ changed how often different output values occur.

Linear Function Density

Suppose that

a, b \in [0, \infty)

, not both

0

, and define

f : [0,1] \to ℝ

f (x) = (1 - x) a + x b .

The normalized density proportional to

f

p (x) = \frac{2 f (x)}{a + b} .

Inversion Sampling for a Linear Function Density

Under the hypotheses of linear function density, if

U \sim unif [0,1]

, then a sample from density

p

X = \frac{a - \sqrt{(1 - U) a^{2} + U b^{2}}}{a - b},

when

a \neq b

. Equivalently, to avoid cancellation,

X = \frac{U (a + b)}{a + \sqrt{(1 - U) a^{2} + U b^{2}}} .

The cumulative distribution function of

p

F (x) = \int_{0}^{x} \frac{2 ((1 - t) a + t b)}{a + b} d t = \frac{x (a (2 - x) + b x)}{a + b} .

Setting

U = F (X)

gives a quadratic equation in

X

. Solving the quadratic and choosing the root in

[0,1]

gives

X = \frac{a - \sqrt{(1 - U) a^{2} + U b^{2}}}{a - b} .

Multiplying numerator and denominator by

a + \sqrt{(1 - U) a^{2} + U b^{2}}

gives the equivalent form.

Change of Variables for Density Functions

Suppose that

X

is an

n

-dimensional random variable with density

p_{X}

, and suppose that

T : ℝ^{n} \to ℝ^{n}

is a differentiable bijection whose Jacobian determinant is nonzero. If

Y = T (X)

, then the density of

Y

p_{Y} (y) = \frac{p_{X} (T^{- 1} (y))}{| \det {[D T (T^{- 1} (y))]}_{E} |} .

For every suitable set

A \subseteq ℝ^{n}

P (Y \in A) = P (X \in T^{- 1} (A)) = \int_{T^{- 1} (A)} p_{X} (x) d x .

By the multivariable change of variables formula, with

y = T (x)

\int_{T^{- 1} (A)} p_{X} (x) d x = \int_{A} \frac{p_{X} (T^{- 1} (y))}{| \det {[D T (T^{- 1} (y))]}_{E} |} d y .

Since this holds for every suitable

A

, the displayed integrand is the density of

Y

One Dimensional Density Transformation

Suppose that

X

has density

p_{X}

, and

T : ℝ \to ℝ

is differentiable and one-to-one with

T^{'} (x) \neq 0

. If

Y = T (X)

, then

p_{Y} (y) = \frac{p_{X} (T^{- 1} (y))}{| T^{'} (T^{- 1} (y)) |} .

This is the one-dimensional case of change of variables for density functions, where the Jacobian determinant is

T^{'}

Polar Coordinate Density Transformation

Under the polar coordinate map

T (r, θ) = (r \cos θ, r \sin θ)

p_{R, Θ} (r, θ) = r p_{X, Y} (r \cos θ, r \sin θ) .

The absolute value of the Jacobian determinant of the polar coordinate map is

r

. Apply change of variables for density functions.

Spherical Coordinate Density Transformation

Under the spherical coordinate map

T (r, θ, ϕ) = (r \sin θ \cos ϕ, r \sin θ \sin ϕ, r \cos θ)

p_{R, Θ, Φ} (r, θ, ϕ) = r^{2} \sin θ p_{X, Y, Z} (T (r, θ, ϕ)) .

The absolute value of the Jacobian determinant of the spherical coordinate map is

r^{2} \sin θ

. Apply change of variables for density functions.

Stratified Sampling

Suppose that a domain

D

is partitioned into disjoint measurable sets

D_{1}, \dots, D_{m}

. Stratified sampling is the practice of estimating an integral over

D

by taking samples separately within each stratum

D_{i}

and combining the resulting estimates.

Stratified Monte Carlo Estimator

Suppose that

D

is partitioned into measurable strata

D_{1}, \dots, D_{m}

, with

v_{i}

denoting the measure of

D_{i}

. If

F_{i}

is a Monte Carlo estimator for

\frac{1}{v_{i}} \int_{D_{i}} f (x) d x,

then the stratified Monte Carlo estimator is

F : = \sum_{i = 1}^{m} v_{i} F_{i} .

Stratified Sampling Does Not Increase Variance

Suppose that

D

is partitioned into measurable strata

D_{1}, \dots, D_{m}

, and samples are allocated proportionally to stratum measure:

n_{i} = v_{i} N

. If

μ_{i}

and

σ_{i}^{2}

are the mean and variance of

f

under the uniform distribution on

D_{i}

, and

Q

is the mean of

f

under the uniform distribution on

D

, then the stratified estimator has variance

Var (F) = \frac{1}{N} \sum_{i = 1}^{m} v_{i} σ_{i}^{2},

while the unstratified estimator has variance

\frac{1}{N} (\sum_{i = 1}^{m} v_{i} σ_{i}^{2} + \sum_{i = 1}^{m} v_{i} {(μ_{i} - Q)}^{2}) .

Hence stratification does not increase variance.

For stratum

D_{i}

, averaging

n_{i}

independent samples gives variance

σ_{i}^{2} / n_{i}

by variance of an average of independent identically distributed random variables. The stratified estimator is

F = \sum_{i} v_{i} F_{i}

, so independence of samples across strata gives

Var (F) = \sum_{i = 1}^{m} v_{i}^{2} Var (F_{i}) = \sum_{i = 1}^{m} v_{i}^{2} \frac{σ_{i}^{2}}{n_{i}} .

Substituting

n_{i} = v_{i} N

gives the displayed stratified variance.

For unstratified sampling, the total variance decomposes into the average within-stratum variance plus the variance of the stratum means:

Var (f (X)) = \sum_{i = 1}^{m} v_{i} σ_{i}^{2} + \sum_{i = 1}^{m} v_{i} {(μ_{i} - Q)}^{2} .

Averaging

N

independent samples divides this by

N

. Since each term

v_{i} {(μ_{i} - Q)}^{2} \geq 0

, the stratified variance is no larger.

Importance Sampling

Importance sampling estimates an integral

\int_{D} f (x) d x

by drawing samples from a density

p

and using the weighted terms

\frac{f (X)}{p (X)} .

The density

p

should be positive wherever

f

is nonzero.

For example, to estimate $\int_{0}^{1} x^{2} d x$ , ordinary uniform sampling uses $X \in [0,1]$ with density $p (x) = 1$ , so each weighted term is $X^{2}$ . Importance sampling could instead use the density $p (x) = 2 x$ , which samples more often near $1$ , where $x^{2}$ is larger. Then each weighted term is $X^{2} / (2 X) = X / 2$ , for samples $X > 0$ .

Ideal Importance Sampling Has Zero Variance for Nonnegative Integrands

Suppose that

f : D \to [0, \infty)

and

0 < \int_{D} f (x) d x < \infty

. If

p (x) = \frac{f (x)}{\int_{D} f (t) d t},

then each importance sampling term satisfies

\frac{f (X)}{p (X)} = \int_{D} f (t) d t .

Hence the Monte Carlo estimator has zero variance.

By the definition of

p

\frac{f (X)}{p (X)} = \frac{f (X)}{f (X) / \int_{D} f (t) d t} = \int_{D} f (t) d t,

whenever

f (X) \neq 0

. When

f (X) = 0

, this value occurs only where the density also vanishes and does not affect the estimator. Thus the estimator is constant almost surely, and a constant random variable has variance

0

Multiple Importance Sampling

Multiple importance sampling estimates an integral by combining samples drawn from several densities

p_{1}, \dots, p_{m}

. If

w_{i}

is a weight function for density

p_{i}

, then a typical weighted term has the form

w_{i} (X) \frac{f (X)}{p_{i} (X)} .

Multiple Importance Sampling Estimator

Suppose that

p_{1}, \dots, p_{m}

are sampling densities and that

X_{i, j}

is sampled from

p_{i}

, for

1 \leq j \leq n_{i}

. The multiple importance sampling estimator is

F : = \sum_{i = 1}^{m} \frac{1}{n_{i}} \sum_{j = 1}^{n_{i}} w_{i} (X_{i, j}) \frac{f (X_{i, j})}{p_{i} (X_{i, j})} .

Multiple Importance Sampling Unbiasedness Conditions

The multiple importance sampling estimator is unbiased for

\int_{D} f (x) d x

if, whenever

f (x) \neq 0

\sum_{i = 1}^{m} w_{i} (x) = 1,

and

w_{i} (x) = 0

whenever

p_{i} (x) = 0

Using linearity of expectation and the definition of density,

E (F) = \sum_{i = 1}^{m} \int_{D} w_{i} (x) \frac{f (x)}{p_{i} (x)} p_{i} (x) d x .

The condition

w_{i} (x) = 0

where

p_{i} (x) = 0

makes the terms well-defined, and the remaining factors cancel:

E (F) = \int_{D} f (x) \sum_{i = 1}^{m} w_{i} (x) d x = \int_{D} f (x) d x .

Balance Heuristic

Suppose that

p_{1}, \dots, p_{m}

are sampling densities and

n_{i}

samples are drawn from

p_{i}

. The balance heuristic assigns the weight

w_{i} (x) : = \frac{n_{i} p_{i} (x)}{\sum_{j} n_{j} p_{j} (x)} .

Balance Heuristic Weights Sum to One

\sum_{j} n_{j} p_{j} (x) \neq 0

, the balance heuristic weights satisfy

\sum_{i = 1}^{m} w_{i} (x) = 1.

Substitute the definition:

\sum_{i = 1}^{m} w_{i} (x) = \sum_{i = 1}^{m} \frac{n_{i} p_{i} (x)}{\sum_{j} n_{j} p_{j} (x)} = \frac{\sum_{i} n_{i} p_{i} (x)}{\sum_{j} n_{j} p_{j} (x)} = 1.

Power Heuristic

Suppose that

p_{1}, \dots, p_{m}

are sampling densities,

n_{i}

samples are drawn from

p_{i}

, and

β \in ℝ^{+}

. The power heuristic assigns the weight

w_{i} (x) : = \frac{{(n_{i} p_{i} (x))}^{β}}{\sum_{j} {(n_{j} p_{j} (x))}^{β}} .

MIS Compensation

Given a density

p

and

δ \geq 0

, the MIS compensated density is

p^{'} (x) : = \frac{\max (0, p (x) - δ)}{\int_{D} \max (0, p (t) - δ) d t},

whenever the denominator is positive.

Russian Roulette Estimator

Suppose that

X

is an unbiased estimator for

Q

, and let

B

be a Bernoulli random variable with

P (B = 1) = q > 0

, independent of

X

. The Russian roulette estimator is

Y : = {\begin{matrix} X / q & B = 1, \\ 0 & B = 0. \end{matrix}

Russian Roulette Estimator is Unbiased

The Russian roulette estimator is unbiased:

E (Y) = Q .

By conditioning on the two values of the Bernoulli random variable

B

E (Y) = P (B = 1) E (X / q) + P (B = 0) E (0) .

Since

P (B = 1) = q

and constants pull out of expectation,

E (Y) = q \frac{1}{q} E (X) = E (X) .

Since

X

is unbiased for

Q

E (X) = Q

Russian Roulette Estimator with a Fallback Value

Suppose that

F

is an estimator,

q \in [0,1)

c \in ℝ

, and

U \sim unif [0,1]

. The Russian roulette estimator with fallback value is

F^{'} : = {\begin{matrix} \frac{F - q c}{1 - q}, & U > q, \\ c, & U \leq q . \end{matrix}

Russian Roulette with a Fallback Value Preserves Expectation

Under the hypotheses of Russian roulette estimator with a fallback value,

E (F^{'}) = E (F) .

Since

P (U \leq q) = q

and

P (U > q) = 1 - q

E (F^{'}) = (1 - q) E (\frac{F - q c}{1 - q}) + q c .

Pulling constants out of expectation gives

E (F^{'}) = E (F) - q c + q c = E (F) .

Splitting Estimator

Suppose that

X_{1}, \dots, X_{n}

are samples from a density

p_{X}

, and for each

X_{i}

, the conditional samples

Y_{i, 1}, \dots, Y_{i, m}

are drawn from a conditional density

p_{Y | X} (y | X_{i})

. The splitting estimator for

\int_{A} \int_{B} f (x, y) d y d x

F : = \frac{1}{n} \sum_{i = 1}^{n} \frac{1}{m} \sum_{j = 1}^{m} \frac{f (X_{i}, Y_{i, j})}{p_{X} (X_{i}) p_{Y | X} (Y_{i, j} | X_{i})} .

Splitting Estimator is Unbiased

Under the hypotheses of the splitting estimator,

E (F) = \int_{A} \int_{B} f (x, y) d y d x .

By linearity of expectation,

E (F) = \frac{1}{n} \sum_{i = 1}^{n} \frac{1}{m} \sum_{j = 1}^{m} E (\frac{f (X_{i}, Y_{i, j})}{p_{X} (X_{i}) p_{Y | X} (Y_{i, j} | X_{i})}) .

For each

i, j

, use the joint density factorization

p (x, y) = p_{X} (x) p_{Y | X} (y | x)

E (\frac{f (X_{i}, Y_{i, j})}{p_{X} (X_{i}) p_{Y | X} (Y_{i, j} | X_{i})}) = \int_{A} \int_{B} \frac{f (x, y)}{p_{X} (x) p_{Y | X} (y | x)} p_{X} (x) p_{Y | X} (y | x) d y d x .

The density factors cancel, so each summand has expectation

\int_{A} \int_{B} f (x, y) d y d x

. Averaging

n m

identical expectations leaves that same integral.