Reading

Vectorization

Vectorization is a programming technique used to avoid explicit loops in order to improve the performance and, often, readability of code.

In its simplest form, vectorization can refer to using built in vectorized functions. That is, rather than looping over a vector to compute its sum, we instead use the vectorized function sum().

Demonstration 1

## [1] TRUE

The vectorized code is not only easier to type and read but is also faster.

Demonstratation 2

Demonstration 2.Comparing the run time of a loop to a vectorized sum.
expr min lq mean median uq max neval
f0(10000) 434794 438305 532478.1 466764.5 559634.5 5152385 1000
f1(10000) 522 755 8996.4 1073.5 4201.5 3536612 1000

Why vectorization is faster has to do with the fact that R is an interpreted language and not a compiled one. The difference in time you see here is primarily due to two factors:

  1. R needs to repeatedly interpret what your code means for each iteration of the loop;

  2. each iteration of the for loop in f0 requires subsetting x.

Monte Carlo Estimation

Vectorization can be particularly useful in Monte Carlo studies where we might otherwise be inclined to use explicit loops. We will look at some examples after an introduction to Monte Carlo estimates.

In statistics and data science we are often interested in computing expectations of random outcomes of various types. When analytical expectations are unavailable or cumbersome to compute, it can be useful to obtain Monte Carlo approximations by simulating a random process and then directly averaging the values of interest.

This works because the sample average is generally a good estimate of the corresponding expectation:

\[ {\bar{\theta}}_{n} := \sum_{i=1}^n X_i / n \to_p \theta := E[X]. \] \[ \bar{\theta}_n \to_p \theta \] In fact, assuming our data are independent and identically distributed (iid) from a distribution with finite variance, we can characterize the rate of convergence of a sample average to its population counterpart using the central limit theorem (CLT),

\[ \sqrt{n}(\bar{\theta}_n - \theta) \to_{d} N(0,\sigma^2) \] where \(\sigma^2 = E[X^2] - E[X]^2\) is the variance of the underlying distribution. This can be useful for constructing approximate confidence intervals for the Monte Carlo error.

Distribution functions

There are functions in R for simulating from many common distributions. Here are a few:

  • rnorm() - Normal
  • runif() - Uniform
  • rt() - the t-distribution
  • rexp() - Exponential
  • rpois() - Poisson

Another useful function in R is sample() for sampling from a finite set of values, i.e. the discrete uniform distribution or any finite probability mass function.

As an aside, you should be aware that each of the distribution families above have corresponding d*, p*, and q* functions for computing densities, percentiles (CDFs), or quantiles (inverse CDFs) for each distribution.

Random Seeds

When we call one of the r* functions to generate random draws from a distribution, R relies on a pseudo-random number generate to generate from \(U(0,1)\) and produce the results. Thus the outcome of these calls depends on the current state of the generator. It is sometimes desirable to reproduce exactly the same pseudo-random sequence. You can do this by fixing the random seed using set.seed() which takes an integer argument. The function RNGkind() can be used to display or set the random number generator.

## [1] "Mersenne-Twister" "Inversion"        "Rejection"
## [1] FALSE
## [1] TRUE

Basic Example 1

As a quick example, let’s use these functions to compute percentiles for t-distributions with various degrees of freedom. Let \(\theta_q\) be the parameter of interest,

\[ \theta_q := F(q) = \int_{-\infty}^{q} f(x) dx = \int_{-\infty}^{\infty} 1[x \le q]f(x) dx \] where \(F(\cdot)\) is the CDF and \(f(\cdot)\) the PDF of a given t-distribution.

In this case, our Monte Carlo estimate of \((\theta_{-1.96}, \theta_{1.96})\) is \(\bar{\theta}=\) (0.0704, 0.9269). The actual values are \((\theta_{-1.96}, \theta_{1.96})\) = (0.0724261, 0.9275739).

Basic Example 2

Suppose we are interested in computing the following integral where \(\phi\) is the standard normal density function:

\[ \int_{-\infty}^\infty [\sin(x)-\cos(x)]^2\phi(x) dx. \] We can recast this as the expectation below,

\[ E[h(X)], \quad h(x) = [\sin(x)-\cos(x)]^2, \quad X \sim N(0,1). \] The following R code provides a Monte Carlo estimate,

## [1] 1.00058

Compare this to an estimate using numerical integration,

## 1 with absolute error < 9.4e-05

These values are fairly close to the analytic solution based on the identity \([\sin(x)-\cos(x)]^2 = 1 - \sin(2x)\) and the symmetry about zero of both \(\sin(\cdot)\) and \(\phi(\cdot)\). Suppose \(X \sim N(0,1)\), then

\[ \begin{align} E[(\sin(X)-\cos(x))^2], &= E[1-\sin(2X)], \\ &=1 - E[\sin(2X)], \\ & = 1 - 0, \\ & = 1. \end{align} \]

The code below will produce a plot that illustrates why \(E[\sin(a*X)] = 0\) for \(X \sim N(0,1)\)

**Figure 1.** *Illustrating why E[sin(X)] = 0 for X ~ N(0,1).*

Figure 1. Illustrating why E[sin(X)] = 0 for X ~ N(0,1).

Examples

The examples in this section are taken directly from Professor Shedden’s 2016 course notes which you can find here.

Simulation Study for Nominal Confidence Intervals

In our first example, we will investigate the coverage probability of nominal 95% confidence intervals when the data does not come from a Normal (Gaussian) distribution.

In the first example, we will assume the data come from an exponential distribution with mean one. The strategy here is to generate many (mcrep) data sets of size n. For each data vector, we then calculate a nominal 95% confidence interval for the mean and check whether this interval contains the true value of one.

## [1] 0.9308

Since coverage is binary with a fixed probability \(p\), the number of intervals that cover one (“successes”) in our study is Binomial(mcrep, \(p\)). We can use this fact to estimate the Monte Carlo error which represents the uncertainty in our estimate from the chosen number of replications.

In this case the estimated coverage is 0.93 (0.926, 0.936).

Broadcasting

We have previously discussed the need to be careful about R recycling values when performing arithmetic on arrays with differing dimensions. The formal name for this is broadcasting and it can be quite useful for vectorization.

Generally, when you perform a simple arithmetic operation “*” on arrays with the same dimensions the operation is applied point-wise so that (A*B)[i,j] = A[i,j]*B[i,j]. However, if the objects A and B have different dimensions, but flattened lengths with one a multiple of the other the values can be broadcast like this:

##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6
##      [,1] [,2]
## [1,]    0    3
## [2,]    0    3
## [3,]    0    3

Compare this to:

##      [,1] [,2]
## [1,]    0    3
## [2,]    0    3
## [3,]    0    3

In this case, if X has dimensions m and n and v length k we see that (X-v)[i,j] = X[i,j] - v[{(j-1)*m+i mod k} + 1].

This is useful for centering each row of a matrix or computing row-wise variances.

We can compare this to direct computation using an implicit loop via the apply function:

## [1] TRUE

Let’s compare the timing of these two approaches:

## Apply: 1.465s, Vectorized: 0.034s, Ratio: 43.1.
## [1] TRUE

Correlation Coefficients

In the next example, we consider the problem of computing the correlation coefficient between many vectors \(x\) and a single outcome \(y\). This can be useful when screening a large set of potential predictors for a relationship with an outcome of interest \(y\).

First, we generate data that has n observations, m predictors, and expected correlation r between each predictor and \(y\).

Now, we can compute the correlations between each row of xmat and y and compare approaches.

## [1] TRUE
## [1] TRUE
## [1] TRUE
## 1: 0.276 s 
## 2: 0.260 s 
## 3: 0.226 s 
## 4: 0.024 s

The fourth approach using linear algebra and broadcasting is by far the most efficient here. All approaches are much more efficient than computing choose(10001, 2) correlations when we only need 10,000.

While we should keep in mind that this was a single trial and not a formal comparison with replicates, a difference of this size is still meaningful. We should also be aware that one of the reasons the other approaches are slower is the time needed to allocate memory for (larger) intermediate objects.

Garbage Collection

The time needed to allocate memory can be influenced by something called garbage collection 🗑 ♻️. In languages like R that bind names to values, garbage collection refers to freeing up memory that is no longer needed because there are no longer any names pointing to it.

Garbage collection is triggered automatically when R needs more memory. As this could happen in the middle of something you are timing, you may get more consistent results if you explicitly force garbage collection prior to starting the timer. However, you rarely need to consider this elsewhere and not everyone agrees it is a best practice for timing comparisons.

##           used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells  916691 49.0    1674622 89.5         NA  1674622 89.5
## Vcells 2892342 22.1   11381007 86.9      16384 11381007 86.9
##           used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells  916737 49.0    1674622 89.5         NA  1674622 89.5
## Vcells 2892396 22.1   11381007 86.9      16384 11381007 86.9
##           used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells  916799 49.0    1674622 89.5         NA  1674622 89.5
## Vcells 2892517 22.1   11381007 86.9      16384 11381007 86.9
##           used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells  916876 49.0    1674622 89.5         NA  1674622 89.5
## Vcells 2892595 22.1   11381007 86.9      16384 11381007 86.9
## 1: 0.219 s 
## 2: 0.236 s 
## 3: 0.146 s 
## 4: 0.016 s

Exercises

  1. Suggest a way to further speed up the broadcasting approach. Test it.
  2. Repeat the exercise above with observations in rows and predictors in columns, i.e. when x has dimension n by m rather than m by n.

Functional programming in Monte Carlo Studies

Recall our original example of computing coverage probabilities of nominal confidence intervals for non-Gaussian data.

If we wanted to carry out similar studies for many distributions, we may wish to write a function whose body carries out the simulation study using configurable parameters.

Here is the original example with a few changes.

Below we incorporate the simulation into a function with parameters for simulation settings. Here we use the special argument ... to pass additional parameters to rgen by name.

Now we can use estimate_nominal_coverage for multiple simulations.

## [1] 0.9309
## [1] 0.9462
## [1] 0.9584

This could be useful, say, for exploring how the mean or another parameter impacts the coverage probability for a particular distribution.

Exercises

  1. Use vectorization to improve the efficiency of estimate_nominal_coverage() and compare your improvemnts to the original.

  2. Use the function estimate_nominal_coverage() defined above to explore how the nominal coverage changes with the sample size n for various distributions.