Math6630_AS_3_Q6_3.R

# Danny Tshitumbu
# ID: 211964756

# Applied Stats 1: Q. 6. 3.

# 6.3.a.

n <- 1e6
X <- rcauchy(n, 0, 1)   
W1 <- exp(-abs(X)^3 / 3) * (pi * (1 + X^2))
Std_W1 <- W1/sum(W1)
sigma1 <- sum(X^2 * Std_W1)
sigma1

## [1] 0.7760992

# 6.3. b.
rm(list = ls()) # Clear all variables to avoid confusion/mistakes in calculation

n2 <- 1e6
U1 <- runif(n2)
X <- -log(U1)
U2 <- runif(n2)
test <- exp(-abs(X)^3 / 3) / exp((2/3) - X)
X_ret <- X[U2 <= test]
sigma2 <- (1/length(X_ret)) * sum(X_ret^2)
sigma2

## [1] 0.7773671

# 6.3.c.

ord_Xret <- sort(X_ret)
n3 <- length(ord_Xret)-1
diff_X <- rep(0, n3)
for (i in 1:n3) {
  diff_X[i] <- ord_Xret[i+1] - ord_Xret[i]
}
h <- (ord_Xret^2)[1:n3]
q <- (exp(-abs(ord_Xret)^3 / 3))[1:n3]
sigma3 <- sum(diff_X * h * q) / sum(diff_X * q)
sigma3

## [1] 0.7765471

# Performance of b
system.time({
  n2 <- 1e6
  U1 <- runif(n2)
  X <- -log(U1)
  U2 <- runif(n2)
  test <- exp(-abs(X)^3 / 3) / exp((2/3) - X)
  X_ret <- X[U2 <= test]
  sigma2 <- (1/length(X_ret)) * sum(X_ret^2)
})

##    user  system elapsed 
##   0.383   0.027   0.414

# Perfomance of c
system.time({
  ord_Xret <- sort(X_ret)
  n3 <- length(ord_Xret)-1
  diff_X <- rep(0, n3)
  for (i in 1:n3) {
    diff_X[i] <- ord_Xret[i+1] - ord_Xret[i]
  }
  h <- (ord_Xret^2)[1:n3]
  q <- (exp(-abs(ord_Xret)^3 / 3))[1:n3]
  sigma <- sum(diff_X * h * q) / sum(diff_X * q)
})

##    user  system elapsed 
##   2.022   0.023   2.069

# Comment: We note that the approach in (c) takes longer, so the estimation using
# rejection sampling seems to be faster than the one described by Philippe and Robert