.Rhistory

#### 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()