set.seed(123)
n <- 10000

### exponential distribution ####
# Generate data
T2 <- rexp(n, rate=2)

# Using the property that for any continuous random variable, applying its survival function 
   # to itself produces a uniform(0, 1) random variable. 
   # Let T ~ Exp(2)
   # Let Y = S_t(T)
   # Then Y ~ Uniform(0, 1)
   # So we can generate Y from a uniform 0, 1 distribution, and apply S_t inverse to it 
   # to get an exponential random variable
Y <- runif(n, min=0, max=1)
T <- -log(Y)/2



# Plot histograms
par(mfrow=c(2,1))
hist(T, 
   breaks=n/100, 
   probability = TRUE)
hist(T2, 
   breaks=n/100, 
   probability = TRUE)

library(lattice)
histogram(~ T + T2, 
   type = "percent",
   xlab = "")


# remember Y is the survival function of T applied to T: S_t(T)
plot(T, Y,
   main = "Survival function for exponentially distributed times with mean 1/2") # plots each time with the associated probability of surviving past that time.

# Plot the cumulative hazard function
H <- -log(Y)
plot(T, H,
   main = "Cumulative hazard function for exponentially distributed times with mean = 1/2")


plot(ecdf(T),
   main = "Empirical cdf of a random variable with Exp(2) dist")
plot(T, 1- Y, main = "Empirical cdf of rv with dist Exp(2)")

Calculate the sample median
median(T)

# Theoretical median
log(2)/2

# Illustrate the memorylessness property of the exponential distribution
t0 <- 0.5
x3 <- T2[T2>t0] - t0
hist(T2, breaks=100)
hist(x3, breaks=100)

### Weibull distribution ###
x <- rweibull(n, shape=2, scale=1)

S <- 1-cumsum(rep(1,n))/n                # empirical survival function

H <- -log(S)
plot(log(sort(x)), log(H), type="o", xlab="ln(t)", ylab="ln(H)")
plot(log(x), log(H), type="o", xlab="ln(t)", ylab="ln(H)")
plot(sort(x), S)

# CDFs
plot(sort(x), 1 - S)
plot(ecdf(x), main = "Empirical cdf")

### log logistic regression ###
alpha <- 1.5
lamda <- 0.01
sigma <- 1/alpha
mu <- -log(lamda)*sigma
y <- rlogis(n, location=mu, scale=sigma)
x <- exp(y)
sum(x>50)/n
mean(x)

## gompertz distribution ##
theta <- 0.01
alpha <- 0.25
y <- runif(n, min=0, max=1)
x <- log(1-log(y)*alpha/theta)/alpha
sum(x>=12)/n
median(x)