Assignment 1: Estimating CDF and PDF

Assignment Objectives

• Develop a clear technical understanding of nonparametric cumulative distribution function (CDF) estimation and various kernel density estimators.

• Translate mathematical formulas into R functions and apply them to solve related problems.

• Create effective visualizations to demonstrate your understanding of key concepts in the following questions.

Question 1: Cumulative Distribution Function (CDF) Estimation

The following failure times (in hours) were observed for 8 electronic components:

23, 45, 67, 89, 112, 156, 189, 245

1. Write an R function implementing the ECDF \(\hat{F}_n(t)\) according to its mathematical definition. Validate your implementation using R’s ecdf() function on the given data, with comparison based on their step functions.

x <- c(23, 45, 67, 89, 112, 156, 189, 245)

my_ecdf_f <- function(x){
  n <- length(x)
  function(t){
    vapply(t, function(tt) sum(x <= tt) / n, numeric(1))
  }
}

F_r <- ecdf(x)
F_mine <- my_ecdf_f(x)

t_grid <- seq(min(x) - 10, max(x) + 10, length.out = 1000)

plot_df <- data.frame(
  t = t_grid,
  R_ecdf = F_r(t_grid),
  My_ecdf = F_mine(t_grid)
)

plot_df_long <- plot_df %>%
  pivot_longer(
    cols = c(R_ecdf, My_ecdf),
    names_to = "Type",
    values_to = "CDF"
  )

Fn.plt <- ggplot(plot_df_long, aes(x = t, y = CDF, color = Type)) +
  geom_step(linewidth = 1) +
  scale_color_manual(values = c("R_ecdf" = "blue", "My_ecdf" = "red")) +
  labs(
    title = "Empirical CDF: My Implementation vs R ecdf()",
    subtitle = paste("Sample size n =", length(x)),
    x = "t",
    y = "F̂(t)"
  ) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.margin = margin(t = 35, r = 20, b = 30, l = 30, unit = "pt")
  )

ggplotly(Fn.plt)

The ECDF computed using my custom ECDF implementation closely matches the estimate produced by R’s ecdf() function, indicating that the implementation is correct.

2. A colleague claims that the probability of failure before 100 hours is 0.5 based on these data. Do you agree? Explain your reasoning using the empirical cumulative distribution function (ECDF).

I agree with the claim because the ECDF evaluated at 100 hours is 0.5, meaning that about 50% of the observed failures occurred before 100 hours. Therefore, the ECDF supports this claim.

Question 2: Density Function Estimation

Consider the following failure times from a mechanical system:

12.3, 14.7, 15.2, 16.8, 18.1, 19.4, 20.6, 22.3, 23.9, 25.4

1. Create a histogram of the data using 3 equally spaced bins. What is the estimated density in each bin? Describe the shape of the histogram’s distribution.

x <- c(12.3, 14.7, 15.2, 16.8, 18.1, 19.4, 20.6, 22.3, 23.9, 25.4)
h <- hist(x, breaks = seq(min(x), max(x), length.out = 4), freq = FALSE)

The estimated density of each bin in order is 0.0687023, 0.0916031, 0.0687023. The distribution is roughly symmetrical and unimodal.

2. Write an R function that computes kernel density estimates using a Gaussian kernel with \(h=2\). Validate your implementation against R’s built-in density() function.

\[ \hat{f}_h(t) = \frac{1}{nh}\sum_{i=1}^n K\left( \frac{t-t_i}{h}\right), \ \ \text{ where } \ \ K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2/2}. \]

my_kde <- function(x,h){
  n <- length(x)
  function(t){
    vapply(t, function(xi) 1/(n*h) * sum(dnorm((xi - x)/h)), numeric(1))
  }
}

F_mine <- my_kde(x, h = 2)
F_r <- density(x, bw = 2)

t_grid <- F_r$x

plot_df <- data.frame(
  t = t_grid,
  R_density = F_r$y,
  My_density = F_mine(t_grid)
)

plot_df_long <- plot_df %>%
  pivot_longer(
    cols = c(R_density, My_density),
    names_to = "Type",
    values_to = "Density"
  )

kde_plt <- ggplot(plot_df_long, aes(x = t, y = Density, color = Type)) +
  geom_line(linewidth = 1) +
  labs(
    title = "KDE: My Implementation vs R density()",
    subtitle = paste("Gaussian kernel, h =", 2, ", n =", length(x)),
    x = "t",
    y = "Density"
  ) +
  theme(plot.title = element_text(hjust = 0.5))

ggplotly(kde_plt)

The KDE computed using my custom Gaussian kernel implementation closely matches the estimate produced by R’s density() function, indicating that the implementation is correct.

3. Write a custom R function that computes kernel density estimates using the Epanechnikov kernel with \(h=2\). Validate your implementation by comparing results with R’s built-in density() function for Gaussian kernel estimation.

\[ \hat{f}_h(t) = \frac{1}{nh}\sum_{i=1}^n K\left( \frac{t-t_i}{h}\right), \ \ \text{ where } \ \ K(u) = \frac{3}{4}(1 - u^2) \ \ \text{ for } \ \ |u| \le 1. \]

my_kde <- function(x,h){
  n <- length(x)
  function(t){
    vapply(t, function(tt) {
      u <- (tt - x)/h
      K <- ifelse(abs(u) <= 1, (3/4)*(1 - u^2), 0)
      1/(n*h) * sum(K)}, numeric(1))
  }
}

F_mine <- my_kde(x, h = 2)
F_r <- density(x, bw = 2, kernel = "epanechnikov")

t_grid <- F_r$x

plot_df <- data.frame(
  t = t_grid,
  R_density = F_r$y,
  My_density = F_mine(t_grid)
)

plot_df_long <- plot_df %>%
  pivot_longer(
    cols = c(R_density, My_density),
    names_to = "Type",
    values_to = "Density"
  )

kde_plt <- ggplot(plot_df_long, aes(x = t, y = Density, color = Type)) +
  geom_line(linewidth = 1) +
  labs(
    title = "KDE: My Implementation vs R density()",
    subtitle = paste("Gaussian kernel, h =", 2, ", n =", length(x)),
    x = "t",
    y = "Density"
  ) +
  theme(plot.title = element_text(hjust = 0.5))

ggplotly(kde_plt)

The KDE computed using my custom Epanechnikov kernel implementation closely matches the estimate produced by R’s density() function, indicating that the implementation is correct. However, the two curves don’t completely overlap this time, I’m assuming the density() function must go above and beyond just the mathematical equation.

4. How does the choice of kernel (Gaussian vs. Epanechnikov) affect the density estimate? For both kernel estimators applied to this dataset, what happens when we select \(h=1.5\) versus \(h=2.5\)?

The Guassian curve is much more smooth than the Epanechnikov curve. Epanechnikov densities are allowed to be zero at the tails, this doesn’t seem to be true for Gaussian densities. For both kernels, the density curves become more rigid when h is decreased to 1.5 and more smooth when its increased to 2.5.