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.

# Failure times (in hours) were observed for 8 electronic components:
x <- c(23, 45, 67, 89, 112, 156, 189, 245)
n <- length(x)

# Implementation of ECDF
ecdf_manual <- function(t, x) {
  mean(x <= t)
}

# Validation of Implementation
ecdf_r <- ecdf(x)

t_vals <- seq(0, 260, by = 1)

plot(t_vals, sapply(t_vals, ecdf_manual, x = x),
     type = "s", col = "blue", lwd = 2, main = "Empirical CDF of Failure Times",ylab = "ECDF", xlab = "Failure Time (hours)")
lines(ecdf_r, col = "red", lwd = 2, lty = 2)
legend("bottomright", legend = c("Manual ECDF", "R's ecdf()"),
       col = c("blue", "red"), lty = c(1,2))

Both ECDF’s return identical values and is continuous from the right at every point and increases by 1/n = 0.125 for each observed failure time

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

# Observation of colleague claim using ECDF 
ecdf_manual(100, x)

[1] 0.5

The ECDF at 100 hours is 0.5, meaning 50% of observed failures occurred before 100 hours.

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.

# Failure times from a mechanical system
y <- c(12.3, 14.7, 15.2, 16.8, 18.1,
       19.4, 20.6, 22.3, 23.9, 25.4)

n <= length(y)

[1] TRUE

# Compute range and bin width
data_range <- range(y)
bin_width <- diff(data_range) / 3

# Define bin boundaries
breaks <- seq(from = data_range[1],
              to   = data_range[2],
              length.out = 4)

breaks

[1] 12.30000 16.66667 21.03333 25.40000

# Count observations in each bin
bin_counts <- table(cut(y, breaks = breaks, include.lowest = TRUE))
bin_counts

[12.3,16.7]   (16.7,21]   (21,25.4] 
          3           4           3 

# Compute estimated density in each bin
bin_densities <- as.numeric(bin_counts) / (n * bin_width)

bin_densities

[1] 0.08587786 0.11450382 0.08587786

# Present results
density_table <- data.frame(
  Bin = levels(cut(y, breaks = breaks, include.lowest = TRUE)),
  Count = as.numeric(bin_counts),
  Density = bin_densities
)

density_table

          Bin Count    Density
1 [12.3,16.7]     3 0.08587786
2   (16.7,21]     4 0.11450382
3   (21,25.4]     3 0.08587786

# Visualization 
hist(y, breaks = breaks, probability = TRUE,
     main = "Histogram with Estimated Bin Densities",
     xlab = "Failure Time")

The estimated density in each histogram bin is computed by dividing the number of observations in that bin by the product of the sample size and the bin width. The densities sum to 1 when multiplied by band width with the middle bin having the highest density where failure times are most concentrated.

The histogram is unimodal and symmetric, having a smooth distribution

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}. \]

# Gaussian Kernel Density Estimation (h = 2)
kde_gaussian <- function(t, x, h) {
  n <- length(x)
  sapply(t, function(tt) {
    mean(dnorm((tt - x) / h)) / h
  })
}

# Validation against density()
t_grid <- seq(min(y) - 2, max(y) + 2, length = 200)

plot(t_grid, kde_gaussian(t_grid, y, h = 2),
     type = "l", lwd = 2, col = "blue",
     xlab = "t", ylab = "Density",
     main = "Gaussian KDE Comparison")

lines(density(y, bw = 2), col = "red", lty = 2)

legend("topright",
       legend = c("Manual KDE", "R density()"),
       col = c("blue", "red"), lty = c(1, 2))

The curves overlap correctly and confirms the correctness of the manuel KDE implementation

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. \]

# Epanechnikov Kernel Density Estimation (h = 2)
kde_epanechnikov <- function(t, x, h) {
  n <- length(x)
  sapply(t, function(tt) {
    u <- (tt - x) / h
    mean(0.75 * (1 - u^2) * (abs(u) <= 1)) / h
  })
}

# Visualization 
plot(t_grid, kde_epanechnikov(t_grid, y, h = 2),
     type = "l", lwd = 2, col = "green",
     main = "Epanechnikov KDE (h = 2)",
     ylab = "Density", xlab = "Failure Time")
lines(density(y, bw = 2), col = "red", lty = 2)
legend("topright", c("Epanechnikov", "Gaussian"),
       col = c("green", "red"), lty = c(1,1))

The Epanechnikov KDE is more compact with limited support

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

Kernel choice

• Gaussian: smoother tails, infinite support
• Epanechnikov: compact support, sharper cutoff
• Differences are small

density(y, bw = 1.5)

Call:
    density.default(x = y, bw = 1.5)

Data: y (10 obs.);  Bandwidth 'bw' = 1.5

       x               y           
 Min.   : 7.80   Min.   :0.000297  
 1st Qu.:13.32   1st Qu.:0.015141  
 Median :18.85   Median :0.055866  
 Mean   :18.85   Mean   :0.045149  
 3rd Qu.:24.38   3rd Qu.:0.071778  
 Max.   :29.90   Max.   :0.077948  

density(y, bw = 2.5)

Call:
    density.default(x = y, bw = 2.5)

Data: y (10 obs.);  Bandwidth 'bw' = 2.5

       x               y            
 Min.   : 4.80   Min.   :0.0001868  
 1st Qu.:11.82   1st Qu.:0.0057773  
 Median :18.85   Median :0.0343512  
 Mean   :18.85   Mean   :0.0355077  
 3rd Qu.:25.88   3rd Qu.:0.0644896  
 Max.   :32.90   Max.   :0.0743308  

Bandwidth effect

• h = 1.5 → less smoothing, more bumps, higher variance
• h = 2.5 → more smoothing, flatter curve, higher bias

Conclusion

• Bandwidth has a much larger impact on the density estimate than kernel choice.