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)