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.

times <- c(23, 45, 67, 89, 112, 156, 189, 245)
uniq.time <- sort(unique(times))  

my.ECDF <- function(indat, outx){
  
  freq.table <- table(indat)                         
  uniq <- as.numeric(names(freq.table))          
  rep.time <- as.vector(freq.table)                   
  cum.rel.feq <- cumsum(rep.time)/sum(rep.time)      
  cum.prob <- NULL
  for (i in 1:length(outx)){
    intvl.id <- which(uniq <= outx[i])      
    cum.prob[i] <- cum.rel.feq[max(intvl.id)] 
  }
  cum.prob             
}

my.ECDF_def <- function(indat, t){
  n <- length(indat)
  sum(indat <= t)/n
}

Fn <- ecdf(times)
grid <- seq(min(times)-10, max(times)+10, by = 1)

plot(Fn, verticals = TRUE, do.points = TRUE,
     main = "ECDF comparison: built-in ecdf() vs my.ECDF()",
     xlab = "Failure time (hours)", ylab = "Empirical CDF")


lines(grid, my.ECDF(indat = times, outx = grid), type = "s", lty = 2, lwd = 2)

legend("bottomright",
       legend = c("ecdf()", "my.ECDF()"),
       lty = c(1,2), lwd = c(1,2), bty = "n")

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).

my.ECDF_def(times, 100)

[1] 0.5

Fn(100)

[1] 0.5

Solution: By using the empirical distribution function, the estimated probability of failure before 100 hours is equal to 0.5. Therefore, I agree with my colleague that the probability of failure before 100 hours is 0.5 based on these small sample data.

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.

fail <- c(12.3, 14.7, 15.2, 16.8, 18.1, 19.4, 20.6, 22.3, 23.9, 25.4)

hisfail<-hist(fail, breaks = 3, freq = FALSE,
      main = "Histogram with 3 equally spaced bins",
      xlab = "Failure time")

hisfail$density

[1] 0.04 0.08 0.06 0.02

Solution:The histogram uses 3 equal-width bins, and the estimated densities are the bin heights shown in the table above.

Solution:Using three equally spaced bins, the histogram suggests an unimodal distribution and approximately symmetric to mildly right-skewed.

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.kerf.gauss <- function(in.data, h, out.x){
  n <- length(in.data)
  den <- numeric(length(out.x))
  for (i in 1:length(out.x)){
    den[i] <- sum(dnorm(out.x[i], mean = in.data, sd = h))/n
  }
  den
}

xx <- seq(min(fail)-3, max(fail)+3, length = 300)

my.den <- my.kerf.gauss(in.data = fail, h = 2, out.x = xx)

kde <- density(fail, bw = 2, kernel = "gaussian")
kde.y <- approx(kde$x, kde$y, xout = xx)$y

plot(xx, my.den, type = "l", lwd = 2,col= "blue",
     main = "Gaussian KDE (h=2): my.kerf.gauss vs density()",
     xlab = "Failure time", ylab = "Density")
lines(xx, kde.y, lty = 2, lwd = 2, col="yellow")
legend("topright", c("my.kerf.gauss()", "density()"),
       lty = c(1,2), lwd = 2, bty = "n")

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.kerf.epan <- function(in.data, h, out.x){
  n <- length(in.data)
  den <- numeric(length(out.x))
  for (i in 1:length(out.x)){
    u <- (out.x[i] - in.data)/h
    K <- ifelse(abs(u) <= 1, 0.75*(1 - u^2), 0)
    den[i] <- sum(K)/(n*h)
  }
  den
}

den.g2 <- my.kerf.gauss(fail, h = 2, out.x = xx)
den.e2 <- my.kerf.epan(fail,  h = 2, out.x = xx)

plot(xx, den.g2, type="l", lwd=2,col="gold",
     main="Kernel comparison (h=2)",
     xlab="Failure time", ylab="Density")
lines(xx, den.e2, lty=2, lwd=2)
legend("topright", c("Gaussian", "Epanechnikov"),
       lty=c(1,2), lwd=2, bty="n")

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\)?

den.g15 <- my.kerf.gauss(fail, h = 1.5, out.x = xx)
den.g25 <- my.kerf.gauss(fail, h = 2.5, out.x = xx)

plot(xx, den.g15, type="l", lwd=2,col="navy",
     main="Gaussian KDE: h=1.5 vs h=2.5",
     xlab="Failure time", ylab="Density")
lines(xx, den.g25, lty=2, lwd=2)
legend("topright", c("h=1.5", "h=2.5"),
       lty=c(1,2), lwd=2, bty="n")

den.e15 <- my.kerf.epan(fail, h = 1.5, out.x = xx)
den.e25 <- my.kerf.epan(fail, h = 2.5, out.x = xx)

plot(xx, den.e15, type="l", lwd=2, col="green",
     main="Epanechnikov KDE: h=1.5 vs h=2.5",
     xlab="Failure time", ylab="Density")
lines(xx, den.e25, lty=2, lwd=2)
legend("topright", c("h=1.5", "h=2.5"),
       lty=c(1,2), lwd=2, bty="n")

Solution:Effect of Kernel Choice (Gaussian vs. Epanechnikov), With the same bandwidth, the Gaussian kernel gives a smoother curve with longer tails, while the Epanechnikov kernel focuses more on nearby observations because it has a limited range. In practice, the difference between kernels is usually small compared to the effect of the bandwidth.

Solution:The bandwidth controls how smooth the density estimate is. A smaller bandwidth (ℎ=1.5) produces a less smooth curve with more local fluctuations, while a larger bandwidth (ℎ=2.5) gives a smoother curve but may hide some details of the data.