绘制经验分布函数以及生存函数

关键词：经验分布函数

经验分布函数：

函数ecdf(x)中的参数x要求是unique（唯一、不重复）值。在这个例子中，y值unique值数量为91。该函数返回的是描述文件中所述的类“ecdf”的函数，简言之这个函数返回的是一个函数。所以，使用返回的函数计算给定数据的经验分布函数(ecdf)。生存函数与分布函数具有如下关系：

S(x)=1-F(x)。使用这个关系计算经验生存函数(empirical survival function)

值得注意的是：在plot函数中为经验生存函数调用的函数stepfun需要一个比向量x（stepfun函数的第一个参数）多一个值的向量y（stepfun函数的第二个参数）.所以在生存函数前加了一个1。这被认为是合理的，因为对于比有序样本（stepfun函数的第二个参数，y）的第一个值小的值，经验生存函数应该是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) 
  }
}

alpha <- 0.2807
lambda <- 0.0030

导入数据y:

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)
length(unique(y)) #unique y的数量为91
length(y) ##y的长度 109

计算给定参数的理论经验分布函数和理论生存函数：

# theoretical cdf and survival function for given parameters: 

cdf_th <- cdf(y = y, alpha = alpha, lambda = lambda)
S_th <- S(y = y, alpha = alpha, lambda = lambda)

验证一下：通过上述自定义函数计算得到的S_th和cdf_th之间是否具有S(x)=1-F(x)这样的关系。即判断S(x)是否等于1-F(x)。

# 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: 

cdf_emp_fun <- ecdf(y) #This  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

其中plot函数中，参数xlab和ylab分表表示图形中横轴和纵轴的标签。lines函数表示在已有的图形上绘制曲线。

绘制经验生存函数和理论生存函数

# 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

 这里stepfun函数中参数x和参数y分别表示该函数中第一个参数，以及第二个参数。

 

参考：data visualization - How can i find empirical survival function using survival function in R? - Cross Validated