source("http://myweb.uiowa.edu/pbreheny/7210/f19/notes/fun.R")

# Fully observed
t <- c(0.1, 0.5, 0.5, 1.6, 2.7)
N <- 99
x <- seq(0.01, 3, len=N)
L1 <- numeric(N)
for (i in 1:N) {
  L1[i] <- prod(dexp(t, rate=x[i]))
}
plotL(x, L1)
abline(h=.15, col='gray75', lwd=2, lty=2)

# Censored at 1
L2 <- L1
for (i in 1:N) {
  L2[i] <- prod(dexp(t[1:3], rate=x[i]))*pexp(1, rate=x[i], low=FALSE)^2
}
L <- cbind(L1, L2)
colnames(L) <- c("Fully observed", "Censored")
plotL(x, L)
abline(h=.15, col='gray75', lwd=2, lty=2)

# Supplemental: Measure in days
L2 <- L1 <- numeric(N)
for (i in 1:N) {
  L1[i] <- prod(dexp(t*365, rate=x[i]/365))
  L2[i] <- prod(dexp(t[1:3]*365, rate=x[i]/365))*pexp(365, rate=x[i]/365, low=FALSE)^2
}
L <- cbind(L1, L2)
colnames(L) <- c("Fully observed", "Censored")
plotL(x, L)
matplot(x, L, col=pal(2), type='l', lty=1, lwd=3, bty='n', xlab=expression(lambda))

# Supplemental: Measure in centuries
L2 <- L1 <- numeric(N)
for (i in 1:N) {
  L1[i] <- prod(dexp(t/100, rate=x[i]*100))
  L2[i] <- prod(dexp(t[1:3]/100, rate=x[i]*100))*pexp(1/100, rate=x[i]*100, low=FALSE)^2
}
L <- cbind(L1, L2)
colnames(L) <- c("Fully observed", "Censored")
plotL(x, L)
matplot(x, L, col=pal(2), type='l', lty=1, lwd=3, bty='n', xlab=expression(lambda))

# Left censored at 0.75
L3 <- L1
for (i in 1:N) {
  L3[i] <- pexp(0.75, rate=x[i])^3 * prod(dexp(t[4:5], rate=x[i]))
}
L <- cbind(L1, L3, L2)
colnames(L) <- c("None", "Left", "Right")
plotL(x, L)
abline(h=.15, col='gray75', lwd=2, lty=2)

# Interval censored #1
L4 <- L1
for (i in 1:N) {
  L4[i] <- (pexp(1, x[i]) - pexp(0, x[i]))^3 *
    (pexp(2, x[i]) - pexp(1, x[i])) *
    (pexp(3, x[i]) - pexp(2, x[i]))
}
L <- cbind(L1, L4, L2)
colnames(L) <- c("None", "Interval", "Right")
plotL(x, L)
abline(h=.15, col='gray75', lwd=2, lty=2)

# Doubly censored
L5 <- L1
for (i in 1:N) {
  L5[i] <- (pexp(1, x[i]) - pexp(0, x[i]))^3 * pexp(1, x[i], low=FALSE)^2
}
L <- cbind(L1, L5, L2)
colnames(L) <- c("None", "Double", "Right")
plotL(x, L)
abline(h=.15, col='gray75', lwd=2, lty=2)

# Left truncation
L6 <- L7 <- L1
for (i in 1:N) {
  L6[i] <- prod(dexp(t[4:5], x[i])/pexp(1, x[i], low=FALSE))
  L7[i] <- prod(dexp(t[4:5], x[i]))
}
L <- cbind(L1, L7, L6)
colnames(L) <- c("Complete", "Ignore", "Truncated")
plotL(x, L[,c(1,3)])
plotL(x, L)

# Right truncation
x <- seq(0.1, 4, len=N)
L7 <- L8 <- L9 <- L1
for (i in 1:N) {
  L7[i] <- prod(dexp(t, rate=x[i]))
  L8[i] <- prod(dexp(t[1:3], x[i])/pexp(1, x[i]))
  L9[i] <- prod(dexp(t[1:3], x[i]))
}
L <- cbind(L7, L9, L8)
colnames(L) <- c("Complete", "Ignore", "Truncated")
plotL(x, L)