STA506_Assignment 1: Estimating CDF and PDF

Question 1: Cumulative Distribution Function (CDF) Estimation

Part A

Since we are estimating an empirical distribution, we can apply the formula here: \[ F_n(x) =\frac{1}{n}\sum_{i=1}^n I(X_i \le x) = \frac{[\ \text{Number of } X_i \le x]}{n} \]

In our case, n = 8. So, for every point, we should get an increase of our cumuluative probability to be 1/8.

x = c(23, 45, 67, 89, 112, 156, 189, 245)
unique.x <- sort(unique(x)) # sort unique values of x

myECDF <- function(indat, outx){
  # outx - a vector of given values
  freq.table <- table(indat)                          # frequency table
  uniq <- as.numeric(names(freq.table))          # unique data values
  rep.x <- as.vector(freq.table)                   # frequencies of the unique data values
  cum.rel.feq <- cumsum(rep.x)/sum(rep.x)       # cumulative relative frequencies: CDF
  cum.prob <- NULL
  for (i in 1:length(outx)){
    intvl.id <- which(uniq <= outx[i])      # identify the index meeting the condition
    cum.prob[i] <- cum.rel.feq[max(intvl.id)] # extract the cumulative prob according to CDF 
  }
  cum.prob          
}

myECDF(x, unique.x)

[1] 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000

We can also verify our estimates via the ecdf() function:

empirical_cdfx <- ecdf(x)
empirical_cdfx(x)

[1] 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000

The results seem to match, so our estimate seems to be correct.

Part B

Yes, I agree with the colleague. 4 data points are less than or equal to 100. Referencing the formula from Part A, we can say the numerator is 4, and our n is still 8, so our cumulative probability should be 4/8, or 1/2. To verify, we can use the ecdf() function from Part A:

empirical_cdfx(100)

[1] 0.5

And we get an output of 0.5. So the colleague is correct.

Question 2: Density Function Estimation

Part A

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

hist.y <-hist(y, breaks = seq(12, 26, length.out = 4), freq = FALSE, main = "Histogram of Failure Times")

hist.y$density

[1] 0.06428571 0.08571429 0.06428571

#seq function taken as suggestion to force 3 bins from stackoverflow.com

According to our graph, the estimates are 0.0643, 0.0857, and 0.0643. The histogram looks symmetrical about the second bin, but we need more bins/data to be sure of the distribution’s true shape.

Part B

Function used to estimate the densities:

mykerf <- function(in.data, h, out.x){
  n <- length(in.data)     # sample size
  den <- NULL              # density vector to store output values
  for (i in 1:length(out.x)){
    den[i] <- sum(dnorm(out.x[i], mean=in.data, sd = h))/n  # kernel density formula
  }
  den  # return density values based on the out.x values
}

These are the values of the gaussian kernel estimates with h = 2:

d1 <- mykerf(y, 2, y)
d1

 [1] 0.03855284 0.06671954 0.07046656 0.07600381 0.07587074 0.07360217
 [7] 0.06993562 0.06400803 0.05652506 0.04237326

To verify, these are the kernel estimates using the density() function:

d2 <- density(y, bw = 2, n = 10)
d2$x #the estimated values of the distribution

 [1]  6.300000  9.088889 11.877778 14.666667 17.455556 20.244444 23.033333
 [8] 25.822222 28.611111 31.400000

d2$y #the values of the corresponding density estimates

 [1] 0.0002260196 0.0060884538 0.0330404547 0.0664300753 0.0762536983
 [6] 0.0711297655 0.0611758565 0.0370947939 0.0068889038 0.0002403998

Part C

d4 <- density(y, bw = 2, kernel = "epanechnikov")

Part D

The different methods use different weighting functions. So different values aren’t weighted the same.

d5 <- density(y, bw = 1.5, kernel = "gaussian", n = 10)
d6 <- density(y, bw = 2.5, kernel = "gaussian", n = 10)
d7 <- density(y, bw = 1.5, kernel = "epanechnikov", n = 10)
d8 <- density(y, bw = 2.5, kernel = "epanechnikov", n = 10)

It seems that the estimated values of the distributions didn’t change, but their density estimates changed.

d5

Call:
    density.default(x = y, bw = 1.5, kernel = "gaussian", n = 10)

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

       x               y           
 Min.   : 7.80   Min.   :0.000297  
 1st Qu.:13.32   1st Qu.:0.011274  
 Median :18.85   Median :0.047917  
 Mean   :18.85   Mean   :0.040735  
 3rd Qu.:24.38   3rd Qu.:0.070350  
 Max.   :29.90   Max.   :0.076519  

d6

Call:
    density.default(x = y, bw = 2.5, kernel = "gaussian", n = 10)

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

       x               y            
 Min.   : 4.80   Min.   :0.0001868  
 1st Qu.:11.82   1st Qu.:0.0042516  
 Median :18.85   Median :0.0262266  
 Mean   :18.85   Mean   :0.0320283  
 3rd Qu.:25.88   3rd Qu.:0.0579508  
 Max.   :32.90   Max.   :0.0740767  

d7

Call:
    density.default(x = y, bw = 1.5, kernel = "epanechnikov", n = 10)

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

       x               y          
 Min.   : 7.80   Min.   :0.00000  
 1st Qu.:13.32   1st Qu.:0.01405  
 Median :18.85   Median :0.04828  
 Mean   :18.85   Mean   :0.04132  
 3rd Qu.:24.38   3rd Qu.:0.06972  
 Max.   :29.90   Max.   :0.07981  

d8

Call:
    density.default(x = y, bw = 2.5, kernel = "epanechnikov", n = 10)

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

       x               y           
 Min.   : 4.80   Min.   :0.000000  
 1st Qu.:11.82   1st Qu.:0.005187  
 Median :18.85   Median :0.027290  
 Mean   :18.85   Mean   :0.032055  
 3rd Qu.:25.88   3rd Qu.:0.056422  
 Max.   :32.90   Max.   :0.074101