# evalrtmc.R
#
# This program implements the LRT for the Gumbel domain of attraction against
# a GPD alternative incooperating measurement error 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=1000

# Set Monte Carlo Replicates
reps=1000

# Set normal error standard deviation  
sigma<-.1

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

# Estimate tau and kappa using model without error
cat("Calculating MLEs for data under GPD without error...\n")
m<-optim(c(log(1/mean(y)),0),gpdnegloglike)
tauhat<-m$par[1]
kappahat<-m$par[2]
cat("tauhat = ",exp(tauhat),"\tkappahat = ",exp(kappahat)-1,"\n")

# Generate Simulated Data for Monte Carlo Integration
cat("Generating simulated data for Monte Carlo Approximation...\n")
x<-rgpd(n*reps,exp(tauhat),exp(kappahat)-1)
dim(x)<-c(n,reps)
x<-apply(x,2,FUN=sort)
cat(reps," data sets generated with sample size ",n,"\n")

# Perform Likelihood Ratio Test

cat("Calculating MLEs with error...\n")

full<-optim(c(tauhat,kappahat),gpdneglogmclike)		
reduced<-optimize(expneglogmclike,c(log(1/(10*mean(y))),log(10/mean(y))))
cat("Full Model: ")
cat("tauhat = ",exp(full$par[1]),"\tkappahat = ",exp(full$par[2])-1,"\n")
cat("Reduced Model: ")
reducedtau<-reduced$minimum
cat("tauhat = ",exp(reduced$minimum),"\tkappahat = ",0,"\n")
cat("l(full)\t\t=\t",-full$value,"\n")
cat("l(reduced)\t\=\t",-reduced$objective,"\n")
lrtts<-(-2*(full$value-reduced$objective))
cat("LRT\t\t\t=\t",lrtts,"\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 with normal error	
	y<-rgpd(n,exp(reducedtau),0)
	newy<- -y
	while(length(newy[newy<0])>0) {
		newy<- -y
		newy<- -newy+rnorm(n,0,sigma)
	}	
	y<-newy
	y<-sort(y)

	# Estimate tau and kappa using model without error
	m<-optim(c(log(1/mean(y)),0),gpdnegloglike)
	tauhat<-m$par[1]
	kappahat<-m$par[2]

	# Generate Simulated Data for Monte Carlo Integration
	x<-rgpd(n*reps,exp(tauhat),exp(kappahat)-1)
	dim(x)<-c(n,reps)
	x<-apply(x,2,FUN=sort)
	full<-optim(c(tauhat,kappahat),gpdneglogmclike)		
        reduced<-optimize(expneglogmclike,c(log(1/(10*mean(y))),log(10/mean(y))))
	teststat[i]<- (-2*(full$value-reduced$objective))
}
cat("\n")
y<-tempdata
cat("\nParametric bootstrap complete...\n\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")