vreme = c(200, 450, 600, 700, 800, 1000, 1000, 1000, 450, 600)
status = c(1, 1, 0, 1, 1, 0, 0, 0, 1, 1)

# Kaplan-Majerova ocena funkcije prezivljavanja u tacki t, kad ima ponavljanja u uzroku
KM_ocena = function(t, x, delta) {
  s = 1 
  
  # sortiramo podatke po vremenu
  sortirano = order(x)
  x = x[sortirano]
  delta = delta[sortirano]
  
  # jedinstvene vrednosti vremena
  vreme_unique = unique(x)
  
  for (xj in vreme_unique) {
    if (xj <= t) {
      d_j = sum(delta[x == xj]) # broj dogadjaja u trenutku xj
      R_j = sum(x >= xj) # broj onih pod rizikom pre trenutka xj
      s = s * (1 - d_j / R_j)
    } else {
      break 
    }
  }
  
  return(s)
}

# vektorizujemo funkciju
KM_ocena_vec = Vectorize(KM_ocena, vectorize.args = "t")

# racunamo ocenu za svaku tacku iz uzorka
KM_vrednosti = KM_ocena_vec(vreme, vreme, status)

# ili 
# sapply(vreme, KM_ocena, x = vreme, delta = status)

sortirano = order(vreme)
vreme = vreme[sortirano]
KM_vrednosti = KM_vrednosti[sortirano]

par(mfrow = c(1, 2))
plot(vreme, KM_vrednosti, type = "p", ylim = c(0, 1), xlab = "Vreme", 
     ylab = "Kaplan-Majerova ocena", col = "blue") 
plot(vreme, KM_vrednosti, type = "s", ylim = c(0, 1), xlab = "Vreme", 
     ylab = "Kaplan-Majerova ocena", col = "blue") 