STAT5003 - Q3.1

# Load necessary libraries
library(ggplot2)
library(rmarkdown)

Simulate 100 Observations

# Simulate 100 observations from a standard normal distribution
set.seed(42)
n <- 100
data <- rnorm(n)

Determine the Best Bandwidth

# Define the true density function (density of standard normal distribution)
true_density <- function(x) {
  dnorm(x)
}

# Define a function to calculate the ISE
calculate_ise <- function(h) {
  # Estimate density
  density_estimate <- density(data, bw = h, n = 512, from = -3, to = 3)
  f_hat <- approx(density_estimate$x, density_estimate$y, xout = density_estimate$x)$y
  
  # Calculate the ISE
  ise <- sum((f_hat - true_density(density_estimate$x))^2) * diff(density_estimate$x[1:2])
  return(ise)
}

# Create a grid of bandwidth values
bandwidths <- seq(0.1, 1.5, length.out = 100)
# Calculate ISE for each bandwidth
ise_values <- sapply(bandwidths, calculate_ise)

# Find the bandwidth with the minimum ISE
optimal_bandwidth <- bandwidths[which.min(ise_values)]
optimal_bandwidth

## [1] 0.3686869

Plot the Results

# Plot the relationship between bandwidth and ISE
ggplot(data = data.frame(bandwidth = bandwidths, ISE = ise_values), aes(x = bandwidth, y = ISE)) +
  geom_line() +
  geom_point(data = data.frame(bandwidth = optimal_bandwidth, ISE = min(ise_values)),
             aes(x = bandwidth, y = ISE), color = "red", size = 3) +
  ggtitle("Bandwidth vs. Integrated Square Error") +
  xlab("Bandwidth") +
  ylab("Integrated Square Error")

# Plot the density estimate with the optimal bandwidth against the true density
density_optimal <- density(data, bw = optimal_bandwidth, n = 512, from = -3, to = 3)
true_density_vals <- true_density(density_optimal$x)

plot(density_optimal, main = "Optimal Bandwidth Density Estimate vs. True Density", xlab = "x", ylab = "Density")
lines(density_optimal$x, true_density_vals, col = "red", lwd = 2)
legend("topright", legend = c("Estimated Density", "True Density"), col = c("black", "red"), lwd = 2)