1.

a.

(Not graded)

i.

### Initialization
n <- 100

### Initialization
out <- 0
### for loop
for (i in 1:n){
  out <- out + 1 / i
}
out

[1] 5.187378

### Vectorized operation
i <- 1:n
sum(1 / i)

[1] 5.187378

ii.

### Set a large n
n <- 1e7

### for loop time
system.time({
  out <- 0

  for (i in 1:n){
    out <- out + 1 / i
  }
  out
})

   user  system elapsed 
   0.17    0.00    0.17

### Vectorized operation time
system.time({
  i <- 1:n
  sum(1 / i)
})

   user  system elapsed 
   0.06    0.02    0.07

b.

i.

### Initialization
n <- 100
m <- 50

### Initialization
out <- 0
### Nested for loop
for (i in 1:n){
  for (j in 1:m){
    out <- out + 2^(- i * j)
  }
}
out

[1] 1.606695

### Vectorized operation 1
### Matrix multiplication
ind <- matrix(1:n) %*% matrix(1:m, 1)
sum(2^(- ind))

[1] 1.606695

### Vectorized operation 2
### Scaler multiplication
ind <- matrix(1:n, nrow = n, ncol = m) * matrix(1:m, nrow = n, ncol = m, byrow = TRUE)
sum(2^(- ind))

[1] 1.606695

### Vectorized operation 3
ind <- outer(1:n, 1:m, "*")
sum(2^(- ind))

[1] 1.606695

ii.

### Initialization
n <- 1e4
m <- 1e3

### Nested for loop time
system.time({
  out <- 0
  
  for (i in 1:n){
    for (j in 1:m){
      out <- out + 2^(- i * j)
    }
  }
  out
})

   user  system elapsed 
   1.16    0.00    1.15

### Vectorized operation 1 time
system.time({
  ind <- matrix(1:n) %*% matrix(1:m, 1)
  sum(2^(- ind))
})

   user  system elapsed 
   0.93    0.00    0.93

### Vectorized operation 2 time
system.time({
  ind <- matrix(1:n, nrow = n, ncol = m) * matrix(1:m, nrow = n, ncol = m, byrow = TRUE)
  sum(2^(- ind))
})

   user  system elapsed 
   1.00    0.02    1.01

### Vectorized operation 3 time
system.time({
  ind <- outer(1:n, 1:m, "*")
  sum(2^(- ind))
})

   user  system elapsed 
   0.92    0.00    0.92

c.

### while loop 1
n <- 1
out <- 0

while (out <= 11.5){
  out <- out + 1 / n
  n <- n + 1
}
### Need to subtract 1 because counts 1 more when stop
(n <- n - 1)

[1] 55425

### Quick check
### Not exceed
sum(1 / (1:(n - 1)))

[1] 11.49999

### Exceed
sum(1 / (1:n))

[1] 11.50001

### while loop 2
n <- 0
out <- 0

while (out <= 11.5){
  n <- n + 1
  out <- out + 1 / n
}
n

[1] 55425

### Quick check
### Not exceed
sum(1 / (1:(n - 1)))

[1] 11.49999

### Exceed
sum(1 / (1:n))

[1] 11.50001

### while loop 3
n <- 1
out <- 0

while (TRUE){
  out <- out + 1 / n
  if (out <= 11.5){
    n <- n + 1
  } else {
    break
  }
}
n

[1] 55425

### Quick check
### Not exceed
sum(1 / (1:(n - 1)))

[1] 11.49999

### Exceed
sum(1 / (1:n))

[1] 11.50001

(Not graded)

### without a loop
n <- 1e5
out <- cumsum(1 / (1:n))
which(out > 11.5)[1]

[1] 55425

d.

### Need double quotes
Q1d <- '"Hello, Iowa. Congratulations to the Iowa hawkers. That was a big win today. I’m thrilled to be back. That was a big win. But I am thrilled to be back especially on such great news as that, that was a big one. You’ve been a great school, a great team, a great tradition, really an amazing job and it all started right here and we’re going to keep it here. Number one, right? We’re going to keep it here. The fairgrounds so they broke the record tonight in the history of the fairgrounds. I don’t know how old it is, but in the history of the fairgrounds, this is the most people they’ve had. So thank you very much."'
writeLines(Q1d)

"Hello, Iowa. Congratulations to the Iowa hawkers. That was a big win today. I’m thrilled to be back. That was a big win. But I am thrilled to be back especially on such great news as that, that was a big one. You’ve been a great school, a great team, a great tradition, really an amazing job and it all started right here and we’re going to keep it here. Number one, right? We’re going to keep it here. The fairgrounds so they broke the record tonight in the history of the fairgrounds. I don’t know how old it is, but in the history of the fairgrounds, this is the most people they’ve had. So thank you very much."

### while loop
i <- 1
while (stringr::str_sub(Q1d, i, i+4) != "great"){
  i <- i + 1
}
i

[1] 171

### Quick check
stringr::str_sub(Q1d, i, i+4)

[1] "great"

### repeat loop
i <- 1
repeat {
  if (stringr::str_sub(Q1d, i, i+4) == "great") {
    break
  } else {
    i <- i + 1
  }
}
i

[1] 171

### Quick check
stringr::str_sub(Q1d, i, i+4)

[1] "great"

2. (Optional)

(Not graded)

### Number of non-zero a
### Set na <- 5 to get the formula shown in the question
### For demonstration, na is set large here to approximate the limit
na <- 100
a <- (seq(1, by=2, length.out=na))^2
b <- c(3, rep(6, na - 1))
denom <- 6 + 0

for (ia in 1:na){
  denom <- rev(b)[ia] + rev(a)[ia] / denom
}
denom

[1] 3.141592

3.

(Not graded)

mymedian <- function(x){
  ### Remove NA values
  x <- na.omit(x)
  
  x.sort <- sort(x)
  ### Check the number of entries in x to decide where the median is
  ifelse(length(x) %% 2 == 0, mean(x.sort[c(length(x)/2, length(x)/2 + 1)]), x.sort[ceiling(length(x)/2)])
}

### Unit test for NA and odd number
x1 <- c(1:5, NA)
testthat::expect_equal(mymedian(x1), median(x1, na.rm=TRUE))

### Unit test for unsorted order and even number
x2 <- c(2, 4, 1, 3, 5, 1)
testthat::expect_equal(mymedian(x2), median(x2, na.rm=TRUE))

4.

(Not graded)

a.

fp <- function(x, p = 2, deriv = FALSE){
  if (deriv) {
    if (p == 1){
      out <- sign(x)
    } else if (p > 1) {
      out <- p * (abs (x))^(p - 1) * sign(x)
    }
  } else {
    out <- (abs(x))^p
  }
  
  out
}

b.

i.

x <- -2: 2
testthat::expect_equal(length(fp(x, p = 1)), length(x))

ii.

x <- -2: 2
p <- 2
testthat::expect_equal(fp(x, p), abs(x)^p)

iii.

x <- -2: 2
p <- 1
testthat::expect_equal(fp(x, p), abs(x)^p)

iv.

x <- -2: 2
p <- 2
testthat::expect_equal(fp(x, p, deriv = TRUE), 2 * x)

v.

x <- -2: 2
p <- 1
testthat::expect_equal(fp(x, p, deriv = TRUE), sign(x))

vi.

x <- numeric(length=0)

p <- 2
testthat::expect_equal(length(fp(x, p)), length(abs(x)^p))
testthat::expect_equal(length(fp(x, p, deriv = TRUE)), length(2 * x))

p <- 1
testthat::expect_equal(length(fp(x, p)), length(abs(x)^p))
testthat::expect_equal(length(fp(x, p, deriv = TRUE)), length(sign(x)))

c.

testthat::test_that("fp works", {
  ### i.
  x <- -2: 2
  testthat::expect_equal(length(fp(x, p = 1)), length(x))
  
  ### ii.
  x <- -2: 2
  p <- 2
  testthat::expect_equal(fp(x, p), abs(x)^p)
  
  ### iii.
  x <- -2: 2
  p <- 1
  testthat::expect_equal(fp(x, p), abs(x)^p)
  
  ### iv.
  x <- -2: 2
  p <- 2
  testthat::expect_equal(fp(x, p, deriv = TRUE), 2 * x)
  
  ### v.
  x <- -2: 2
  p <- 1
  testthat::expect_equal(fp(x, p, deriv = TRUE), sign(x))
  
  ### vi.
  x <- numeric(length=0)
  
  p <- 2
  testthat::expect_equal(length(fp(x, p)), length(abs(x)^p))
  testthat::expect_equal(length(fp(x, p, deriv = TRUE)), length(2 * x))
  
  p <- 1
  testthat::expect_equal(length(fp(x, p)), length(abs(x)^p))
  testthat::expect_equal(length(fp(x, p, deriv = TRUE)), length(sign(x)))
})

Test passed

5.

### While loop 1 (Recommended)
X <- cars$speed

### Step 0
x <- 0
k <- 1
tol <- 1e-5
Kmax <- 1000

g.prime <- function(x, X){
  mean(2 * (x - X))
}

### Step 1 + 2
while ((abs(g.prime(x, X)) >= tol) && (k <= Kmax)){
  x <- x - 0.99^k * g.prime(x, X)
  k <- k + 1
}

### Output
x

[1] 15.4

### Unit test
testthat::expect_equal(x, mean(X), tolerance = tol)

### While loop 2
X <- cars$speed

### Step 0
x <- 0
k <- 1
tol <- 1e-5
Kmax <- 1000

g.prime <- function(x, X){
  mean(2 * (x - X))
}

### Step 1 + 2
while (TRUE){
  x <- x - 0.99^k * g.prime(x, X)
  k <- k + 1
  
  if ((abs(g.prime(x, X)) < tol) || (k > Kmax)){
    break()
  }
}

### Output
x

[1] 15.4

### Unit test
testthat::expect_equal(x, mean(X), tolerance = tol)

6.

a.

grad.des <- function(X, tol, Kmax){
  ### Step 0
  x <- 0
  k <- 1
  
  g.prime <- function(x, X){
    mean(2 * (x - X))
  }
  
  ### Step 1 + 2
  while ((abs(g.prime(x, X)) >= tol) && (k <= Kmax)){
    x <- x - 0.99^k * g.prime(x, X)
    k <- k + 1
  }
  
  list(est = x, abs.g.prime = abs(g.prime(x, X)), num.iter = k - 1)
}

X <- cars$speed
tol <- 1e-5
Kmax <- 1000

### Output
grad.des(X = X, tol = tol, Kmax = Kmax)

$est
[1] 15.4

$abs.g.prime
[1] 4.678836e-06

$num.iter
[1] 36

### Unit test
testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$est, mean(X), tolerance = tol)

b.

### Unit tests
### Unit test for a different X
X <- -10:5
tol <- 1e-5
Kmax <- 1000
testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$est, mean(X), tolerance = tol)

### Unit test for tolerance
X <- -10:5
tol <- 5
Kmax <- 1000
testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$abs.g.prime <= tol, TRUE)

### Unit test for tolerance (simple way)
### When set tol to be very large, the while loop would not run
X <- -10:5
tol <- 1000000
Kmax <- 1000
testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$est, 0)

### Unit test for Kmax
X <- -10:5
tol <- 1e-5
Kmax <- 3
testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$num.iter, Kmax)

### Unit test for Kmax (simple way)
### When set Kmax=0, the while loop would not run
X <- -10:5
tol <- 1e-5
Kmax <- 0
testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$est, 0)

### Group unit tests
testthat::test_that("grad.des works", {
  ### Unit test for a different X
  X <- -10:5
  tol <- 1e-5
  Kmax <- 1000
  testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$est, mean(X), tolerance = tol)

  ### Unit test for tolerance
  X <- -10:5
  tol <- 5
  Kmax <- 1000
  testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$abs.g.prime <= tol, TRUE)
  
  ### Unit test for tolerance (simple way)
  ### When set tol to be very large, the while loop would not run
  X <- -10:5
  tol <- 1000000
  Kmax <- 1000
  testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$est, 0)
  
  ### Unit test for Kmax
  X <- -10:5
  tol <- 1e-5
  Kmax <- 3
  testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$num.iter, Kmax)
  
  ### Unit test for Kmax (simple way)
  ### When set Kmax=0, the while loop would not run
  X <- -10:5
  tol <- 1e-5
  Kmax <- 0
  testthat::expect_equal(grad.des(X = X, tol = tol, Kmax = Kmax)$est, 0)
})

Test passed

7. (Optional)

(Not graded)

grad.des.2 <- function(X, tol, Kmax, p){
  ### Step 0
  x <- 0
  k <- 1
  n <- length(X)
  
  ### Only needs to modify the g.prime
  g.prime <- function(x, X){
    1 / n * sum(fp(x = x - X, p = p, deriv = TRUE))
  }
  
  ### Step 1 + 2
  while ((abs(g.prime(x, X)) >= tol) && (k <= Kmax)){
    x <- x - 0.99^k * g.prime(x, X)
    k <- k + 1
  }
  
  list(est = x, abs.g.prime = abs(g.prime(x, X)), num.iter = k - 1)
}

dat <- c(1, 2, 3, 5, 1000)
tol <- 1e-5
Kmax <- 1000

### Unit tests
testthat::expect_equal(grad.des.2(X = dat, tol = tol, Kmax = Kmax, p = 1)$est, median(dat), tolerance = tol)
testthat::expect_equal(grad.des.2(X = dat, tol = tol, Kmax = Kmax, p = 2)$est, mean(dat), tolerance = tol)

8.

### Gradient descent algorithm
###
### Input:
###   X:            A data vector
###   tol:          A small value for the tolerance level
###   Kmax:         A large value for maximum number of iterations
###   p:            Order of absolute value in the target function
###
### Target function: g (x) = mean((abs(x - X)) ^ p)
###
### Returns a list contains the following:
###   est:          A numeric value that minimizes the target function
###   abs.g.prime:  A numeric value equals to the absolute value of the derivative of the target function
###   num.iter:     The number of iterations runs
###
### Example usages:
###   p = 1:        Approximate the median of X
###   p = 2:        Approximate the mean of X

Homework 5

STAT 486/586 – Introduction to Statistical Computing

Yuhang (Tom) LIN

Due: March 1, 2022

1.

a.

i.

ii.

b.

i.

ii.

c.

d.

2. (Optional)

3.

4.

a.

b.

i.

ii.

iii.

iv.

v.

vi.

c.

5.

6.

a.

b.

7. (Optional)

8.