#### Ordered statistics - Time to extinction with the Dodo dodosightingtimes <- c(1662, 1638, 1631, 1628, 1628, 1611, 1607, 1602, 1601, 1598) compute.lambda <- function(dates, v){ n <- length(dates) lambda <- matrix(data=NA, nrow=n, ncol=n) for(i in 1:n){ for(j in 1:n){ lambda[i,j] <- (gamma(2*v + i) * gamma(v + j))/(gamma(v + i) * gamma(j)) } } return(lambda) } ## optimal linear estimation function OLE <- function(data, alpha){ ## sort the data and define a parameter for length (here k) obs <- rev(sort(data)) k <- length(obs) ## --------------------------------------------------------------- ## compute v, e, lambda, and a ## --------------------------------------------------------------- ## estimate the shape parameter of the joint Weibull distribution ## 4th equation in Roberts and Solow 2003 v <- (1/(k-1)) * sum(log((obs[1] - obs[k])/(obs[1] - obs[2:(k-1)]))) ## define a vector of k 1’s e <- matrix(rep(1,k), ncol=1) ## Λ is the symmetric k x k matrix with typical element lambda <- compute.lambda(obs, v) ## make the Λ matrix symmetric lambda <- ifelse(lower.tri(lambda), lambda, t(lambda)) ## vector of weights is given by: a = (e^t * Λ^-1 * e)^-1 * Λ^-1*e a <- as.vector(solve(t(e) %*% solve(lambda) %*% e)) * solve(lambda) %*% e ## --------------------------------------------------------------- ## calculate confidence intervals ## --------------------------------------------------------------- SL <- (-log(1-alpha/2)/length(obs))^-v SU <- (-log(alpha/2)/length(obs))^-v lowerCI <- max(obs) + ((max(obs)-min(obs))/(SL-1)) ## lower confidence interval upperCI <- max(obs) + ((max(obs)-min(obs))/(SU-1)) ## upper confidence interval ## --------------------------------------------------------------- ## calculates the time of extinction ## --------------------------------------------------------------- extincttime <- sum(t(a)%*%obs) ## --------------------------------------------------------------- ## put the results into a data from ## --------------------------------------------------------------- res <- data.frame(Estimate=extincttime, lowerCI=lowerCI, upperCI=upperCI) ##return the results return(res) } OLE(data=dodosightingtimes, alpha=0.05) q()