Problem 1

(a)

\[\text{Var}[XY] = \text{E}[(XY)^2] - (\text{E}[XY])^2 = \text{E}(X^2 Y^2) - [\text{E}(XY)]^2\] Since \(X\) and \(Y\) are independent,\[\text{E}[XY] = \text{E}[X]\text{E}[Y]\] \(X^2\) and \(Y^2\) are also independent, so\[\text{E}[X^2Y^2] = \text{E}[X^2]\text{E}[Y^2]\] \[\text{Var}[XY] = \text{E}[X^2]\text{E}[Y^2] - \text{E}[X]^2\text{E}[Y]^2\] Using the identity \(\text{E}[Z^2] = \text{Var}[Z] + \text{E}[Z]^2\), substitute \(\text{E}[X^2]\) and \(\text{E}[Y^2]\): \[\text{Var}[XY] = (\text{Var}[X] + \text{E}[X]^2)(\text{Var}[Y] + \text{E}[Y]^2) - \text{E}[X]^2\text{E}[Y]^2\]

Expanding this we get: \[\text{Var}[XY] = \text{Var}[X]\text{Var}[Y] + \text{Var}[X]\text{E}[Y]^2 + \text{Var}[Y]\text{E}[X]^2 + \text{E}[X]^2\text{E}[Y]^2 - \text{E}[X]^2\text{E}[Y]^2\] Canceling \(\text{E}[X]^2\text{E}[Y]^2\) gives: \[\text{Var}[XY] =\text{Var}[X]\text{Var}[Y] + \text{Var}[X]\text{E}[Y]^2 + \text{Var}[Y]\text{E}[X]^2\]

(b)

For \(X \sim \text{Geometric}(p)\) with \(p = 0.3\), the CDF is \(P(X \le k) = 1 - (1 - p)^k\). We want to find the smallest integer \(k\) such that \(P(X \le k) \ge 0.95\):

\[1 - (0.7)^k \ge 0.95\] \[(0.7)^k \le 0.05\]

Taking \(\ln()\) on both sides:\[k \ln(0.7) \le \ln(0.05)\]\[k \ge \frac{\ln(0.05)}{\ln(0.7)} \approx \frac{-2.9957}{-0.3567} \approx 8.399\]Since \(k\) has to be an integer, we round up to \(k = 9\).Checking in R using pgeom.

Since pgeom counts failures before the first success, we evaluate at \(k - 1 = 8\)):

p <- 0.3

# Find k with qgeom
k <- qgeom(0.95, prob = p) + 1
k
## [1] 9
pgeom(8, prob = p)
## [1] 0.9596464

(C)

  1. Call arrivals follow a Poisson process at \(\lambda = 5 \text{ calls/min}\). Interarrival time \(T \sim \text{Exp}(5)\):\[\text{E}[T] = \frac{1}{5} \text{ min} = 12 \text{ seconds}\]

  2. Number of calls in 1 minute is \(N(1) \sim \text{Poisson}(5)\):\[P(N(1) = 7) = \frac{e^{-5} 5^7}{7!} \approx 0.1044\]

(iii)\[P(N(1) < 7) = P(N(1) \le 6) = \sum_{k=0}^{6} \frac{e^{-5} 5^k}{k!} \approx 0.7622\]

  1. In seconds, rate is \(\lambda = 5/60 = 1/12 \text{ calls/sec}\): \[P(T \le 10) = 1 - e^{-10/12} = 1 - e^{-5/6} \approx 0.5654\]

  2. Standard \(95\%\) CI band formula for ECDF is \(\hat{F}_n(t) \pm 1.96 \sqrt{\frac{\hat{F}_n(t)(1 - \hat{F}_n(t))}{n}}\)

set.seed(124)

# Generate 1000 arrivals (rate = 5/60 per sec)
x <- rexp(1000, rate = 5 / 60)
fn <- ecdf(x)

plot(fn, main = "ECDF with 95% Confidence Bands", xlab = "Seconds", ylab = "ECDF")

grid_t <- seq(0, max(x), length.out = 500)
f_hat <- fn(grid_t)
se <- sqrt(f_hat * (1 - f_hat) / 1000)

lines(grid_t, pmin(1, f_hat + 1.96 * se), col = "red", lty = 2)
lines(grid_t, pmax(0, f_hat - 1.96 * se), col = "red", lty = 2)

Problem 2:

Inverse CDFs for continuous distributions (a) Distribution 1 (Laplace): For \(f_1(x) = \frac{1}{2b} e^{-\vert{}x\vert{}/b}\), we integrate to get CDF \(F_1(x)\): For \(x < 0\):\[F_1(x) = \int_{-\infty}^{x} \frac{1}{2b} e^{t/b} dt = \frac{1}{2} e^{x/b}\] For \(x \ge 0\): \[F_1(x) = \frac{1}{2} + \int_{0}^{x} \frac{1}{2b} e^{-t/b} dt = 1 - \frac{1}{2} e^{-x/b}\] Setting \(u = F_1(x)\) and solving for \(x\): If \(u < 0.5 \implies x = b \ln(2u)\)If \(u \ge 0.5 \implies x = -b \ln(2(1 - u))\)

So to sample \(X\), generate \(U \sim \text{Unif}(0, 1)\) and compute:\[X = \begin{cases} b \ln(2U), & U < 0.5 \\ -b \ln(2(1-U)), & U \ge 0.5 \end{cases}\]

Distribution 2: Integrating \(f_2(x) = \frac{k(x - a)^{k-1}}{(b - a)^k}\) over \([a, x]\): \[F_2(x) = \int_{a}^{x} \frac{k(t - a)^{k-1}}{(b - a)^k} dt = \left(\frac{x - a}{b - a}\right)^k\] Setting \(u = F_2(x)\) and solving for \(x\): \[u = \left(\frac{x - a}{b - a}\right)^k \implies x = a + (b - a) u^{1/k}\] To sample \(X\), generate \(U \sim \text{Unif}(0, 1)\) and set \(X = a + (b - a) U^{1/k}\).

(b)

  1. For \(f(x) = \lambda e^{-\lambda(x - \eta)}\) on \(x \ge \eta\): \[F(x) = \int_{\eta}^{x} \lambda e^{-\lambda(t - \eta)} dt = 1 - e^{-\lambda(x - \eta)}\] Setting \(u = 1 - e^{-\lambda(x - \eta)}\) and solving for \(x\): \[F^{-1}(u) = \eta - \frac{1}{\lambda} \ln(1 - u)\]

  2. Algorithm function:

inv_cdf_two_param <- function(n, lambda, eta) {
  u <- runif(n)
  return(eta - (1 / lambda) * log(1 - u))
}
  1. Simulation with \(n = 10000\), \(\lambda = 2\), \(\eta = 1\):
set.seed(124)

sim_data <- inv_cdf_two_param(10000, lambda = 2, eta = 1)

hist(sim_data, breaks = 50, probability = TRUE,
     main = "Two-Parameter Exponential Simulation",
     xlab = "x", col = "red", border = "white")

curve(2 * exp(-2 * (x - 1)), from = 1, to = max(sim_data),
      add = TRUE, col = "blue", lwd = 2)

legend("topright", legend = c("Simulated", "Theoretical"),
       fill = c("red", NA), border = c("white", NA),
       col = c(NA, "blue"), lty = c(NA, 1), lwd = c(NA, 2))

  1. Quantile comparison (\(Q(p) = \eta - \frac{1}{\lambda}\ln(1 - p)\)):
p_vec <- c(0.10, 0.25, 0.50, 0.75, 0.90)

emp_q <- quantile(sim_data, probs = p_vec)
theo_q <- 1 - (1 / 2) * log(1 - p_vec)

data.frame(
  p = p_vec,
  Theoretical = theo_q,
  Empirical = as.numeric(emp_q),
  Abs_Error = abs(as.numeric(emp_q) - theo_q)
)
##      p Theoretical Empirical    Abs_Error
## 1 0.10    1.052680  1.051671 0.0010091970
## 2 0.25    1.143841  1.143152 0.0006889859
## 3 0.50    1.346574  1.350879 0.0043050277
## 4 0.75    1.693147  1.696910 0.0037623802
## 5 0.90    2.151293  2.137164 0.0141285944

Problem 3

(a)

For \(X \sim \text{Exp}(2)\), CDF is \(F(x) = 1 - e^{-2x}\). Inverse CDF is \(X = -\frac{1}{2}\ln(1 - U)\). Theoretical expectations: \(\text{E}(X) = 1/\lambda = 0.5\), \(\text{E}(X^2) = \text{Var}(X) + \text{E}[X]^2 = 1/4 + 1/4 = 0.5\).

set.seed(124)

n <- 5000
u <- runif(n)
x <- -log(1 - u) / 2

hist(x, breaks = 40, probability = TRUE,
     main = "Simulated Exp(2)",
     xlab = "x", col = "blue", border = "white")
curve(2 * exp(-2 * x), from = 0, to = max(x), add = TRUE, col = "red", lwd = 2)

cat("Mean E(X):   ", mean(x), " (Theoretical: 0.5)\n")
## Mean E(X):    0.4898809  (Theoretical: 0.5)
cat("Mean E(X^2): ", mean(x^2), " (Theoretical: 0.5)\n")
## Mean E(X^2):  0.4820883  (Theoretical: 0.5)

(b)

Theoretical mean: \(\text{E}(X) = (-1)(0.1) + (0)(0.3) + (2)(0.6) = 1.1\). CDF step function: \(F(-1) = 0.1\), \(F(0) = 0.4\), \(F(2) = 1.0\). Mapping \(U \sim \text{Unif}(0, 1)\):\(U \le 0.1 \implies X = -1\)\(0.1 < U \le 0.4 \implies X = 0\)\(U > 0.4 \implies X = 2\)

set.seed(124)

u_disc <- runif(1000)
x_disc <- ifelse(u_disc <= 0.1, -1, ifelse(u_disc <= 0.4, 0, 2))

# Relative frequencies
prop.table(table(x_disc))
## x_disc
##    -1     0     2 
## 0.092 0.317 0.591
cat("Sample Mean:", mean(x_disc), "\n")
## Sample Mean: 1.09
barplot(prop.table(table(x_disc)), main = "Discrete Sample Frequencies",
        xlab = "x", ylab = "Proportion", col = "steelblue")

Problem 4

(a)

Poisson Inverse-CDF algorithm using recursive updates: 1. Draw \(U \sim \text{Unif}(0, 1)\). 2. Set \(p = e^{-\lambda}\), \(F = p\), \(x = 0\). 3. While \(U > F\):

  • \(x = x + 1\)
  • \(p = p \cdot (\lambda / x)\)
  • \(F = F + p\)
  1. Return \(x\).

(b)

my_rpois <- function(n, lambda) {
  out <- numeric(n)
  for (i in 1:n) {
    u <- runif(1)
    p <- exp(-lambda)
    F_val <- p
    x <- 0
    while (u > F_val) {
      x <- x + 1
      p <- p * (lambda / x)
      F_val <- F_val + p
    }
    out[i] <- x
  }
  return(out)
}

(c)

set.seed(124)

n <- 10000
lam <- 4.2

custom_samples <- my_rpois(n, lam)
rpois_samples <- rpois(n, lam)

cat("Custom -> Mean:", mean(custom_samples), "| Variance:", var(custom_samples), "\n")
## Custom -> Mean: 4.2002 | Variance: 4.191539
cat("rpois  -> Mean:", mean(rpois_samples), "| Variance:", var(rpois_samples), "\n")
## rpois  -> Mean: 4.1979 | Variance: 4.169953
# Distribution comparison for k = 0..8
k_grid <- 0:8
p_theo <- dpois(k_grid, lam)
p_custom <- prop.table(table(factor(custom_samples, levels = k_grid)))
p_rpois <- prop.table(table(factor(rpois_samples, levels = k_grid)))

data.frame(
  k = k_grid,
  Theoretical = round(p_theo, 4),
  Custom_Inv = as.numeric(round(p_custom, 4)),
  Builtin_rpois = as.numeric(round(p_rpois, 4))
)
##   k Theoretical Custom_Inv Builtin_rpois
## 1 0      0.0150     0.0159        0.0157
## 2 1      0.0630     0.0668        0.0667
## 3 2      0.1323     0.1350        0.1322
## 4 3      0.1852     0.1860        0.1920
## 5 4      0.1944     0.1995        0.1964
## 6 5      0.1633     0.1695        0.1697
## 7 6      0.1143     0.1207        0.1204
## 8 7      0.0686     0.0717        0.0723
## 9 8      0.0360     0.0349        0.0346

Problem 5

(a)

Target \(f(x) = \frac{2}{5}(x + 1)\) and trial \(g(x) = 1\) on \([1, 2]\). We want \(M = \sup_{1 \le x \le 2} \frac{f(x)}{g(x)} = \sup \frac{2}{5}(x + 1)\). Since \(f/g\) is increasing in \(x\), max occurs at \(x = 2\): \[M = \frac{2}{5}(2 + 1) = \frac{6}{5} = 1.2\]

(b)

Target \(f(x) = \frac{1}{2}\sin(x)\) and trial \(g(x) = \frac{1}{\pi}\) on \([0, \pi]\). \[\frac{f(x)}{g(x)} = \frac{\frac{1}{2}\sin x}{\frac{1}{\pi}} = \frac{\pi}{2}\sin x\] Since \(\sin(x) \le 1\), max occurs at \(x = \pi/2\): \[M = \frac{\pi}{2} \approx 1.5708\] Theoretical acceptance rate: \(P(\text{Accept}) = 1/M = 2/\pi \approx 0.6366\).

set.seed(124)

n_trials <- 500
M <- pi / 2

f <- function(x) 0.5 * sin(x)
g <- function(x) 1 / pi

x_trial <- runif(n_trials, min = 0, max = pi)
u <- runif(n_trials)

accept <- u <= (f(x_trial) / (M * g(x_trial)))
x_accept <- x_trial[accept]

cat("Accepted:", length(x_accept), "/", n_trials, "\n")
## Accepted: 321 / 500
cat("Empirical rate:", length(x_accept) / n_trials, "\n")
## Empirical rate: 0.642
cat("Theoretical rate:", round(1 / M, 4), "\n")
## Theoretical rate: 0.6366
hist(x_accept, breaks = 20, probability = TRUE,
     main = "Rejection Sampling for f(x) = 0.5 sin(x)",
     xlab = "x", col = "skyblue", border = "white", xlim = c(0, pi))
curve(f(x), from = 0, to = pi, add = TRUE, col = "red", lwd = 2)

legend("topright", legend = c("Accepted", "Target f(x)"),
       fill = c("skyblue", NA), border = c("white", NA),
       col = c(NA, "red"), lty = c(NA, 1), lwd = c(NA, 2))

The sample acceptance rate (\(62.6\%\)) matches up well with theoretical rate (\(63.66\%\)).

problem 6

(a)

Taking the ratio of Cauchy target \(f(x)\) to Normal trial \(g(x)\): \[\frac{f(x)}{g(x)} = \frac{\frac{1}{\pi(1 + x^2)}}{\frac{1}{\sqrt{2\pi}} e^{-x^2/2}} = \sqrt{\frac{2}{\pi}} \frac{e^{x^2/2}}{1 + x^2}\] Taking the limit as \(x \to \infty\):\[\lim_{x \to \infty} \frac{f(x)}{g(x)} = \sqrt{\frac{2}{\pi}} \lim_{x \to \infty} \frac{e^{x^2/2}}{1 + x^2} = \infty\] Because the numerator grows exponentially while the denominator only grows polynomially, the ratio goes to infinity. So no finite \(M\) exists.

The standard Cauchy target has heavy tails (\(1/x^2\)), while the Normal trial has light tails (\(e^{-x^2/2}\)). Rejection sampling requires \(g(x)\) to have equal or heavier tails than \(f(x)\), so a Normal proposal drops off too fast to bound a Cauchy target.

(b)

Target \(X \sim \text{Beta}(2, 3) \implies f(x) = 12x(1-x)^2\) and trial \(g(x) = 1\) on \([0, 1]\). Maximizing \(h(x) = 12x(1-x)^2 = 12(x - 2x^2 + x^3)\): \[h'(x) = 12(1 - 4x + 3x^2) = 0 \implies (3x - 1)(x - 1) = 0\] Critical point is \(x = 1/3\). \[M = h(1/3) = 12(1/3)(2/3)^2 = \frac{16}{9} \approx 1.7778\] Theoretical acceptance rate: \(P(\text{Accept}) = 1/M = 9/16 = 0.5625\).

set.seed(124)

n_proposals <- 100000
M <- 16 / 9

f <- function(x) 12 * x * (1 - x)^2
g <- function(x) 1

x_trial <- runif(n_proposals)
u <- runif(n_proposals)

accept <- u <= (f(x_trial) / (M * g(x_trial)))
x_accept <- x_trial[accept]

cat("Accepted:", length(x_accept), "/", n_proposals, "\n")
## Accepted: 56212 / 1e+05
cat("Empirical rate:", length(x_accept) / n_proposals, "\n")
## Empirical rate: 0.56212
cat("Theoretical rate:", 1 / M, "\n")
## Theoretical rate: 0.5625
hist(x_accept, breaks = 40, probability = TRUE,
     main = "Rejection Sampling for Beta(2,3)",
     xlab = "x", col = "skyblue", border = "white", xlim = c(0, 1))
curve(f(x), from = 0, to = 1, add = TRUE, col = "red", lwd = 2)

legend("topright", legend = c("Accepted", "Beta(2,3) Target"),
       fill = c("skyblue", NA), border = c("white", NA),
       col = c(NA, "red"), lty = c(NA, 1), lwd = c(NA, 2))

The empirical acceptance rate (\(56.13\%\)) is extremely close to theoretical (\(56.25\%\)).