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)