Chapter 6 Review (Chapter 1-5)

6.1 Simulation

Basic vectorized comparison

X <- c(1, 1, 5)
Y <- c(5, 5, 1)
X > Y # FALSE FALSE TRUE
## [1] FALSE FALSE  TRUE
sum(X > Y) # TRUE = 1, FALSE = 0
## [1] 1
mean(X > Y) # 1/3
## [1] 0.3333333

The code mean(X > Y) is computing \[\begin{equation*} \frac{1}{3}( I(X_1 > Y_1) + I(X_2 > Y_2) + I(X_3 > Y_3)), \end{equation*}\] where \(I(\cdot)\) is the indicator function, that is, \(I(X_1 > Y_1) = 1\) if \(X_1 > Y_1\) and \(I(X_1 > Y_1) = 0\) if \(X_1 \leq Y_1\).

Example: Let \(X \sim N(mean = 0, sd = 2)\) and \(Y \sim Exp(rate = 3)\). Estimate \(P(X > Y)\) using simulation.

n <- 10000 # number of simulations
X <- rnorm(n, 0, 2) 
Y <- rexp(n, 3)
mean(X > Y) 
## [1] 0.4309

The code mean(X > Y) is computing \[\begin{equation*} \frac{1}{n} \sum^n_{i=1} I(X_i > Y_i), \end{equation*}\] which is an approximation of \(P(X>Y)\).

Theory

Why the sample mean \(\frac{1}{n} \sum^n_{i=1} I(X_i > Y_i)\) approximates the required probability \(P(X>Y)\)?

Recall the strong law of large numbers (SLLN): Let \(X_1,\ldots,X_n\) be independent and identically distributed random variables with mean \(\mu:=E(X)\). Let \(\overline{X}_n:= \frac{1}{n} \sum^n_{i=1}X_i\). Then \[\overline{X}_n \stackrel{a.s.}{\rightarrow} \mu.\] Note: a.s. means almost surely (probability = 1). The above convergence means \(P(\lim_{n\rightarrow \infty} \overline{X}_n = \mu) = 1\). Note that (the expectation of an indicator random variable is the probability that the corresponding event will happen) \[ P(X>Y) = E(I(X>Y)).\] To apply SLLN, we just need to recognize the underlying random variable is \(I(X>Y)\) (which follows a Bernoulli distribution with parameter \(P(X>Y)\)). Then, with probability \(1\), the sample mean \[ \frac{1}{n} \sum^n_{i=1} I(X_i > Y_i) \rightarrow P(X>Y).\] The quantity on the LHS is what we compute in mean(X > Y) (note that we are using vectorized comparison). This is the reason why we can use mean(X > Y) to estimate \(P(X>Y)\).

Ex 1: Let X ~ N(mean = 2, sd =1), Y ~ Exp(rate = 2), Z ~ Unif(0, 4) (continuous uniform distribution on [0,4]). Estimate \(P(\max(X,Y)>Z)\).

Recall the difference between pmax and max:

# always try with simple examples when you test the usage of functions
x <- c(1, 2, 3)
y <- c(0, 5, 10)
max(x, y)
## [1] 10
pmax(x, y)
## [1]  1  5 10
n <- 100000
x <- rnorm(n, 2, 1)
y <- rexp(n, 2)
z <- runif(n, 0, 4)
mean(pmax(x, y) > z) 
## [1] 0.51502

Ex 2 Let X ~ N(mean = 2, sd = 1), Y ~ Exp(rate = 2), Z ~ Unif(0, 4) (continuous uniform distribution on [0,4]). Estimate \(P(\min(X,Y)>Z)\).

n <- 100000
x <- rnorm(n, 2, 1)
y <- rexp(n, 2)
z <- runif(n, 0, 4)
mean(pmin(x, y) > z) # what is the difference between pmin and min?
## [1] 0.11361

Ex 3 Let X ~ N(mean = 2, sd = 1), Y ~ Exp(rate = 2), Z ~ Unif(0, 4) (continuous uniform distribution on [0,4]). Estimate \(P(X^2 Y>Z)\).

n <- 100000
x <- rnorm(n, 2, 1)
y <- rexp(n, 2)
z <- runif(n, 0, 4)
mean(x ^ 2 * y > z)
## [1] 0.41158

Ex 4 Person \(A\) generates a random variable \(X \sim N(2, 1)\) and Person \(B\) generates a random variable \(Z \sim Unif(0, 4)\). If \(X < Z\), Person \(A\) will discard \(X\) and generate another random variable \(Y \sim Exp(0.5)\). Find the probability that the number generated by \(A\) is greater than that by \(B\).

n <- 100000
greater <- rep(0, n)
for (i in 1:n) {
  X <- rnorm(1, 2, 1)
  Z <- runif(1, 0, 4)
  if (X< Z) {
    Y <- rexp(1, rate = 0.5)
    greater[i] <- Y > Z
  } else {
    greater[i] <- 1 # 1 means A's no > B's no
  }
}

mean(greater)
## [1] 0.63698

Remark: you may find that the following code gives you almost the same answer. Why?

n <- 100000
x <- rnorm(n, 2, 1)
y <- rexp(n, 0.5)
z <- runif(n, 0, 4)
mean(pmax(x, y) > z)
## [1] 0.63904

Ex 5 Let X ~ N(mean = 2, sd = 1) and Y ~ Exp(rate = 2). Estimate \(E(min(X,Y))\).

n <- 100000
x <- rnorm(n, 2, 1)
y <- rexp(n, 2)
mean(pmin(x, y))
## [1] 0.445328

The code mean(pmin(x, y)) computes \[\begin{equation*} \frac{1}{n} \sum^n_{i=1} \min(X_i, Y_i), \end{equation*}\] which approximates \(E(\min(X,Y))\) by the SLLN.

6.2 Matrix

Ex6: Write a function called matrix_times_vector to compute \(X Y\), where \(X\) is a matrix and \(Y\) is a vector. The output should be a vector.

matrix_times_vector = function(X, Y) {
  as.vector(X %*% Y)
}

# e.g.
X <- matrix(1:12, 3, 4)
Y <- 1:4
X
##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12
Y
## [1] 1 2 3 4
matrix_times_vector(X, Y)
## [1] 70 80 90

Note:

  1. X %*% Y will return a matrix. We can use as.vector to change it into a vector.
  2. It is common to see the error non-conformable arguments. This is because the dimensions of your matrices/vectors do not match.
  • If you have a \(n\times p\) matrix \(A\) and \(m \times q\) matrix \(B\), you can do the matrix multiplication \(AB\) only if \(p = m\). In R, if this is not the case, there will be an error.
  • Similarly, if you have a vector \(d\) of length \(m\). You can do the matrix multiplication \(A d\) only if \(p = m\).

Ex7: Write a function called matrix_times_vector2 to compute \(X Y\), where \(X\) is a matrix and \(Y\) is a vector. The output should be a vector. However, you should check if the dimensions of the inputs are appropriate before you perform the calculation. Display an error message The dimensions do not match if this is not the case.

matrix_times_vector2 = function(X, Y) {
  p <- ncol(X)
  m <- length(Y)
  if (p == m) {
    return(as.vector(X %*% Y))   
  } else {
    cat("The dimensions do not match")
  }
}

# e.g.
X <- matrix(1:12, 4, 3)
Y <- 1:4
matrix_times_vector2(X, Y)
## The dimensions do not match

X <- matrix(1:12, 3, 4)
Y <- 1:4
matrix_times_vector2(X, Y)
## [1] 70 80 90

Example: if A is matrix and x is a vector, you can find the sums of the elements in A and x by sum(A) and sum(x) respectively.

x <- 1:10
sum(x)
## [1] 55

A <- matrix(1:10, 5, 2)
sum(A)
## [1] 55

Ex 8:

Find \(\sum^{5}_{x=1}\sum^{4}_{y=1} \frac{x}{x+y}\) without any loops.

(x <- matrix(1:5, nrow = 4, ncol = 5, byrow = TRUE)) # define and display at the same time using (...)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    2    3    4    5
## [2,]    1    2    3    4    5
## [3,]    1    2    3    4    5
## [4,]    1    2    3    4    5
(y <- matrix(1:4, nrow = 4, ncol = 5, byrow = FALSE))
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    1    1    1    1
## [2,]    2    2    2    2    2
## [3,]    3    3    3    3    3
## [4,]    4    4    4    4    4
x / (x + y) # vectorized operation
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.5000000 0.6666667 0.7500000 0.8000000 0.8333333
## [2,] 0.3333333 0.5000000 0.6000000 0.6666667 0.7142857
## [3,] 0.2500000 0.4000000 0.5000000 0.5714286 0.6250000
## [4,] 0.2000000 0.3333333 0.4285714 0.5000000 0.5555556
sum(x / (x + y))
## [1] 10.72817

6.3 Basic Operation

Example Let v be a vector of integers. Write a one-line R code to compute the product of all the even integers in v.

To illustrate how to solve the question step by step:

v <- -10:10 # begin writing your code by setting some integers
v %% 2 # find the remainder, if the remainder is 0, it is an even number
##  [1] 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0
v %% 2 == 0  # check which elements is 0
##  [1]  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE
## [16] FALSE  TRUE FALSE  TRUE FALSE  TRUE
v[v %% 2 == 0] # select the elements which are "TRUE"
##  [1] -10  -8  -6  -4  -2   0   2   4   6   8  10
prod(v[v %% 2 == 0]) # find the product
## [1] 0

# the result is 0. v may not be a good example to check if the code is correct
# let's change to some other vector
v <- 2:8
prod(v[v %% 2 == 0])
## [1] 384

2 * 4 * 6 * 8 #check 
## [1] 384

# now, the final answer is prod(v[v %% 2 == 0])

Ex 9a: Given that x = 1:100. Write a one-line R code to copmute \[\begin{equation*} S := 1^2 - 2^2 + 3^2 - \ldots + 99^2 - 100^2. \end{equation*}\]

x <- 1:100
sum((x[x %% 2 == 1]) ^ 2) - sum((x[x %% 2 == 0]) ^ 2)
## [1] -5050

Ex 9b (optional) Can you do Ex 9a without any program or calculator?

Ex 10: Write a function with inputs n and p to compute \[\begin{equation*} S(n, p) := 1^p - 2^p + 3^p -\ldots + (-1)^{n} (n-1)^p + (-1)^{n+1} n^p. \end{equation*}\]

snp <- function(n, p) {
  x <- 1:n
  return(sum((x[x %% 2 == 1]) ^ p) - sum((x[x %% 2 == 0]) ^ p))
}

6.4 Some basic plots in R

Example Write a function called my_summary_plot with input being two numeric vectors x,y and outputs a \(2 \times 2\) multi-frame plot with (i) histogram of x, (ii) histogram of y, (iii) scatter plot of y versu x, and (iv) boxplots of x and y.

# when you write a function, you can begin with some sample x and y
x <- rnorm(100, 0, 1)
y <- rnorm(100, 0, 1)
par(mfrow = c(2, 2))
hist(x)
hist(y)
plot(x, y)
boxplot(x, y)


# after that, you wrap the code into a function without defining x, y
my_summary_plot <- function(x, y){
  par(mfrow = c(2, 2))
  hist(x)
  hist(y)
  plot(x, y)
  boxplot(x, y)
}

# after you write the function, you should try if it will work or not
rm(list = ls()) # let's remove everything in the memory
my_summary_plot <- function(x, y){
  par(mfrow = c(2, 2))
  hist(x)
  hist(y)
  plot(x, y)
  boxplot(x, y)
}

x <- rnorm(100, 0, 1)
y <- rnorm(100, 0, 1)
my_summary_plot(x, y)


# try another example
x <- rnorm(100, 0, 1)
y <- 5 * x + rnorm(100, 0, 1)
my_summary_plot(x, y)