Lab 3: Sampling and Monte Carlo Integration

Lab 3: Sampling and Monte Carlo Integration#

This lab is adapted from AM207 (Harvard). Data files are in the book’s data/ folder. Worked solutions are collapsed: try each exercise before expanding them.

In this lab

Monte Carlo Estimation of a General Integral
Practice: Implementing Monte Carlo Integration

Arbitrary Definite Integrals#

We formulated Monte Carlo estimation as a way to approximate the expected value of functions of random variables. However, given a function of deterministic variables, \(f(x)\), we can rewrite the integral

\[ \int_{\Omega} f(x) dx, \]

where \(x\in\mathbb{R}^D\) and \(\Omega \subseteq \mathbb{R}^D\), like an expectation

\[ \int_{\Omega} f(x) dx = \int_{\Omega} g(x)h(x) dx, \]

if we can factor \(f\) as \(f(x) = g(x)h(x)\).

The rewritten integral looks a lot like the expectation of \(g(X)\) where the pdf of \(X\) is \(h(X)\). The only problem is that \(h(x)\) is not necessarily a pdf!

Rewriting an Arbitrary Definite Integral as an Expectation#

We want to rewrite the following expression as an expectation:

\[ \int_{\Omega} f(x) dx = \int_{\Omega} g(x)h(x) dx. \]

Let \(V\) denote the integral of \(h\) over \(\Omega\):

\[ V = \int_{\Omega} h(x)dx. \]

Note that, then, \(\int_{\Omega} \frac{1}{V}h(x)dx = 1\). Hence, \(\frac{1}{V}h(x)\) is a valid pdf over \(\Omega\).

Then we can write:

\[ \int_{\Omega} f(x) dx = \int_{\Omega} \left(Vg(x)\right)\left(\frac{1}{V}h(x)\right)dx. \]

We recognize the above integral as an expectation, and hence can estimated using Monte Carlo methods:

\[ \mathbb{E}_{X\sim \frac{1}{V}h(X)}[Vg(X)] \approx V\sum_{s=1}^S g(X_s),\; X_s\sim \frac{1}{V}h(X). \]

How Do You Factor the Integrand \(f\)?#

Once we factor \(f(x)\) as \(f(x) = g(x)h(x)\), then we can rewrite the integral of \(f\) (over \(\Omega\)) as an expectation

\[ \int_{\Omega} f(x) dx = \int_{\Omega} \left(Vg(x)\right)\left(\frac{1}{V}h(x)\right)dx = \mathbb{E}_{X\sim \frac{1}{V}h(X)}[Vg(X)]. \]

But how does one choose \(g\) and \(h\)? Two simple things to keep in mind:

since we need to sample from \(h\), we need to choose \(h\) such that it is easily approximated by distributions from which we know how to sample.
since the variance of the Monte Carlo estimator depends on \(g\), we should choose \(g\) to be relatively flat (so that the output of \(g\) has low variance).

Exercise: Monte Carlo Integration#

Use Monte Carlo methods to estimate the following integral:

\[ \int_{0}^{1} \frac{3x^2}{x^3 - x^2 - 11x - 4} dx \]

Hint: write the integrand \(f\) as \(f(x) = g(x)h(x)\) where \(h(x)\) defines a pdf over \([0, 1]\) and is easy to sample from.

Solution

The integrand \(f(x) = \frac{3x^2}{x^3 - x^2 - 11x -4}\) naturally decomposes as

\[ f(x) = \frac{1}{x^3 - x^2 - 11x -4}\cdot 3x^2 \]

Out of the two factors, it looks like \(3x^2\) is easier to integrate (and differentiate) so it’s a good candidate for a pdf. Thus, we set

\[ h(x) = 3x^2,\quad g(x) = \frac{1}{x^3 - x^2 - 11x -4}. \]

Before writing a sampler for \(h(x)\), we need to check that it defines a pdf on \([0, 1]\)! To this end, we compute:

\[ \int_0^1 h(x) dx = \left.x^3\right\vert_{0}^{1} = 1 \]

So we see that \(h(x)\), when restricted to \([0, 1]\) defines a pdf!

But how do we sample from \(h(x)\) (after all, \(h(x)\) doesn’t look like the pdf of a well-known distribution)? Since the CDF \(H(a) = \int_{0}^a h(x) dx = a^3\) is an invertible function, \(h\) is a natural candidate for inverse CDF sampling. The inverse of the CDF of \(h\) is

\[ H^{-1}(a) = \sqrt[3]{a}. \]

Now, we write rewrite the integral we want to compute as an expectation

\[ I = \int_{0}^{1} \frac{3x^2}{x^3 - x^2 - 11x - 4} dx = \mathbb{E}_{X\sim h(X)}\left[ g(x)\right]. \]

The naive Monte Carlo estimator for \(I\) is

\[ \widehat{I} = \frac{1}{S}\sum_{s=1}^S \frac{1}{X_s^3 - X_s^2 - 11X_s -4},\; X_s\sim h(X). \]

We implement this estimator below:

#implementation of inverse CDF sampling
def inverse_cdf_sampling(CDF_inv, S):
    U = np.random.uniform(0, 1, size=S)
    samples = CDF_inv(U)
    return samples

#implementation of naive monte carlo
def monte_carlo(sampler, g, S, N=1):
    estimates = []
    for n in range(N):
        #sample from sampler
        samples = sampler(S)
        #evaluate g at each sample and average
        estimates.append(g(samples).mean())
    return estimates
    
#definition of h
h = lambda x: 3 * x**2
#definition of the inverse CDF of h
H_inv = lambda x: x**(1./3.)
#definition of g
g = lambda x: 1. / (x**3 - x**2 - 11 * x - 4.)
#define our sampler for the estimator
sampler = partial(inverse_cdf_sampling, H_inv)

#total number of samples
S = 100
#total number of trials
N = 500

#perform monte carlo estimate N number of times
estimates = monte_carlo(sampler, g, S, N=N)

print('mean of {} number of monte carlo estimates: {}'.format(N, np.mean(estimates)))
print('variance of {} number of monte carlo estimates: {}'.format(N, np.std(estimates)**2))

mean of 500 number of monte carlo estimates: -0.0840213302768032
variance of 500 number of monte carlo estimates: 3.4118616765217585e-06

The sample variance of the naive Monte Carlo is already small, but we can still consider some variance reduction techniques. In particular, we can consider importance sampling.

Importance Sampling

If we want to compute the importance sampling estimate we need an importance distribution for \(h\). Since \(h\) is a distribution on \([0, 1]\), it’s natural to consider the beta distribution as an importance distribution.

#definte the beta pdf
beta_pdf = sp.stats.beta(5, 1).pdf

#visualize both the target and proposal distributions
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
x = np.linspace(0, 1, 100)
ax.plot(x, h(x), color='black', label='h(x)')
ax.plot(x, beta_pdf(x), color='red', label='beta pdf, alpha={}, beta={}'.format(5, 1))
ax.set_title('h(x) vs the importance distribution')
ax.legend(loc='best')
plt.show()

Implementing the importance sampling estimate is straight-forward:

\[ I \approx \widehat{I}_{\text{IS}} = \frac{1}{S}\sum_{s=1}^S \frac{h(X_s)}{q(X_s)} g(X_s),\; X_s \sim Beta(5, 1) \]

We just need to sample from the Beta, compute the importance weights \(\frac{h(X_s)}{q(X_s)}\) and compute \(g(X_s)\).

#implement importance sampling MC
def importance_sampling_MC(sampler, g, h, S, N=1):
    estimates = []
    weights_over_trials = []
    for n in range(N):
        #sample from sampler
        samples = sampler(S=S)
        #evaluate the importance weights
        weights = np.multiply(h(samples), sampler(x=samples)**-1)
        #multiply weights to g(x) and average
        estimate = (weights * g(samples)).mean()

        weights_over_trials.append(weights)
        estimates.append(estimate)
        
    return estimates, np.array(weights_over_trials)

#defines a function that can either sample from a beta distribution or get it's pdf values
def beta_sampler(S=None, x=None):
    if S is not None:
        samples = np.random.beta(5, 1, size=S)
        return samples
    elif x is not None:
        likelihood = beta_pdf(x)
        return likelihood
    else:
        return None
    
    
#perform monte carlo estimate N number of times
estimates, weights = importance_sampling_MC(beta_sampler, g, h, S, N=N)

print('mean of {} number of importance sampling monte carlo estimates: {}'.format(N, np.mean(estimates)))
print('variance of {} number of importance sampling monte carlo estimates: {}'.format(N, np.std(estimates)**2))

mean of 500 number of importance sampling monte carlo estimates: -0.08496357680015376
variance of 500 number of importance sampling monte carlo estimates: 0.00018620123953831943

The variance of importance sampling Monte Carlo is higher than that of the naive Monte Carlo! Weren’t we promised that importance sampling was a variance reduction method? Looking at the discussion of the design of \(q\) (“Nitty Gritty of Importance Sampling”) in Chapter 7, we realize that an important criteria was that:

the importance distribution \(q\) look like the nominal distribution \(h\), otherwise the importance weights become heterogenous (polarized)
the importance distribution is roughly proportional to \(h|g|\).

Let’s plot both the weights from one of the trials and \(h|g|\) to diagnose the problem.

n = 20
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.hist(weights[n], bins=30, color='red', alpha=0.4)
ax.set_title('histogram of importance weights from trial N={}'.format(20))
plt.show()

We see that the majority of the weights are below 1.0, so these samples have a “small effect” on the overall estimate. We might conjecture that maybe this is why the variance of \(\widehat{I}_{\text{IS}}\) is higher (try confirming this conjecture by computing the effective sample size!).

fig, ax = plt.subplots(1, 2, figsize=(20, 5))
x = np.linspace(0, 1, 100)
ax[0].plot(x, h(x) * np.abs(g(x)), color='black')
ax[0].set_title('h|g|')

ax[1].plot(x, beta_pdf(x), color='red', label='beta pdf, alpha={}, beta={}'.format(5, 1))
ax[1].set_title('importance distribution')
plt.show()

From this visualization we see that our importance distribution doesn’t look very much like the product \(h|g|\). Maybe this is the source of the higher variance of our importance sampling Monte Carlo estimate.

To test our hypothesis, you can try to pick an importance distribution that looks more proportional to \(h|g|\) and see if using that distribution will decrease the variance of \(\widehat{I}_{\text{IS}}\).