piecewise exp firing time sampling

pwexp.R

rpwexp <- function(n, failRates) {
  end_time <- cumsum(failRates$duration)
  n_rate <- nrow(failRates)
  t <- 0
  tau <- Tk <- rep(0, n)
  Pk <- rexp(n = n)
  notfired <- rep(TRUE, n)
  for (i in 1:n_rate) {
    delta_tk <- (Pk[notfired] - Tk[notfired]) / failRates$rate[i]
    tau[notfired] <- t + delta_tk
    if (i < n_rate) {
      notfired[notfired] <- tau[notfired] >= end_time[i]
      t <- t + failRates$duration[i]
      Tk <- Tk + failRates$rate[i]*failRates$duration[i]
    }
    if (!any(notfired)) {
      break
    }
  }
  return(tau)
}

failRates <- data.frame(rate = c(1, 3, 0.5, 0.1, 2, 1.5, 0.1, 0.9), duration = c(.5, .5, 1, 0.25, 0.5, 1, 2, 1))

x <- rpwexp(n = 1e4, failRates = failRates)
plot(sort(x), (1e4:1)/10001,
     log = "y", main = "PW Exponential simulated survival curve",
     xlab = "Time", ylab = "P{Survival}", col = adjustcolor("red",alpha.f = 0.15), pch = 16)