Simulation of Random Variables
Attribution
Much of this material is from Chapter 2 of Robert, Casella, and Casella (1999) and Chapter 4 of Owen (2013)
Recap
From last week we had two important theorems that will aid in computing intractable integrals and in assessing our computational precision
Weak Law of Large Numbers
This says that for \(X_i \overset{\mathrm{iid}}{\sim} F, i = 1, \dots, S\), and \(\hat{\mu}_S = \frac{1}{S}\sum_{i=1}^S X_i\), for all \(\epsilon > 0\):
\[ \lim_{S\to\infty} P(\abs{\hat{\mu}_S - \Exp{X}} < \epsilon) = 1 \]
The shorthand for this is \(\hat{\mu}_S \overset{p}{\to} \Exp{X}\)
Central limit theorem
This says that the following convergence in distribution holds, for \(\Var{X_i} = \sigma^2\):
\[ \sqrt{S}(\hat{\mu}_S - \Exp{X}) \overset{d}{\to} N(0, \sigma^2) \] This is often informally stated as the following:
\[ \hat{\mu}_S \overset{\mathrm{asy}}{\sim} N(\Exp{X}, \frac{\sigma^2}{S}) \]
Motivation for today’s lecture
The preceding theorems suggest a strategy to solve intractable integrals, assuming that \(F(x)\) admits a density with respect to the Lebesgue measure, \(f(x) dx\): \[ \Exp{h(X)} = \int h(x) f(x) dx \]
- Generate \(S\) iid draws, \(X_i\) from \(F\)
- Set \(\hat{\mu}_S = \frac{1}{S} \sum_{i=1}^S h(X_i)\)
- As \(S\to\infty\), \(\hat{\mu}_S \overset{p}{\to} \Exp{h(X)}\)
Now, that’s all well and good, but how do we generate iid draws of \(X_i\)? This is what today’s lecture is about.
Pseudo-random number generators
We’ll take for granted that we can generate a suitably random sequence of independent uniformly-\([0,1]\)-distributed random variables. “Suitably random independent” RVs means that given the draws \(U_1, \dots, U_n\) we have no way of predicting \(U_{n+1}\).
We call them “pseudo-random numbers” because, given \(u_0\) and a transformation function \(D\), the sequence of draws given by \((u_i) = (D^{i}(u_0))\) is exactly reproducible given the initial value (also known as the seed), and \(D\). They’re uniform because the set of values \(u_1, \dots, u_n\) is statistically indistinguishable from the set of independent uniformly distributed \(V_1, \dots, V_n\). As with anything, whether we can distinguish between psuedo random numbers and the theoretical distribution depends on what statistical tests we use.
Uniformly distributed random variables are important tools in generating non-uniformly distributed random numbers. Consider the task of generating a Bernoulli RV \(X\) with probability \(\theta\). With a uniform \(U\), \(X = \ind{U \leq \theta} = \ind{1 - U > 1 - \theta} = \ind{U > 1 - \theta}\).
This is actually a special case of a more general algorithm, the inverse CDF method. Recall the definition of the CDF, \(F(x)\) for a random variable \(X\): \[ F(x) = P(X \leq x) \] First let’s do the continuous case and we’ll come back to the discrete case. In the continuous case, we have that \(X = F^{-1}(U)\): \[ \begin{aligned} P(X \leq x) & = P(F^{-1}(U) \leq x) \\ & = P(F(F^{-1}(U)) \leq F(x)) \\ & = P(U \leq F(x)) \\ & = F(x) \end{aligned} \]
Thus, \(X\) is distributed according the RV with CDF \(F(x)\).
Example 1 (Exponential from uniform) We want to generate \(X \sim \text{Exp}(1)\) from \(U \sim \text{Unif}(0,1)\).
The CDF of the exponetial is \[ P(X \leq x) = 1 - e^{-x} \] Inverting this function gives: \[ x = -\log(1 - u) \]
So for \(U \sim \text{Unif}(0,1)\):
\[ X = -\log(1 - U) \implies X \sim \text{Exp}(1) \]
The discrete case is a little more tricky because \(F^{-1}(u)\) isn’t unique. In the Bernoulli example, we have: \[ F(x) = \begin{cases} 0 & x < 0 \\ 1 - \theta & 0 \leq x < 1 \\ 1 & x \geq 1 \\ \end{cases} \]
The inverse of this function isn’t unique, because any \(x\) between \(0\) and \(1\) (inclusive of \(0\) and exclusive of \(1\)) will give \(1-\theta\). In order to deal with this we’ll define the following function which is unique: \[ F^{-}(u) = \inf(x \mid F(x) \geq u) \]
In the Bernoulli example, what would the \(F^{-}(u)\) be for \(u = 0.2\) when \(\theta = 0.2\)?
Let’s tackle the discrete case, where \(X = F^{-}(U)\). First let’s work through the following identities: \[ F(F^{-}(u)) \geq u \quad \text{and} \quad F^{-}(F(x)) \leq x \] Which implies that the sets \[ \{(u,x): F^{-}(u) \leq x \} = \{(u,x): F(x) \geq u \} \] To see why, start with the inequality on the left: \[ \begin{aligned} F^{-}(u) & \leq x \\ F(F^{-}(u)) & \leq F(x) \\ u \leq F(F^{-}(u)) & \leq F(x) \end{aligned} \] where the last line follows from the first inequality and we are ultimately left with \(u \leq F(x)\) Now for the inequality on the right: \[ \begin{aligned} F(x) & \geq u \\ F^{-}(F(x)) & \geq F^{-}(u) \\ x \geq F^{-}(F(x)) & \geq F^{-}(u) \end{aligned} \] ultimately leading to \(x \geq F^{-}(u)\)
Thus, when we write our string of equalities, we can use the set equivalence: \[ \begin{aligned} P(X \leq x) & = P(F^{-}(U) \leq x) \\ & = P(F^{-}(U) \leq x) \\ & = P(U \leq F(x)) \\ & = F(x) \end{aligned} \]
This method is only useful if \(F^{-}(x)\) is easily computable.
Transformations
Another way to generate draws from a target distribution is to start with draws from one distribution and transform them somehow into the target distribution.
Example 2 (Gamma from exponentials) We saw how to get an exponential from a uniform. Suppose we now have a set of independent and identically distributed \(X_i \sim \text{Exp}(1), i = 1, \dots, d\).
Then we can generate a gamma random variable \(W\) with shape parameter \(d\) and scale \(1\) by setting \[ W = \textstyle\sum_{i=1}^d X_i \] Furthermore, if we want \(W \sim \text{Gamma}(d, \beta)\), we set \[ W = \beta\textstyle\sum_{i=1}^d X_i \]
Other transformations are useful when dealing with univariate normal random variables.
Multivariate normal random vectors
Let’s say we have \[X_i \sim \text{Normal}(0, 1), i = 1, \dots, d\] and we have a matrix \(A \in \R^{k \times d}\) with \(k \leq d\). If we arrange each \(X_i\) into a vector \(\mathbf{x}\) with the \(i^\mathrm{th}\) element \(X_i\) \[A X \sim \text{MultiNormal}(0, AA^T)\].
We can show this with the MGF. The MGF of a univariate standard normal is \(e^{t^2/2}\), and the MGF of a multivariate normal vector \(\mathbf{y}\) with variance covariance matrix \(\Sigma\) is \[ \Exp{e^{\mathbf{t}^T \mathbf{y}}} = e^{\frac{1}{2} \mathbf{t}^T\Sigma\mathbf{t} } \] Now we want to calculate the MGF of the vector \(\mathbf{y} = \mathbf{A} \mathbf{x}\)
\[ \begin{aligned} \Exp{e^{\mathbf{t}^T \mathbf{A} \mathbf{x}}} & = \Exp{e^{\sum_{i=1}^k t_i \sum_{j=1}^d A_{ij} x_j}} \\ & = \Exp{e^{\sum_{j=1}^d \sum_{i=1}^k t_i A_{ij} x_j}} \\ & = \prod_{j=1}^d \Exp{e^{\sum_{i=1}^k t_i A_{ij} x_j}} \\ & = \prod_{j=1}^d e^{(\sum_{i=1}^k t_i A_{ij})^2 / 2 } \\ & = e^{\sum_{j=1}^d (\sum_{i=1}^k t_i A_{ij})^2 / 2 } \\ & = e^{\sum_{j=1}^d \mathbf{t}^T A_{[:,j]}A^T_{[:,j]}\mathbf{t} / 2 } \\ & = e^{\mathbf{t}^T\lp \sum_{j=1}^d A_{[:,j]}A^T_{[:,j]}\rp\mathbf{t} / 2 } \\ & = e^{\mathbf{t}^T\mathbf{A}\mathbf{A}^T / 2 } \\ \end{aligned} \]
We can use this result to generate \(Y \sim \text{MultiNormal}(0, \Sigma)\) if we can decompose \(\Sigma\) into \(A A^T\).
One such decomposition is the matrix square root:
Definition 1 (Matrix square root) For an \(n\times n\) symmetric positive semidefinite matrix \(\Sigma\), we define the matrix square root as: \[ \Sigma = AA \]
In particular, we may set \(A = U \Lambda^{1/2} U^T\) where \(U\) are the eigenvectors of \(A\) with associated eigenvalues.
It is clear that \(\Sigma = A A\) from the fact that \(U\) has orthonormal columns (i.e. \(U^T U = I\)).
The procedure to generate multivariate normal random vectors from univariate normal draws is as follows:
- Store the eigendecomposition of \(\Sigma = U \Lambda U^T\)
- Draw \(X_i \sim \text{Normal}(0,1), i = 1,\dots, d\)
- Set \(Y= U \Lambda^{1/2} U^T X\)
The matrix square root is not the only decomposition that we can use. Any decomposition such that \(\Sigma = A A^T\) is acceptable. The most common and performant decomposition is called the Cholesky decomposition.
Definition 2 (Cholesky decomposition) For an \(n\times n\) symmetric positive definite matrix \(\Sigma\), we define the Cholesky decomposition as: \[ \Sigma = LL^T \] where \(L\) is the Cholesky factor and is a lower triangular matrix with a strictly positive diagonal.
Then, in order to generate \(Y \sim \text{MultiNormal}(0, \Sigma)\), we can replace the left multiplication by the matrix square root above by the Cholesky factor.
Any algorithm to decompose a matrix will be on the order of \(n^3\), but the Cholesky decomposition has a small constant multiplied by \(n^3\) while something like an \(LU\) decomposition, or an eigendecomposition have larger constants.
Furthermore, once you have the Cholesky factor computed, computing the inverse can be done in \(n^2\) operations. To see why consider the triangular system of equations \(L x = b\), with solution \(x = L^{-1}b\): \[ \begin{array}{rcrcrcrcrcl} L_{1,1}x_1 & & & & & & & = & b_1 \\ L_{2,1}x_1 & + & L_{2,2}x_2 & & & & & = & b_2 \\ \vdots & & \vdots & & \ddots & & & & \vdots \\ L_{m,1}x_1 & + & L_{m,2}x_2 & + & \cdots & + & L_{m,m}x_m & = & b_m \end{array} \] The first equation has solution \(x_1 = b_1 / L_{1,1}\), while the second equation can be solved for \(x_2\) using the prior solution for \(x_1\), and so on. This code was produced by Claude after I fed in an image from Wikipedia
This is called forward substitution, or a forward solve. The R function for this is forwardsolve.
If we instead had an upper triangular system of equations, we could use the same idea, but starting from the \(n^\mathrm{th}\) row instead of the first row. This is called backward substitution or, in R-code, backsolve.
Consider evaluating quadratic piece of the multivariate normal density: \[ (y - \mu)^T \Sigma^{-1}(y - \mu) \] Using the Cholesky decomposition \(\Sigma = L L^T\), we can evaluate the quadratic form quickly:
\[ \begin{aligned} (y - \mu)^T (L L^T)^{-1}(y - \mu) & = (y - \mu)^T L^{-T} L^{-1}(y - \mu) \\ & = (L^{-1} (y - \mu))^T L^{-1} (y - \mu)) \end{aligned} \] Thus we need to calculate \(x = L^{-1}(y - \mu)\), which we can do efficiently by solving the triangular system of equations.
Generating draws from mixtures
Often we can write a random variable as the marginal distribution of a mixture distribution. For example, we can write negative binomial random variable as an appropriate mixture of Poisson RVs with Gamma distributed rates.
More generally, we can use a somewhat trivial fact about densities to motivate our procedure. Suppose we want \(X \sim F\) where \(f(x)\) is the density corresponding to \(F\). We must find a bivariate density \(g(x,y)\) such that: \[ f(x) = \int_{\mathcal{Y}} g(x,y) dy \] and, crucially, \(g\) is easier to simulate from than directly simulating from \(f\). Often, our mixture distributions will take a more explicit form:
\[ f(x) = \int_{\mathcal{Y}} h_1(x \mid y) h_2(y) dy \] where \(h_1(x \mid y)\) is a conditional density that is easy to sample from.
Example 3 (Student’s t distribution) We can show that \[ \begin{aligned} Y \sim \chi^2_{\nu} \quad X \mid Y \sim \text{Normal}(0, \frac{\nu}{Y}) \\ \end{aligned} \] Then \[ X \sim t_\nu \]
Example 4 (Negative binomial) We can show that \[ \begin{aligned} Y \sim \text{Gamma}(n, (1-p)/p) \quad X \mid Y \sim \text{Pois}(Y) \\ \end{aligned} \] Then \[ X \sim \text{NegBinom}(n, p) \]