How can i find empirical survival function using survival function in R?

Question

The survival function is given by:

S(y; α, λ) = (α/α−1)* (1 − α^(−e^(−λy ))), if α is not equal 1

   =   e^(−λy) if α =1

y = 1 4 4 7 11 13 15 15 17 18 19 19 20 20 22 23 28 29 31 32 36 37 47 48 49 50 54 54 55 59 59 61 61 66 72 72 75 78 78 81 93 96 99 108 113 114 120 120 120 123 124 129 131 137 145 151 156 171 176 182 188 189 195 203 208 215 217 217 217 224 228 233 255 271 275 275 275 286 291 312 312 312 315 326 326 329 330 336 338 345 348 354 361 364 369 378 390 457 467 498 517 566 644 745 871 1312 1357 1613 1630

alpha = 0.2807

lambda = 0.0030

I used this coding but it's returning 91 values, so i can't plot it with my data of 109 values.What am i doing wrong here?

alpha=0.2807

lambda=0.0030

cdf <-if(alpha!=1){

((alpha^(1-exp(-lambda*y)))-1)/ (alpha-1)}else{

1-(exp(-lambda*y))} cdf

S <- if(alpha!=1){

(alpha/(alpha-1))(1-alpha^(-exp(-lambday)))}else{

exp(-lambda*y) } S

#The fitted APE survival function

plot(DATA,S,"l")

#The empirical survival function

ecdf(S)

RESULT :

ecdf(S)

Empirical CDF

Call: ecdf(S)

x[1:91] = 0.0037469, 0.0039439, 0.0085506, ..., 0.97909, 0.99472

Answer 1

Inside the function ecdf only unique values of the argument x are used. In your case length(unique(y)) returns 91.

I've attached a reproducible example.

What I've done: The function ecdf returns a function of class 'ecdf' as it is stated in the description file. So I've used the returned function to compute the empirical cdf of the given data. The survival function is the counterpart of the cdf, i.e. for a random variable X $$ S(x) = 1 - F(x) .$$ I've used this relationship to calculate the empirical survival function.

At the end of the reprex I've created two visualizations to compare the theoretical and empirical distributions.

One more note: The function stepfun which is called inside the plot for the empirical survival function needs a vector y that has one more value than the vector x. So I put a 1 in front of the survival function. I think that is reasonable, because for values of y smaller than the first value of the ordered sample the empirical survival function should be 1.

# functions: 
cdf <- function(y, alpha,lambda) {
  if (alpha != 1) {
    apecdf <-((alpha^(1-exp(-lambda * y))) - 1)/ (alpha - 1)
  } else {
    apecdf<- 1-(exp(-lambda * y))
  }
  return(apecdf)
} 

S <- function(y, alpha, lambda) {
  if (alpha != 1) {
    (alpha / (alpha - 1)) * (1 - alpha^(-exp(-lambda * y)))
  } else {
    exp(-lambda * y) 
  }
}
    
# given: 
y <- c(1, 4, 4, 7, 11, 13, 15, 15, 17, 18, 19, 19, 20, 20, 22, 23, 28, 29, 31, 32, 
       36, 37, 47, 48, 49, 50, 54, 54, 55, 59, 59, 61, 61, 66, 72, 72, 75, 78, 78, 
       81, 93, 96, 99, 108, 113, 114, 120, 120, 120, 123, 124, 129, 131, 137, 145, 
       151, 156, 171, 176, 182, 188, 189, 195, 203, 208, 215, 217, 217, 217, 224, 228, 
       233, 255, 271, 275, 275, 275, 286, 291, 312, 312, 312, 315, 326, 326, 329, 330, 
       336, 338, 345, 348, 354, 361, 364, 369, 378, 390, 457, 467, 498, 517, 566, 644, 
       745, 871, 1312, 1357, 1613, 1630)
alpha <- 0.2807
lambda <- 0.0030

# theoretical cdf and survival function for given parameters: 
S_th <- S(y = y, alpha = alpha, lambda = lambda)
cdf_th <- cdf(y = y, alpha = alpha, lambda = lambda)

# implementation check, S(x) = 1 - F(x) where F(.) is the cdf: 
all.equal(S_th, 1 - cdf_th)
#> [1] TRUE

# empirical cdf and survival function: 
## argument x of is specified as a numeric vector of observations; ie. y in your case: 

cdf_emp_fun <- ecdf(x = y) #  is a function
cdf_emp <- cdf_emp_fun(y) # values

# S as the counterpart of the cdf: 
S_emp <- 1 - cdf_emp

# Plotting the cdf: 
plot(cdf_emp_fun, main = "Empirical and Theoretical CDF", xlab = "y", ylab = "F(y)") # empirical cdf
lines(x = y, y = cdf_th, col = "red", type = "l") # theoretical cdf

# Plotting the survival functions: 
plot(stepfun(x = y, y = c(1, S_emp)), main = "Empirical and Theoretical Survival Function", xlab = "y", ylab = "S(y)") # empirical S
lines(x = y, y = S_th, col = "red", type = "l") # theoretical S

^{Created on 2020-07-16 by the reprex package (v0.3.0)}