dodoextinctiontime.R

#### 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
	
	## ---------------------------------------------------------------
	## compute time of extinction
	## ---------------------------------------------------------------
	extincttime <- sum(t(a)%*%obs)
	
	## ---------------------------------------------------------------
	## return a dataframe
	## ---------------------------------------------------------------
	res <- data.frame(Estimate=extincttime, lowerCI=lowerCI, upperCI=upperCI)
	
	##return the results
	return(res)
}
OLE(data=dodosightingtimes, alpha=0.05)