1 Q1

below, i just try and get an idea of whats going on by visualizing :

set.seed(123)
# experiment 
exp <- rexp(n=1000, rate=1/3)
par(mfrow= c(1,2))
hist(exp)
hist(sqrt(exp)) # Rayleigh distribution
abline(v=mean(sqrt(exp)), col = "red")

1.1 a)

Anyways, we should get some value close to :

mc2 <- 
  function(N, R){
  e <- rexp(n=N, rate=R)
  mean(sqrt(e))
}

est <- mc2(1000, 1/3); est
## [1] 1.531248

1.1.1 est. variance of estimate

m <- 1000
r <- rexp(n=m, rate=1/3) |> sqrt() 
# empirical calculation
V <- (1/(m-1))*sum((r - mean(r))^2) # m-1 for unbiased est. 
theoretical <- 3-(sqrt(3*pi)*(1/2))^2; theoretical
## [1] 0.6438055
V
## [1] 0.6819743
SE <- sqrt(V/m)
CI95 <- est + c(-1,1) * SE * qnorm(0.975); CI95
## [1] 1.480064 1.582432

1.2 b)

truth <- sqrt(3*pi) / 2
B <- 1000
N <- 1000

est <- lower <- upper <- numeric(B) # three numeric vectors, each of length B
cover <- logical(B)

for(i in 1:B){

  r <- sqrt(rexp(N, rate = 1/3))

  est[i] <- mean(r)

  SE <- sqrt(var(r) / N)

  lower[i] <- est[i] - qnorm(0.975) * SE
  upper[i] <- est[i] + qnorm(0.975) * SE

  cover[i] <- (lower[i] <= truth) & (truth <= upper[i])
}

mean(cover)
## [1] 0.958
  • so, yeh – 95% of our ci are in like we expected

1.2.1 Visualize 1000 repeated experiments

plot(est,
     ylim = range(lower, upper),
     pch = 16)

segments(1:B, lower,
         1:B, upper,
         col = ifelse(cover, "black", "red"))

abline(h = truth,
       col = "blue",
       lwd = 2)

2 Q2

2.1 a)

Let the parameter of interest (the center of the standard Cauchy distribution) be

\[ \theta = 0. \]

For a sample of size \(n\), the level-\(k\) trimmed mean is obtained by removing the \(k\) smallest and \(k\) largest observations.

The Monte Carlo algorithm is as follows:

  1. Choose the sample size \(n\), trimming level \(k\), and the number of Monte Carlo replications \(B\).

  2. For each replication \(b=1,2,\ldots,B\):

    • Generate a random sample \[ X_1^{(b)}, X_2^{(b)}, \ldots, X_n^{(b)} \overset{\text{iid}}{\sim} \text{Cauchy}(0,1). \]
    • Sort the sample: \[ X_{(1)}^{(b)} \le X_{(2)}^{(b)} \le \cdots \le X_{(n)}^{(b)}. \]
    • Compute the level-\(k\) trimmed mean: \[ T_k^{(b)} = \frac{1}{n-2k} \sum_{i=k+1}^{n-k} X_{(i)}^{(b)}. \]
    • Compute the squared error: \[ \left(T_k^{(b)}-\theta\right)^2 = \left(T_k^{(b)}\right)^2, \] since \(\theta=0\).

  3. Estimate the mean squared error (MSE) by averaging the squared errors over all \(B\) replications:

\[ \widehat{\mathrm{MSE}}(T_k) = \frac{1}{B} \sum_{b=1}^{B} \left(T_k^{(b)}-\theta\right)^2 = \frac{1}{B} \sum_{b=1}^{B} \left(T_k^{(b)}\right)^2. \]

repeat for each trimming level \(k\) to compare the estimated MSEs.

2.2 b)

set.seed(2)

# Simulation settings
n <- 30                  # sample size
B <- 10000               # number of Monte Carlo replications
k_values <- 1:5          # trimming levels
theta <- 0               # center of the standard Cauchy distribution

# Store the estimated MSE for each value of k
mse_estimates <- numeric(length(k_values))

# Monte Carlo simulation
for (j in seq_along(k_values)) {
  
  k <- k_values[j]
  trimmed_means <- numeric(B)
  
  for (b in 1:B) {
    
    # Generate a sample from the standard Cauchy distribution
    x <- rcauchy(n, location = 0, scale = 1)
    
    # Sort the sample
    x_sorted <- sort(x)
    
    # Remove the k smallest and k largest observations
    x_trimmed <- x_sorted[(k + 1):(n - k)]
    
    # Compute the level-k trimmed mean
    trimmed_means[b] <- mean(x_trimmed)
  }
  
  # Estimate the MSE relative to the center theta = 0
  mse_estimates[j] <- mean((trimmed_means - theta)^2)
}

# Display the estimated MSE values
results <- data.frame(
  k = k_values,
  Estimated_MSE = mse_estimates
)

results

2.3 c)

The simulation results are

results

For the fixed sample size \(n=30\), the estimated MSE decreases as the trimming level \(k\) increases. indicating removing more extreme observations makes the estimator more stable, because the Cauchy distribution has very heavy tails.

However, these results alone do not establish consistency. Consistency requires that, for a fixed trimming level \(k\),

\[ T_{k,n} \xrightarrow{P} 0 \qquad \text{as } n \to \infty, \]

or equivalently that the estimation error becomes arbitrarily small as the sample size increases.

For a level-\(k\) trimmed mean, only a fixed number \(k\) of observations is removed from each tail. As \(n\) increases, the proportion of observations removed is

\[ \frac{2k}{n} \longrightarrow 0. \]

Because the Cauchy distribution has extremely heavy tails, the remaining extreme order statistics can still have a substantial effect on the trimmed mean. Therefore, for fixed \(k\), the level-\(k\) trimmed mean is not generally a consistent estimator of the Cauchy center.

Thus, the implementation shows that increasing \(k\) reduces the MSE when \(n=30\), but it does not agree with the claim of consistency for a fixed level \(k\). To investigate consistency directly, the simulation should be repeated for increasing values of \(n\) while holding \(k\) fixed and checking whether the estimated MSE approaches zero.

MSE should approach zero as \(n\) increases

3 Q3

3.0.1 (a)

Let :

\[ \theta=\int_0^2 x e^{-x/4}\,dx. \]

If \(U \sim \operatorname{Uniform}(0,2)\), then

\[ \theta = 2\,E\left[Ue^{-U/4}\right]. \]

Therefore, the crude Monte Carlo estimator is

\[ \widehat{\theta}_{MC} = \frac{2}{N}\sum_{i=1}^{N}U_i e^{-U_i/4}, \]

where \(U_1,\ldots,U_N\) are independent random variables from the \(\operatorname{Uniform}(0,2)\) distribution.

set.seed(123)

# Number of Monte Carlo samples
N <- 10000

# Generate random values from Uniform(0, 2)
u <- runif(N, min = 0, max = 2)

# Evaluate the integrand
g_u <- u * exp(-u / 4)

# Crude Monte Carlo estimate of the integral
theta_hat_mc <- 2 * mean(g_u)

theta_hat_mc
## [1] 1.43893

3.0.2 (b)

Let

\[ g(u)=u e^{-u/4}, \]

where \(U\sim\operatorname{Uniform}(0,2)\). Since \(2-U\) also has a \(\operatorname{Uniform}(0,2)\) distribution, the antithetic pair is

\[ U_i \qquad \text{and} \qquad 2-U_i. \]

The antithetic variate estimator is

\[ \widehat{\theta}_{AV} = \frac{1}{N} \sum_{i=1}^{N} \left[ g(U_i)+g(2-U_i) \right]. \]

Equivalently,

\[ \widehat{\theta}_{AV} = \frac{1}{N} \sum_{i=1}^{N} \left[ U_i e^{-U_i/4} + (2-U_i)e^{-(2-U_i)/4} \right]. \]

set.seed(123)

# Number of antithetic pairs
N <- 10000

# Generate Uniform(0, 2) random variables
u <- runif(N, min = 0, max = 2)

# Construct the antithetic variables
u_antithetic <- 2 - u

# Evaluate the integrand for each pair
g_u <- u * exp(-u / 4)
g_u_antithetic <- u_antithetic * exp(-u_antithetic / 4)

# Antithetic variate Monte Carlo estimate
theta_hat_av <- mean(g_u + g_u_antithetic)

theta_hat_av
## [1] 1.444835

3.0.3 (c)

Let

\[ g(u)=u e^{-u/4}, \]

where \(U\sim\operatorname{Uniform}(0,2)\). Then

\[ \theta=2E[g(U)]. \]

Use \(U\) as the control variate because its expected value is known:

\[ E[U]=1. \]

The control variate estimator is

\[ \widehat{\theta}_{CV} = 2\left[ \overline{g(U)} - \widehat{c}\left(\overline{U}-1\right) \right], \]

where the estimated optimal control coefficient is

\[ \widehat{c} = \frac{\widehat{\operatorname{Cov}}(g(U),U)} {\widehat{\operatorname{Var}}(U)}. \]

set.seed(123)

# Number of Monte Carlo samples
N <- 10000

# Generate Uniform(0, 2) random variables
u <- runif(N, min = 0, max = 2)

# Evaluate the integrand
g_u <- u * exp(-u / 4)

# Known expected value of the control variate
E_u <- 1

# Estimate the optimal control coefficient
c_hat <- cov(g_u, u) / var(u)

# Control variate Monte Carlo estimate
theta_hat_cv <- 2 * mean(g_u - c_hat * (u - E_u))

theta_hat_cv
## [1] 1.444799

3.0.4 (d) Comparison with the Theoretical Value and Estimated Variances

The theoretical value of the integral is

\[ \begin{aligned} \theta &=\int_0^2 xe^{-x/4}\,dx \\ &=\left[-4(x+4)e^{-x/4}\right]_0^2 \\ &=16-24e^{-1/2}. \end{aligned} \]

The estimated variance of each Monte Carlo estimator is computed as the sample variance of its summands divided by \(N\).

# Theoretical value
theta_theoretical <- 16 - 24 * exp(-1/2)

# Recover the crude Monte Carlo estimate from objects in the global environment
theta_hat_mc <- 2 * mean(g_u)

# Estimated variance of the crude Monte Carlo estimator
var_mc <- var(2 * g_u) / N

# Antithetic-pair contributions
av_contributions <- g_u + g_u_antithetic

# Estimated variance of the antithetic variate estimator
var_av <- var(av_contributions) / N

# Control-variate contributions
cv_contributions <- 2 * (g_u - c_hat * (u - E_u))

# Estimated variance of the control variate estimator
var_cv <- var(cv_contributions) / N

# Comparison table
comparison <- data.frame(
  Method = c(
    "Crude Monte Carlo",
    "Antithetic variates",
    "Control variate",
    "Theoretical value"
  ),
  Estimate = c(
    theta_hat_mc,
    theta_hat_av,
    theta_hat_cv,
    theta_theoretical
  ),
  Estimated_Variance = c(
    var_mc,
    var_av,
    var_cv,
    NA
  ),
  Absolute_Error = c(
    abs(theta_hat_mc - theta_theoretical),
    abs(theta_hat_av - theta_theoretical),
    abs(theta_hat_cv - theta_theoretical),
    0
  )
)

comparison

The theoretical value is

theta_theoretical
## [1] 1.443264
comparison[2,]

The antithetic variate and control variate estimators both produced estimates very close to the theoretical value and had substantially smaller estimated variances than the crude Monte Carlo estimator. In this simulation, the control variate method had the smallest estimated variance, making it the most efficient of the three approaches.