library(ggplot2)
library(tidyr)
lcg <- function(X_0, a, c, m, length) {
  stream <- vector("numeric", length)
  stream[1] <- X_0
  
  for (i in 1:(length)) {
    stream[i + 1] <- (a * stream[i] + c) %% m
  }
  
  stream <- stream[-1]
  stream <- stream / m
  
  return(stream)
}

Banks et al. (2010). Discrete-Event System Simulation (5th Ed.).

Chapter 8, Problem 5:

Inverse-transform technique

\(F(x) = \begin{cases} 0, & x \leq -3 \\ 1/2 + x/6, & -3 < x \leq 0 \\ 1/2 + x^2 / 32, & 0 < x \leq 4 \\ 1, & x > 4 \end{cases}\)

\(f(x) = \begin{cases} 0, & x \leq -3 \\ 1/6, & -3 < x \leq 0 \\ x/16, & 0 < x \leq 4 \\ 0, & x > 4 \end{cases}\)

\(F(X) = R\)

\(1/2 + X/6 = R,\ -3 < X \leq 0\)

\(X/6 = R - 1/2,\ -3 < X \leq 0\)

\(X = 6R - 3,\ -3 < X \leq 0\)

\(1/2 + X^2 / 32 = R,\ 0 < X \leq 4\)

\(X^2 / 32 = R - 1/2,\ 0 < X \leq 4\)

\(X^2 = 32R - 16,\ 0 < X \leq 4\)

\(X = 4 \sqrt{2R - 1},\ 0 < X \leq 4\)

\(X = \begin{cases} 6R - 3, & 0 \leq R \leq 1/2 \\ 4 \sqrt{2R - 1}, & 1/2 < R\leq 1 \end{cases}\)

rvg_it <- function(n) {
  # for testing
  # R <- runif(1000)

  # use random seed
  R <- lcg(sample(2^31 - 1, 1), 7^5, 0, 2^31 - 1, 1000)
    x <- sapply(R, function(r) 
    if (r >= 0 & r <= 0.5) {
      6 * r - 3 
    } else if (r > 1/2 & r <= 1) {
      4 * sqrt(2 * r - 1)
    }
  )
  return(x)
}

x <- rvg_it(1000)

pdf <- function (x) {
  out <- sapply(x, function(el) 
    if (el <= -3) {
      0
    } else if (el > -3 & el <= 0) {
      1 / 6
    } else if (el > 0 & el <= 4) {
      el / 16
    } else {
      0
    }) 
  return(out)
}

ggplot(data.frame(x = x), aes(x)) +
  geom_histogram(aes(y = ..density..), bins = 50) +
  stat_function(fun = pdf, colour = "red")

times <- data.frame(method = c(rep("Inverse-Transform", 1000), rep("Acceptance-Rejection", 1000)),
                    runtime = numeric(2000))

for(i in 1:1000) {
  t0 <- Sys.time()
  rvg_it(1000)
  t1 <- Sys.time()
  times$runtime[i] <- t1 - t0
}

Acceptance-rejection technique

rvg_ar <- function(n) {
  # random number stream continously and uniformly distributed on [-3, 4]
  R <- lcg(sample(2^31 - 1, 1), 7^5, 0, 2^31 - 1, n * 10) * 7 - 3
  
  # random number stream continously and uniformly distributed on [0, 1]
  y <- lcg(sample(2^31 - 1, 1), 7^5, 0, 2^31 - 1, n * 10)
  
  # get densities of values in R 
  p <- sapply(R, function(r) 
    if (r <= -3) {
      0
    } else if (r > -3 & r <= 0) {
      1 / 6
    } else if (r > 0 & r <= 4) {
      r / 16
    } else {
      0
    }
  )
  x <- vector("numeric", n)
  j <- 1
  for (i in 1:length(R)) {
    if(y[i] < p[i]) {
      x[j] <- R[i]
      j <- j + 1
    }
    if (j == n + 1) {
      break
    }
  }
  return(x)
}

x <- rvg_ar(1000)

ggplot(data.frame(x = x), aes(x)) +
  geom_histogram(aes(y = ..density..), bins = 50) +
  stat_function(fun = pdf, colour = "red")

for(i in 1001:2000) {
  t0 <- Sys.time()
  rvg_ar(1000)
  t1 <- Sys.time()
  times$runtime[i] <- t1 - t0
}

Runtime comparison

ggplot(times, aes(runtime, group = method, color = method)) +
  geom_line(stat = "density")