• Gain an intuitive understanding of how a process operates

• Distribution or properties of distribution often analytically unavailable

• e.g. we've already seen that the posterior distribution \(f(\theta|D)\) involves an integral which is, typically, impossible to calculate exactly.

• Even when \(f(\theta|D)\) is available, finding some property of it such as an expectation, involves another large integral.

• Integration is hard!

Basic idea is obtain samples of values \(\{\theta^{(1)}, \theta^{(2)}, \theta^{(3)},\ldots,\theta^{(n)} \}\) from distribution of interest and use this sample to estimate properties of distribution.

• Example: The mean of the distribution \(\pi(\theta)\) is \[E[\theta] = \int_{\theta \in \Omega} \theta \pi(\theta) d\theta.\] This can be approximated by \[E[\theta] \approx \bar{\theta} = \frac1 n \sum_1^n \theta^{(i)}\] where \(\theta_i \sim \pi(\theta)\).

Say \(\theta\) is a univariate random variable.

• Suppose we can simulate values of \(\theta\)…
• but cannot find the distribution function analytically.
• How could we study it?

• Let's draw 1000 samples, \(\{\theta^{(1)}, \theta^{(2)}, \theta^{(3)},\ldots,\theta^{(1000)} \}\).

• The first 20 look like 4.359, 4.42, 5.323, 2.814, 1.476, 7.14, 2.255, 3.211, 3.085, 2.31, 4.549, 3.078, 2.14, 0.802, 4.558, 8.176, 6.257, 2.385, 1.263, 1.455 etc.

Now make a histogram to get an idea of the shape of the distribution:

Let's approximate the mean of the distribution by the average of the sampled values:

\(E[\theta] \approx \bar{\theta} = \frac1 {1000} \sum_1^{1000} \theta^{(i)} =\) 3.8390958

• The sample \(\{\theta^{(1)}, \theta^{(2)}, \theta^{(3)},\ldots,\theta^{(n)} \}\) is random so estimates based on it will vary with different samples.

• As the number of samples \(n\) increases, the approximation gets better.

• In general, if want to estimate \(E[g(\theta)]\) by \(\bar g (\theta) = \frac1 n \sum_1^n g(\theta^{(i)})\), \(\bar g\) is normally distributed with mean \(E[g]\) and variance \(\mbox{var}(\bar g) = \mbox{var}(g)/ n\).

• When stated formally, this is known as a central limit theorem (CLT).

• So with lots of samples, arbitrarily good approximations to quantities of interest can be made.

• Interactive example of the CLT here

• Simulation relies on a ready supply of random numbers
• Computers are, by design, determinisitc
• So instead of using truly random numbers, use pseudo-random numbers
• Indistinguishable from real random numbers for most purposes
• Since sequences of pseudo-random numbers are repeatable, often better than truly random numbers for scientific simulations
• Generate pseudo-random numbers that are uniformly distributed between 0 and 1, i.e. \(u \sim U(0,1)\)
• Then transform \(u\) to have the desired distribution.

• have a very long period so that no repetitions of the cycle occur
• be efficient in memory and speed
• repeatable so that simulations can be reproduced
• portable between implementations
• have efficient jump-ahead methods so that a number of virtual streams of random numbers can be created from a single stream for use in parallel computations; and,
• have good statistical properties approximating the uniform distribution on \([0,1]\).

We omit details of how RNGs actually work (see notes for short description)

Important: Standard random number generators used as the default in many languages (e.g. java.util.Random) are not adequate for scientific simulation studies.

I recommend instead the Mersenne Twister which is implemented in most languages. It is the default generator in Python.

In Python to repeart a sequence of random bits, choose a seed a and call random.seed(a) to set generator in same state

To sample from \(X\) a discrete r.v. with \(P(X = x_i) = p_i\) 1. generate \(u \sim U(0,1)\) 2. if \(u< p_1\), set \(X = x_1\) otherwise if \[\sum_{j = 1}^{i-1}p_j < u < \sum _{j = 1}^i p_j\] set \(X = x_i\).

Suppose \(X\) can take values A, B, C, D or E with the probabilities

##   value probability
## 1     A        0.10
## 2     B        0.30
## 3     C        0.15
## 4     D        0.35
## 5     E        0.10

Suppose we draw \(u \sim U(0,1)\) and get \(u =\) 0.503.

Then \(x\) takes the value C since \[\Pr(A) + \Pr(B) = 0.1+0.3 = 0.4 \leq u = 0.503 \leq 0.55 = \Pr(A) + \Pr(B) + \Pr(C)\].

##   value probability
## 1     A        0.10
## 2     B        0.30
## 3     C        0.15
## 4     D        0.35
## 5     E        0.10

Drawing a uniform rv of 0.503 produces a C.

So generate 20 realisations of the rv \(X\) by generating the uniform rv \(U\) and transforming: for u

##  [1] 0.601 0.223 0.376 0.048 0.474 0.274 0.197 0.765 0.197 0.435 0.103
## [12] 0.260 0.691 0.794 0.800 0.808 0.528 0.375 0.740 0.250

Produces \(x\): D, B, B, A, C, B, B, D, B, C, B, B, D, D, D, D, C, B, D, B

Suppose we want to simulate from \(X\) which has the cumulative distribution function (cdf) \(F(X)\).

Remember, if \(X\) has probability distribution function \(f(X)\) then the cdf is \(F(X) = \int_{-\infty}^x f(t) \, dt\)

Result (Inversion method)

If \(U \sim U(0,1)\), then \(X = F^{-1}(U)\) produces a draw from \(X\) where \(F^{-1}\) is the inverse of \(X\).

So when we can find the inverse of the cdf is known, sampling from the distribution is easy.

If \(X \sim \mbox{Exp}(\lambda)\) then \(F(x) = 1-e^{-\lambda x}\) and so \(F^{-1}(u) = - \frac{\log(1-u)}{\lambda}\)

So generate exponential rv \(x\) by generating the uniform rv \(u\) and set \(x = - \log(1-u)/\lambda\)

So to draw 20 exponential rvs, draw 20 uniform rvs and transform them:

\(u\): 0.012, 0.057, 0.187, 0.238, 0.251, 0.295, 0.298, 0.351, 0.397, 0.464, 0.473, 0.479, 0.485, 0.536, 0.602, 0.646, 0.721, 0.761, 0.915, 0.938

\(x\): 0.01, 0.06, 0.21, 0.27, 0.29, 0.35, 0.35, 0.43, 0.51, 0.62, 0.64, 0.65, 0.66, 0.77, 0.92, 1.04, 1.28, 1.43, 2.47, 2.78