4.2.1 Inversion Method

Suppose \(X\) is a random variable with distribution function \(F_X\). Then, the distribution’s quantile function is defined as the generalized inverse \(F_X^{\leftarrow}\) of \(F_X\) where \(F_X^{\leftarrow}\) is defined as \[\begin{align*} F_X^{\leftarrow}(u) = \inf \{ x \in \mathbb{R}\, \vert\, F_X(x) \geq u \} \end{align*}\] for all \(u \in [0, 1)\).

The inversion method relies on the fact that the random variable \(F_X^{\leftarrow}(U)\), where \(U \sim\) \(U(0, 1)\), has distribution function \(F_X\), i.e. follows the same distribution as \(X\). Consequently, we have found a method that justifies the following algorithm:

1. Compute quantile function \(F_X^{\leftarrow}\) of desired distribution
2. Draw sample \(u\) from standard-uniform distribution
3. \(x = F_X^{\leftarrow}(u)\) can now be thought of as a sample from distribution of \(X\)

This can be easily simulated. Let us do that and compare the corresponding histogram to the theoretical density. For the histogram we will have to adjust the default’s binwidth in order to catch the behavior close to zero better.

set.seed(123)
lambda <- 2

simu_data <- tibble(
  u = runif(100000),
  x = -log(1 - u) / lambda
)

density_data <- tibble(
  t = seq(0.01, max(simu_data$x), 0.01),
  density = dexp(t, rate = lambda) # density of exponential
)

simu_data %>% 
  ggplot() +
  geom_histogram(aes(x, y = ..density..), binwidth = 0.01) +
  geom_line(data = density_data, aes(t, density), col = "red")
Copy

Similarly, we can also simulate realizations from discrete random variables.

Again let us use this and compare the histogram to the real distribution.

set.seed(123)
p <- 0.3

density_data <- tibble(
    x = c(0, 1), 
    y = c(1- p, p)
)

tibble(
  u = runif(100000),
  x = between(u, 1-p, 1) * 1 # * 1 converts logical to double
) %>% 
  ggplot() +
  geom_bar(aes(x, y = ..prop..)) +
  geom_point(
    data = density_data, 
    aes(x, y),
    col = "red",
    size = 2
  )
Copy

4.2.2 Acceptance-Rejection Method

Imagine that you use the previous methods to simulate independent realizations \(u_i, x_i\), \(i = 1, \dots, N\), \(N \in \mathbb{N}\), of the random variables \(U \sim U(0,1)\) and \(X \sim \text{Exp}(\lambda)\), \(\lambda > 0\), respectively. Imagine further that you now use these realisations to plot the points \((x_i, u_i f_X(x_i))\) and \(f_X\), where \(f_X\) is the density of \(X\). This may look like this.

As you can see, these these points fill up the area between the x-axis and the density function of \(X\). What you may not see immediately but what is true if you think about it, is the fact that the amount of points for a given value of \(x\) is related to \(f(x)\). Obviously, this comes from the fact that we simulated realisations \(x_i\) of \(X\).

Similarly, if we scale the y-coordinates of the points by a factor \(k > 1\), i.e. plot \((x_i, k u_i f_X(x_i))\), then the points will fill the area under \(kf_X(x)\). Thus, we could use this fact to generate realizations of a different random variable \(Y\) whose density \(f_Y\) we can bound from above by \(kf_X\). We do this by filling the area under \(kf_X\) with the help of realizations of \(X\) as above but rejecting points that do not also fall within the area under \(f_Y\). We check the latter by checking whether \(\frac{k u_i f_X(x_i)}{f_Y(x_i)} < 1\). This procedure is known as acceptance-rejection method.

Imagine that we want to simulate a random variable \(Y\) that follows a Gamma distribution with parameters 3/2 and 1, i.e. \[\begin{align*} f_Y(y) = \frac{2}{\sqrt{\pi}} \sqrt{y} e^{-y} \end{align*}\] for \(y > 0\).

Figure 4.2: Density of Gamma distribution (in green) and rejection area (in red)

Now, to implement the acceptance-rejection method, we will have to use if-conditions to check whether a simulated point is in the rejection area. In R, the general structure of if-conditions is given as follows

if (<condition>) {
  <statement if condition is true>
} else {
  <statement if condition is false>
}
Copy

For short if-conditions and corresponding statements R allows it to write everything in one line. Personally, I dislike that style but it may be worth noting this style of writing is possible if you so desire. Further, we will need while-loops to keep generating random variables until we can accept one. In R, the general structure of while-loops is given as follows

while (<condition>) {
  <statement while condition is true>
}
Copy

Here, we will use while and if in conjunction to loop until we can accept a realization and then set the looping condition to FALSE.

set.seed(123)
lambda <- 2 / 3
k <- 3 * sqrt(3 / (2 * pi)) * exp(-1/ 2)
# Boolean for rejection, initial state needs to be TRUE
rejected <- TRUE
while (rejected) {
  # Inversion method for Exp(2 / 3) realization
  u <- runif(1)
  x <- -log(1 - u) / (lambda)
  
  # independent realization of x
  u <- runif(1)
  # Check whether accepted or rejected
  f_X <- lambda * exp(- lambda * x)
  f_Y <- 2 / sqrt(pi) * sqrt(x) *  exp(-x)
  if (k * u * f_X / f_Y < 1) {
    # Accepted so set rejected to FALSE
    rejected <- FALSE
  }
}
x
#> [1] 0.5086263
Copy

Similarly, we can simulate a lot of realizations:

set.seed(123)
N <- 10000
xs <- numeric(N)
for (i in seq_along(xs)) {
  # Boolean for rejection, initial state needs to be TRUE
  rejected <- TRUE
  while (rejected) {
    # Inversion method for Exp(2 / 3) realization
    u <- runif(1)
    x <- -log(1 - u) / (lambda)
    
    # independent realization of x
    u <- runif(1)
    # Check whether accepted or rejected
    f_X <- lambda * exp(- lambda * x)
    f_Y <- 2 / sqrt(pi) * sqrt(x) *  exp(-x)
    if (k * u * f_X / f_Y < 1) {
      # Accepted so set rejected to FALSE
      rejected <- FALSE
    }
  }
  xs[i] <- x
}

# See if realizations look like gamma density
x <- seq(0, max(xs), 0.05)
f_Y <- 2 / sqrt(pi) * sqrt(x) *  exp(-x)
ggplot() +
  geom_histogram(aes(x = xs, y = ..density..), binwidth = 0.1) +
  geom_line(aes(x, f_Y), col = "red", size = 1) +
  theme_bw()
Copy

This looks pretty decent. So, we can assume that the acceptance-rejection method delivered the results we were wishing for.

However, notice that depending on what density \(f_X\) and constant \(k\) we use to bound the density \(f_Y\),it can take a lot of trials before we can accept one realization. For instance, imagine that in our previous example we would have used a larger value \(k\). Immediately, this would lead to a larger rejection area in Figure 4.2.

Therefore, a good choice of \(k\) and \(f_X\) is crucial for the efficiency of the acceptance-rejection method. Here, however, we did not have any other means to simulate the realizations from the gamma distribution as we do not have a closed form of the quantile function. Another popular distribution for which this is the case is the normal distribution which we want to talk about next.

4.2.3 (Correlated) Normal Random Variables

Let us start by mentioning that it is numerically possible to use the inversion method to simulate normal random variables even though we do not have a closed formula for the quantile function. Similarly, one may construct an algorithm using the acceptance-rejection method. Though, neither is what we want to do here.

Here, we will simulate independent normal random variables through the so-called Box-Muller transformations. For independent realizations \(u\) and \(v\) from the standard uniform distribution these transformations are defined as \[\begin{align*} x &= \sqrt{-2 \log (u)} \cos(2 \pi v), \\ y &= \sqrt{-2 \log (u)} \sin(2 \pi v) \end{align*}\]

Using arguments based on the polar decomposition of a vector \((x, y)\) it may be shown that here, \(x\) and \(y\) may be regarded as two independent realizations from the standard normal distribution \(\mathcal{N}(0, 1)\)⁴⁷. Let us check this empirically.

set.seed(123)
N <- 10000
# Simulate Realizations
simu_dat <- tibble(
  u = runif(N),
  v = runif(N),
  x = sqrt(-2 * log(u)) * cos(2 * pi * v),
  y = sqrt(-2 * log(u)) * sin(2 * pi * v)
) %>% 
  pivot_longer(
    cols = -c(u, v),
    names_to = "coordinate",
    values_to = "value"
  ) 

# Compute theoretical density
theo_dat <- tibble(
  x = seq(-3, 3, 0.1),
  f = dnorm(x) # standard normal density
)
  
# Plot histograms
ggplot() +
  geom_histogram(
    data = simu_dat,
    aes(
      x = value, 
      y = ..density..,
      fill = coordinate
    ),
    show.legend = F
  ) +
  geom_line(
    data = theo_dat,
    aes(x, f),
    col = "black", 
    size = 1
  ) +
  facet_wrap(vars(coordinate)) +
  theme_bw()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Copy

This is how standard normal distributions are supposed to look. But are they also independent? To find out, we will need to look at both coordinates simultaneously. Theoretically, if the realizations can be thought of as stemming from two independent normal random variables, then the vector of the vectors \((x_i, y_i)\) should correspond to realizations of a multivariate normal distribution.

Figure 4.3: Bivariate standard normal density for \(\rho \in \{-0.8, 0, 0.8\}\)

As Figure 4.3 suggests, our simulated vectors \((x_i, y_i)\) should be concentrated circularly around the origin. Let’s check that.

simu_dat <- simu_dat %>% 
  select(coordinate, value) %>% 
  # Create unique identifier column for pivot_wider
  mutate(i = rep(1:N, each = 2)) %>% 
  pivot_wider(
    id_cols = 3,
    names_from = coordinate,
    values_from = value
  )

simu_dat %>% 
  ggplot() +
  geom_point(aes(x, y))  +
  coord_equal() # to see circle better
Copy

This looks pretty circular but there is also a lot of overplotting in this picture. Let’s use geom_bin2d() as a two-dimensional version of a histogram to see what is going on close to the origin.

simu_dat %>% 
  ggplot() +
  geom_bin2d(aes(x, y))  +
  coord_equal() # to see circle better
Copy

This looks like it is supposed to if the vectors \((x_i, y_i)\) are truly realizations of a two-dimensional standard normal vector with correlation \(\rho = 0\). So, for now let’s accept this as proof that our simulation procedure works.

Now that we can simulate independent (multivariate) standard normal random variables we can easily transform the realizations such that they can be considered as coming from an arbitrary normal distribution. In the univariate case this is an easy linear transformation of the form \(y = \sigma x + \mu\), where \(\mu \in \mathbb{R}\) and \(\sigma > 0\) are the desired expectation and standard deviation respectively.

mu <- 5
sigma <- 3

theo_dat <- tibble(
  x = seq(mu - 3 * sigma, mu + 3 * sigma, 0.1),
  f = dnorm(x, mean = mu, sd = sigma)
)

simu_dat %>% 
  mutate(x = sigma * x + mu) %>% 
  ggplot() +
  geom_histogram(aes(x = x, y = ..density..)) +
  geom_line(
    data = theo_dat,
    aes(
      x = x, 
      y = f
    ),
    col = "red",
    size = 1
  ) +
  theme_bw()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Copy

In the multivariate case, this is a bit more complicated. As we have learned in the quick refresher on multivariate normal vectors \(X\), we can rewrite \(X\) as \(X = \mu + AZ\). Here, \(Z\) is a standard normal random vector (which we can simulate now), \(\mu\) is the vector of desired means and \(A\) is a matrix and related to the desired covariance matrix \(\Sigma\) via \(AA^T = \Sigma\).

Thus, we can transform the mean of a random standard normal vector by simply adding the vector \(\mu\). As for the desired covariance structure, we need to decompose \(\Sigma\) into \(AA^T\) and multiply the resulting matrix \(A\) with \(Z\). This may seem complicated but luckily, this decomposition is famously known as cholesky decomposition and R has the built-in chol() command to compute it.

Consequently, we can easily implement the transformation using matrices via matrix(), applying chol() and performing the corresponding matrix multiplication via %*%. Let’s try this to generate realizations from a two-dimensional random normal vector \((X_1, X_2)\) with \(\mathbb{E}X_i = 0\), \(\text{Var}(X_i) = 1\), \(i = 1, 2\), and correlation \(\rho = 0.8\).

# Mean of components
mu <- c(0, 0)
# Standard deviations of components
sds <- c(1, 1)
# Correlation 
rho <- 0.8
# Desired covariance matrix
Sigma <- matrix(
  c(sds[1]^2, 
    rho * prod(sds), 
    rho * prod(sds), 
    sds[2]^2),
  nrow = 2,
  ncol = 2
)
# Cholesky decomposition
A <- chol(Sigma)

N <- 100000
x <- numeric(N)
y <- numeric(N)
# Compute standard normal realisations
set.seed(123)
tib_Z <- tibble(
  u = runif(N),
  v = runif(N),
  x = sqrt(-2 * log(u)) * cos(2 * pi * v),
  y = sqrt(-2 * log(u)) * sin(2 * pi * v)
)

# Transform each realisation
# And save each coordinate into corresponding vector
for (i in seq_along(x)) {
  z <- c(tib_Z[[i, 'x']], tib_Z[[i, 'y']])
  x[i] <-  (mu + A %*% z)[1]
  y[i] <-  (mu + A %*% z)[2]
}
Copy

Now, we can check whether our realizations look similar to the right pane in Figure 4.3. If so, we can assume that our simulation worked.

tibble(
  x = x,
  y = y
) %>% 
  ggplot() +
  geom_bin2d(aes(x, y)) +
  scale_fill_distiller(palette = "Spectral") +
  theme_bw() 
Copy

4.2.4 R Commands

Naturally, R already comes with routines to simulate the most common distributions. So while it is interesting to learn how simulations can be implemented, often we can use R’s own procedures. Except, of course, in the exercises of this chapter. There, you will have to generate realizations from certain distributions yourself.

Nevertheless, let us briefly talk about R’s simulation procedures here. All of R’s simulation start with r (for realization) followed by the distribution’s abbreviation, e.g. runif(), rexp(), rbinom(), rnorm(), etc.

To find out how to adjust the parameters of a specific distribution, e.g. \(\lambda > 0\) of \(\text{Exp}(\lambda)\), you will have to look into the function’s documentation. What you will find out, if you do so, is that on the documentation page there are multiple functions described. In fact, for each distribution R has not only a function with r prefix but it also uses the prefixes q (for quantile), d (for density) and p (for probability meaning the distribution function).

Sometimes you will have to be careful with R’s own functions though. That is because some distributions such as e.g. the exponential, gamma or geometric distribution are defined with different parameterization in the literature.

For instance, in the case of a random variable \(X\) following a geometric distribution with parameter \(p \in (0, 1)\), there are two options. \[\begin{align*} \mathbb{P}(X = k) &= p(1-p)^k, \quad k = 0, 1, \dots\\ \mathbb{P}(X = k) &= p(1-p)^{k-1}, \quad k = 1, 2, \dots\\ \end{align*}\] Thus, you will have to check the documentation if the command you want to use really refers to the distribution you want to simulate.

4.3.1 For-loops and if-conditions

1. Use a for-loop to print all numbers between 1 and 100 that are not divisible by 3 and 4.
2. Do the same thing without a for-loop. Hint: A vector consisting of TRUEs and/or FALSEs can be used as indices to access a vector’s elements. The numbers also do not have to be printed on separate lines

4.3.2 Efficiency of Vector Creation

4.3.3 Binomial distribution

Write a function my_binomial(n, size, prob) that generates n realizations from a binomial distribution with size size and probability prob. Use it to simulate 100000 realizations of a binomial distribution with parameters \(10\) and \(0.4\) and plot a histogram based on these realizations together with the theoretical probabilities.

You may only use the runif() for simulating realizations of random variables but you may use dbinom() and pbinom() to compute probabilities.

Hint: You may have to tweak the histogram options to get the result you expect.

4.3.4 Geometric Distribution

Write a function my_geom(n, prob) that generates n realizations from a geometric distribution with probability prob. In this exercise, the density of the geometric distribution is defined as \[\begin{align*} \mathbb{P}(X = k) &= p(1-p)^{k-1}, \quad k = 1, 2, \dots, \end{align*}\] where \(p\) is the probability of success. Use my_geom() to simulate 1000 realizations of a geometric distribution with parameter \(0.4\) and plot the empirical distribution function of your sample against the theoretical distribution function.

Here, you may use rgeom() but you need to make sure that your realizations correspond to the above density. Also, the function ecdf() may be helpful.

4.3.5 Cauchy Distribution

Write a function my_Cauchy(n, location, scale) that generates n realizations from a Cauchy distribution with location parameter location and scale parameter scale. For information on the Cauchy distribution refer to Wikipedia.

Use your function to simulate 100000 realizations of a Cauchy distribution with location 2 and scale 5 and plot an estimate of the density based on your simulated sample against the theoretical density. You may use dcauchy() and runif() but not rcauchy() here.

Hint: It is possible to transform a standard uniform variable into a standard Cauchy variable (location 0 and scale 1). Also, a standard Cauchy variable can be transformed to an arbitrary Cauchy distribution.

Hint 2: You may have to filter out realizations that are extremely large in order to estimate the density nicely with geom_density(). Here it could be useful to stick to the window \([-50, 50]\).

4.3.6 Normal distribution

Write a function my_Gaussian(n, mu, sigma) that based on an acceptance-rejection method and the standard Cauchy density generates n realizations from a Gaussian distribution with location parameter mu and scale parameter sigma.

Note that it is enough to generate standard Gaussian random realizations via acceptance-rejection and then transform those realizations accordingly. Further, the parameter \(k\) of the acceptance-rejection method does not have to be computed explicitly. You can simply compare the densities visually and choose a \(k\) that fits.

Use your function to generate 100000 realizations of a Gaussian random variable with mean 5 and standard deviation 2 and recreate the following plot.

Hint: The package patchwork may help with arranging the plots.