1.1 Simulation of random variables

# Generating random numbers
import numpy as np
np.random.rand(20)

array([0.79871369, 0.84146444, 0.07895613, 0.87024115, 0.01566885,
       0.14597181, 0.89502404, 0.86038388, 0.09117039, 0.25184403,
       0.66376079, 0.5814277 , 0.40918816, 0.12313052, 0.5817434 ,
       0.98322438, 0.11492539, 0.07491476, 0.29077235, 0.77603614])

1.2 Generating discrete distributions

1.3 Examples

• Roll a dice

# Generating the roll of n dice
import numpy as np
n=10
np.ceil(6*np.random.rand(n))

array([3., 3., 4., 2., 1., 3., 5., 4., 4., 2.])

1.4 Inversion method

1.5 Rejection method

The next method permits to simulate X provided we can simulate Y (X and Y should have the same support).

1.6 Example

Suppose X has pdf f(x)=20x(1-x)^3 for 0 \le x\le 1. Since X is supported on [0,1] we pick Y uniform on [0,1] with g(y)=1. Then C = \max \frac{f(x)}{g(x)} = \max_{x\in [0,1]} 20x(1-x)^3 = \frac{135}{64} So generate two random numbers U_1 and U_2 and if U_1 \le \frac{256}{27}U_2 (1-U_2)^3 set X=U_2. The average proportion of accepted values is \frac{64}{135}=.4747..

import numpy as np
import matplotlib.pyplot as plt
def accept_reject(N):   # Generate N samples
    n_accept=0
    x_list = [] 
    while n_accept < N:
        a=np.random.rand(2)
        if a[0] < (256/27)*a[1]*(1-a[1])**3 :
            n_accept += 1
            x_list.append(a[1])
    return x_list
plt.figure(figsize=(5,3))  
plt.hist(accept_reject(100000), bins=100, density= 'true')
t = np.arange(0., 1., 0.02)
plt.plot(t, 20*t*(1-t)**3, 'r--' )
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Histogram vs the exact pdf')
plt.show()  

1.7 Normal RV with the rejection method

• We use the fact that we know how to generate an exponential random variable (we will take \lambda=1)

1.8 Box-Muller algorithm

1.9 Generating a random vector on the unit ball

1.10 Generating random permutations

2.1 Law of large numbers

We start by recalling some basic limit theorems from probability theory.

• X_1, X_2, \cdots independent and identically distributed random variables (IID). Think of repeating a random experiment n times and recording the results. We denote by S_n = X_1 + X_2 + \cdots + X_n the sum.

• The Law of Large numbers shows that if X_i have a finite mean \mu then the empirical average \frac{S_n}{n} converges to \mu.

2.2 The Monte-Carlo method

• The (simple) Monte-Carlo method is a probabilistic algorithm using sums of independent random variables the law of large numbers to estimate a (deterministic) quantity \mu \in \mathbb{R} (or \mathbb{R}^d).
  The basic idea is to express \mu as the expectation of some random variable \mu = E [ h(X) ] and then use the law of large numbers to build up an estimator for \mu.

• An interesting part is that there are, in general, many ways to find the random variables h(X) and we will play with this, e.g. for importance sampling.

• Conversely in many problems the random variable h(X) is given but the expection is too diffcult to compute, so we rely on the LLN to compute \mu.

2.3 Buffons needles

2.4 Computing \pi

• Use the rejection method: Enclose a disk of radius 1 in a square of side length 2 and consider the following Bernoulli random variable X.

  • Generate 2 independent vectors V_1,V_2 uniformly distributed on [-1,1].

  • If V_1^2 + V_2^2 \le 1, set X=1, otherwise set X=0.

X is Bernoulli with p=\frac{\textrm{area of the disk}}{\textrm{ area of the square}}=\frac{\pi}{4} so \frac{4}{N}\sum_{k=1}^N X_i \to \pi with probability 1

import numpy as np
N = 1000000

dart_positions = 2 * np.random.rand(N, 2) - 1  # generates numbers in [-1, 1]
Ncircle = [0] # start the count with 0 to make our loop-logic easier

for x,y in dart_positions:
    if  x**2 + y**2 < 1:
        Ncircle.append(Ncircle[-1] + 1) # another dart has fallen in the circle
    else:
        Ncircle.append(Ncircle[-1])  # the dart fell outside of the circle - Ncircle is unchanged

running_estimate = []
for n_total, n_circle in enumerate(Ncircle[1:]): # skip the inital 0-count
    # n_total starts at 0, so we need to add 1
    running_estimate.append(4 * n_circle / (n_total + 1))

np.transpose(running_estimate[-20:])

array([3.14173569, 3.14173655, 3.14173741, 3.14173827, 3.14173913,
       3.14173998, 3.14174084, 3.1417417 , 3.14174256, 3.14173942,
       3.14174028, 3.14174113, 3.14173799, 3.14173885, 3.14173971,
       3.14174057, 3.14174143, 3.14174228, 3.14174314, 3.141744  ])

2.5 Monte-Carlo to compute integrals

• Goal is to computes integral: \displaystyle I\,=\, \int_{S} h( {\bf x}) d {\bf x} for some function h:S\subset \mathbb{R}^d \to R.

2.6 Goals

• Performance gurantees: How many sample should we generate to obtain a given precision for the computation of \mu = E[h(X)]

  • Use central limit theorems type results to obtain asymptotic confidence intervals.

    • Simple and useful but somewhat unsatisfying.

  • Use concentration inequalities to obtain non-asymptotic confidence intervals.

    • True performance guarantess but harder to get and often too pessimistic.

2.7 Central limit theorem

• IID random variables X_1, X_2, \cdots *independent with \mu=E[X_i] and \sigma^2 = \textrm{Var}(X_i).

2.8 Confidence intervals (version 1)

• We build build confidence interval for a Monte-Carlo estimator I_N=\frac{1}{N}\sum_{k=1}^N h(X_k), since by the Central limit theorem \frac{\sqrt{N}}{\sigma} ( I_N - \mu) is asymptotically normal.

2.9 Central limit theorem with Slutsky theorem

• Building an estimator for the variance

2.10 Confidence intervals (version 2)

• We can get an asymptotic confidence interval using Theorem 2.5

P\left( \mu \in \left[ I_N - z_\alpha \frac{V_N}{\sqrt{N}}\,,\, I_N + z_\alpha \frac{V_N}{\sqrt{N}} \right] \right) \preceq \alpha \,. or

• Take home message.

  • Theory: build the Monte-Carlo estimator with the smallest possible variance long\rightarrow Variance Reduction

  • For a precision \varepsilon ne needs N \sim \frac{1}{\varepsilon^2} samples.

  • Practice: With your code always compute both I_N and V_N!

2.11 Examples

• Compute \pi with accept/reject. Estimator I_N=\frac{4}{N} \sum_{k=1}^N {\bf 1}_{\{V_1^2 + V_2^2 \le 1\}} where V_i, i=1,2 is uniform on [-1,1] so we get \textrm{Var}(I_N) \,=\, \frac{16}{N}\frac{\pi}{4}\left(1-\frac{\pi}{4}\right)= \pi(4-\pi) = \frac{1}{N} 2.6967 To get \pi with a precision of \varepsilon=10^{-3} and 99% confidence we need at least N = 2.576 \times 2.6967 \times 10^6= 6.946,699 samples.

• Let h(x)= \sqrt{1 - x^2} for x \in [0,1]. Then \int_0^1 h(x) dx = \frac{\pi}{4} and a estimator for \pi is I_N=\frac{16}{N} \sum_{k=1}^N h(U_k) with variance \textrm{Var}(I_N) \,=\, \frac{16}{N} \left( E[h(U)^2] -E[h(U)]^2 \right)= \frac{1}{N}\left(\frac{2}{3}-\left(\frac{\pi}{2}\right)^2\right)= \frac{1}{N} 0.797 Better, roughly 4 times less sample (plus using only half as many random numbers).

3.1 Antithetic random variables and common random numbers

• The idea is to use dependent random variables or dependent random numbers to decrease the variance. Recall the formula for variance: for random variables X and Y we have \textrm{Var}(X+Y) = \textrm{Var}(X)+ \textrm{Var}(Y)+ 2 \textrm{Cov}(X,Y) \,. If X and Y are negatively correlated, i.e. \textrm{Cov}(X,Y) < 0 then we can possibly reduce the variance.

3.2 A simple example

• Compute \mu = E\left[ e^U \right] use MC (of course the answer is e-1…).

• Since h(x)=e^x is increasing and g(x)=e^{1-x} is decreasing we can use antithetic RV.

• We have the covariance \textrm{Cov}\left(e^U, e^{1-U}\right)=E\left[e^U e^{1-U}\right]-E\left[e^U\right] E\left[e^{1-U}\right]=e -(e-1)^2=-0.2342

• We also have \textrm{Var}\left(e^U\right)= E\left[e^{2U}\right]- E\left[e^U\right]^2 = \frac{e^{2}-1}{1}-(e-1)^2=0.2420

• To compare:

  • Independent random numbers U_1,U_2 give the variance \textrm{Var}\left( \frac{e^{U_1}+e^{U_2}}{2} \right) = \frac{ \textrm{Var}\left(e^{U}\right)}{2}= 0.1210

  • Antithetic random variables U, 1-U give the variance \textrm{Var}\left( \frac{e^{U_1}+e^{U_2}}{2} \right) = \frac{ \textrm{Var}\left(e^{U}\right)}{2} + \frac{\textrm{Cov}\left(e^U, e^{1-U}\right) }{2}=.0039

• A variance reduction of 96.7 percent (not bad…)

3.3 Reliability function

• A system consists of n components, each of which is either works or fails. We set x_i= \left\{ \begin{array}{cl} 1 & \textrm{ if component } i \textrm { works} \\ 0 & \textrm{ if component } i \textrm { fails} \end{array} \right. \,. We call \mathbf{x}=(x_1, \cdots, x_n) the state vector describing the state of the system.

3.4 Stucture function

General principled way to compute the structure function

• A path A is a subset A \subset \{1,\cdots, n\} such that h(\mathbf{x})=1 if x_i=1 for all i\in A, that is it a subset of components which will ensure that the system works.

• A minimal path A is a path A such that no subset of A is itself a path, that is it a minimal subset of components which will ensure that the system works.

The theorem is of limited usability since it requires enumeration of all minimal paths. If the system is moderately big, this is very difficult and/or tedious.

3.5 Monte-Carlo for structure functions

• If we assume a probabilistic model, which is pretty reasonable. Think for example, of a routing network, where each node have a certain probability of failure.

• The reliability function g(U_1, \cdots, U_n)= h\left( \mathbf{1}_{\{U_1\le p_1\}}, \cdots, \mathbf{1}_{\{U_n\le p_n\}} \right) is an increasing function so \textrm{Cov}(g(U_1, \cdots, U_n),g(1-U_1, \cdots, 1-U_n) ) <0 and we can reuse each random number and reduce the variance.

3.6 Importance sampling

• Consider continuous RVs but the same ideas works for discrete RVs. Suppose X has pdf f(x) and we want to compute \mu= E[h(X)] = \int h(x) f(x) dx. One can always use the simple MC algorithm I_N =\frac{1}{N} \sum_{k=1}^N h(X_k).

3.7 Examples

• Suppose we want to sample the mean of standard normal RV X using an importance sampling estimator (it may look like a trivial example but wait…).
  So we take h(x)=x and let us pick Y to be a normal RV with \mu=0 and variance \sigma^2 E\left[h(x)^2 \frac{f(x)}{g(x)}\right] = \int_{-\infty}^\infty x^2 \frac{ (e^{-x^2/2}/\sqrt{2 \pi})^2}{e^{-x^2/2\sigma}/\sqrt{2 \pi} \sigma } \, dx = \sigma \int_{-\infty}^\infty \frac{1}{\sqrt{2 \pi}} x^2e^{-x^2( 2 - \sigma^{-2})/2} \, dx = \left\{ \begin{array}{cl} \frac{\sigma}{(2-\sigma^{-2})^{3/2}} & \textrm{ if } \sigma^2 > 1/2 \\ +\infty & \textrm{otherwise} \end{array} \right. An elementary computation shows that the optimal choice is \sigma^2 = 5/2.

3.8 Tilted distribution

• Tilted distributions are a natural tool for importance sampling. If f(x) is the pdf of X consider the random variable X^{(t)} with the new distribution f_t(x) = \frac{e^{tx}}{M(t)}f(x) \quad \quad \textrm{ where } M(t)=E[e^{tX}]=\int e^{tx} f(x) dx that is M(t) is the moment generating function of the random variable X. We have M'(t)=E[Xe^{tX}], M''(t)=E[X^2e^{tX}], \cdots

3.9 Rare event simulation

• Suppose X_1, \cdots X_n are IID random variables with joint densities f(x_1)f(x_2)\cdots f(x_n) and S_n=X_1+\cdots+ X_n the sum. We want to estimate the (high-dimensional integral) P\left( \frac{S_n}{n} \ge a\right) = E\left[1_{\{S_n \ge n a\}}\right] = \int 1_{\{S_n \ge na\}} f(x_1)f(x_2)\cdots f(x_n) dx_1 \cdots dx_n If a > E[X_1] this probability is very small (exponentially small in n, see Chapter 4) and thus the simple Monte-Carlo will fail miserably.

3.10 A reliability example

4.1 Chernov bounds

We start with a basic but fundamental bound

4.2 Examples of Chernov bounds

• If X is normal with mean \mu and varaince \sigma^2, then M(t)=e^{\mu t + \frac{\sigma^2 t^2}{2}} and, for a = \mu + \varepsilon, Chernov gives P(X \ge \mu + \varepsilon) \le \inf_{t \ge 0} \frac{e^{\mu t + \frac{\sigma^2 t^2}{2}} }{e^{t(\mu + \varepsilon)}} = e^{ -\sup_{t\ge 0}\left\{ t\varepsilon - \frac{\sigma^2 t^2}{2}\right\} } = e^{-\frac{\varepsilon^2}{2\sigma^2}} This is sharp: one can show (see exercises) that \left(\frac{1}{\varepsilon}-\frac{1}{\varepsilon^3\sigma^2} \right)\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{\varepsilon^2}{2\sigma^2}} \le P(X \ge \mu + \varepsilon) \le \frac{1}{\varepsilon}\frac{1} {\sqrt{2\pi}\sigma} e^{-\frac{\varepsilon^2}{2\sigma^2}}

4.3 Hoeffding’s bound

• Chernov bounds are very sharp but requires knowledge of the mgf M_X(t).

4.4 Hoeffding’s theorem

4.5 McDiarmid theorem

A very useful extension of Hoeffding’s applies to suitable function h(X_1, X_2, \cdots, X_n) other than just the sum.

We say that h(x_1, x_2, \cdots, x_n) satisifies the bounded difference property if there exist constants c_k such that for all x_k, x_k' |h(x_1,\cdots, x_k, \cdots, x_n) - h(x_1,\cdots, x_k', \cdots, x_n)| \le c_k

4.6 Confidence interval for variance estimator

We show how to build confidence interval for the unbiased variance estimator V_N(X_1, \cdots, X_N)=\frac{1}{N-1} \sum_{k=1}^N\left(X_i - \frac{S_N}{N}\right)^2 = \frac{1}{N-1}\sum_{i=1}^N X_i^2 - \frac{S_N^2}{(N-1)N} If we change X_1 into \widehat{X}_1 then S_N(X_1, \cdots, X_N ) - S_N(\widehat{X}_1, \cdots, X_N ) = X_1 - \widehat{X}_1 and so V_N(X_1, \cdots, X_N)- V_N(\widehat{X}_1, \cdots, X_N) = \frac{X_1^2 - \widehat{X}_1^2}{N-1} - \frac{X_1 - \widehat{X}_1}{N-1} \left(\frac{S_N(X_1, \cdots, X_N )}{N} + \frac{S_N(\widehat{X}_1, \cdots, X_N )}{N} \right) So if we assume |X_i|\le c we have the bound

|V_N(X_1, \cdots, X_N)- V_N(\widehat{X}_1, \cdots, X_N)| \le \frac{5c^2}{N-1} and thus by McDiarmid we get P( V_N - \sigma^2 \ge \varepsilon ) \, \le \, e^{ - \frac{(N-1)^2}{N}\frac{\varepsilon^2}{50c^4}} \,\le\, e^{ -(N-2)\frac{\varepsilon^2}{50c^4}} and this decay exponentially in N again and can be used for a confidence interval for the variance which is very similar to the one for the mean.

4.7 Bernstein bound

In Hoeffding’s bound we use, in an essential way, a bound on the variance. If the variance is small then one should expect the bound to be poor. Bernstein bound, which can be used if we have some a-priori knowledge about the variance, improves.

4.8 Bernstein confidence interval

To obtain a confidence interval we need to solve an optimization problem: after some straightforward computation we find

\sup_{t \ge 0} \left\{ \varepsilon t - \frac{a t^2}{2(1-bt)} \right\} = \frac{a}{b^2}h\left( \frac{b\epsilon}{a}\right) \quad \textrm{ where } h(u) = 1 +u - \sqrt{1 + 2u} and we note that the h^{-1}(z)= z + \sqrt{2z}. We obtain the exact same bound for the left tail (by symmetry) and thus
to obtain a confidence interval we have \delta = 2 e^{-\frac{a}{b^2}h\left( \frac{b\varepsilon}{a}\right)} \iff \varepsilon = b \ln\left( \frac{2}{\delta}\right) + \sqrt{2a \ln\left( \frac{2}{\delta}\right)}

(\textrm{Bernstein}) \quad \frac{c}{3N} \ln\left( \frac{2}{\delta}\right) + \sqrt{2 \ln\left( \frac{2}{\delta}\right)} \frac{\sigma}{\sqrt{N}} \quad \textrm{ versus} \quad \sqrt{2\ln\left(\frac{2}{\delta}\right)} \frac{\sigma_{max}}{\sqrt{N}} \quad (\textrm{Hoeffding})

You should note that this bound can be substantially better than Hoeffding’s bound is \sigma \le \sigma_{max}

4.9 Comparing the bounds for Bernoulli

Hoeffdings amount to bounding the variance by \frac14. Bernstein which incorporporates knowledge of (or a bound on) the variance is much better than Hoeffdings as soon as the variance is not too big (0.1 in the epicture below.)