Q1) Theory

Q2) Application of Inv. CDF

(a-b)

(c) Simulate

set.seed(123)
lambda <- 3
eta <- 1
n <- 10000

# Real Exponential Distribution
real_samples <- eta + rexp(n, rate = lambda)
real_density <- density(real_samples)
plot(real_density,
     main = "Density Plot of Exponential(3,1)",
     col = "red")

# Inverse CDF Samples
F_inv <- function(u) {
  eta - log(1 - u) / lambda
}

u <- runif(n, min = 0, max = 1)
x <- F_inv(u)

invCDF_density <- density(x)
lines(invCDF_density, col = "blue")

legend("topright",
       legend = c("Real Density", "Inverse CDF"),
       col = c("red", "blue"),
       lty = c(1, 1))

  • As we can see our inverse CDF method does a pretty good job of est. the Real Density

(d) Quantiles

# Theoretical probabilities
p <- (1:n - 0.5) / n

# Theoretical quantiles
theoretical_q <- eta - log(1 - p) / lambda

# Q-Q Plot
plot(theoretical_q,
     sort(x),
     xlab = "Theoretical Quantiles",
     ylab = "Sample Quantiles",
     main = "Q-Q Plot for Two-Parameter Exponential")

abline(0, 1, col = "red")

  • as we can see, there is minor deviation from our diagnol line on our quantile plot – indicating a good fit

Q3) Application of Inv. CDF

(a-b)

(c) Simulate

set.seed(123)
alpha <- 4
beta <- 1
n <- 10000

# Inverse CDF Samples
F_inv <- function(u) {
  beta * (u)^(-1 / alpha)
}

u <- runif(n)
x <- F_inv(u)

# Estimated density of simulated samples
invCDF_density <- density(x)
plot(invCDF_density,
     main = "Density Plot of Pareto(4,1)",
     col = "blue")

# Theoretical Pareto density
curve(
  ifelse(x >= beta,
         alpha * beta^alpha / x^(alpha + 1),
         0),
  from = beta,
  to = max(x),
  add = TRUE,
  col = "red",
  lwd = 2
)

legend("topright",
       legend = c("Inverse CDF", "Theoretical Density"),
       col = c("blue", "red"),
       lty = c(1, 1))

  • we see values below 1 because of the density() function, this is not an actual artifact of our generated-data :
range(x)
## [1]  1.000015 11.122712
  • ie. there arent actually values below 1, it just looks like that.

(d) QQ-PLT

# Theoretical probabilities
p <- (1:n - 0.5) / n

# Theoretical quantiles
theoretical_q <- beta * (1 - p)^(-1 / alpha)

# Q-Q Plot
plot(theoretical_q,
     sort(x),
     xlab = "Theoretical Quantiles",
     ylab = "Sample Quantiles",
     main = "Q-Q Plot for Pareto(4,1)")

abline(0, 1, col = "red")

  • as we can see, its a strong fit with the exception of one outlier