# evalrt.R
#
# This program implements the LRT for the Gumbel domain of attraction against
# a GPD alternative as implemented in
# Beisel et. al., "Testing the extreme value domain of attraction for 
# distributions of beneficial fitness effects." Genetics (In Press).
#
# Input:  Fitness values relative to a reference genotype (fitnessdata.txt)
# Output: MLE for tau and kappa, p-value based on parametric bootstrap

source("evalrtfunc.R")

# Set number of bootstrap replicates to perform
bootsize=10000

# Read in data set from file
y<-scan(file = "fitnessdata.txt",comment.char = "#")
y<-sort(y)

# Compute sample size n
n<-length(y)

cat("n=",n,"\n")
print(y)

# Perform Likelihood Ratio Test

cat("Calculating MLEsr...\n")

reducedtau<-mean(y)
reduced<-expnegloglike(log(reducedtau))
full<-optim(c(log(reducedtau),0),gpdnegloglike)		
cat("Full Model: ")
cat("tauhat = ",exp(full$par[1]),"\tkappahat = ",exp(full$par[2])-1,"\n")
cat("Reduced Model: ")
cat("tauhat = ",reducedtau,"\tkappahat = ",0,"\n")
cat("l(full)\t\t=\t",-full$value,"\n")
cat("l(reduced)\t\=\t",-reduced,"\n")
lrtts<-(-2*(full$value-reduced))
cat("LRT\t\t=\t",lrtts,"\n\n")

#Generate Parametric Bootstrap of Test Statistic
cat("Generating parametric bootstrap of test statistic...\n")

#Save Data to temp
tempdata<-y

progesswidth<-floor(bootsize/20)
cat("0%")
for(i in 1:((bootsize/progesswidth)-3)){
	cat("-")
}
cat("100%\n")

teststat<-rep(0,bootsize)
for(i in 1:bootsize) {
	if(i%%progesswidth==0) cat("*")
	# Generate data set of size n
	y<-rgpd(n,reducedtau,0)
	y<-sort(y)
	# Perform Likelihood Ratio Test
	boottau<-mean(y)
	reduced<-expnegloglike(log(boottau))
	full<-optim(c(log(boottau),0),gpdnegloglike)
	teststat[i]<- (-2*(full$value-reduced))
}
cat("\n")
y<-tempdata
cat("\nParametric bootstrap complete...\n")
f<-ecdf(teststat)
cat("Testing null hypothesis that data belongs to Gumbel domain (exponential tail)...\n")
pvalue<-1-f(lrtts)
cat("p-value\t\t=\t",pvalue,"\n")