S_a <- function(a, v, Id, t)   # V = Y - A   t <- Y
{
  Tt <- cbind(time = t, address = c(1: length(t)))
  Aa <- cbind(time = a, address = rep(0, length(a)))
  Vv <- cbind(time = v, address = rep(0, length(v)))
  u <- rbind(Tt, Aa, Vv)
  u <- u[order(u[,1]), ]
  Q <- colSums(outer(a, u[,1], function(x,y) as.numeric(x <= y)) + Id * outer(v, u[,1], function(x,y) as.numeric(x <= y)))
  dQ <- Q - c(0, Q[-length(Q)])
  K <- colSums(outer(a, u[,1], function(x,y) as.numeric(x >= y)) + outer(v, u[,1], function(x,y) as.numeric(x >= y)))
  Sa_avt <- cumprod(1 - dQ /K)          #value at points a_i, v_i and t_i
  u <- cbind(u, Sa_avt)[order(u[,2]),]
  Sa_hat <- u[which(u[,2] > 0),3]            # Sa_hat is not ordered from 0 to 1, but with the same order of the input Y
  return(Sa_hat)
}

#估计F
F_t<- function(a, y, Id, t, is.order)    ####   is.order = 1 represent that Fhat is order from 0 to 1
{
  n <- length(a)
  Tt <- cbind(time = t, address = c(1: length(t)))
  Yy <- cbind(time = y, address = rep(0, length(y)))
  u <- rbind(Tt, Yy)
  u <- u[order(u[,1]), ]
  R <- colSums(outer(y, u[,1], function(x,y) as.numeric(x >= y ))) - n * S_a(a, y-a, Id, u[,1])
  N <- colSums(Id * outer(y, u[,1], function(x,y) as.numeric(x <= y )))
  dN <- N - c(0, N[-length(N)])
  F_ty <- 1 - cumprod(1 - dN/ R)
  if(is.order == 1)
  {
    u <- cbind(u, F_ty)
    F_t <- u[which(u[,2] > 0),3]
  }
  else
  {
    u <- cbind(u, F_ty)[order(u[,2]),]
    F_t <- u[which(u[,2] > 0),3]
  }
  return(F_t)
}

这是我分布函数，我想要求使F大于等于p的最小的t。我用for循环让t一点一点的加然后计算，效率太低了