dir <- "/home/orzechma/Documents/SeniorThesis/discLog/"
library(numbers)
library(nortest)
#require(graphics)

# Read Data
dat <- read.csv(file=paste(dir, "discLogData.csv", sep=""), header=TRUE, stringsAsFactors=FALSE)
attach(dat)

#n <- dat$n

#true <- ((1 + n) * (sqrt(pi)*gamma(n/2)*n - 2*gamma(1/2 + n/2))) / (2*gamma(1/2 + n/2)*n)
trueTL <- n*sqrt(pi)*exp(lgamma(n/2) - lgamma(n/2 + 1/2))/4 + sqrt(pi)*exp(lgamma(n/2) - lgamma(n/2 + 1/2))/2 - 1 - 1/n
trueCL <- n*sqrt(pi)*exp(lgamma(n/2) - lgamma(n/2 + 1/2))/4
trueRL <- n*sqrt(pi)*exp(lgamma(n/2) - lgamma(n/2 + 1/2))/2 - 1/n + sqrt(pi)*exp(lgamma(n/2) - lgamma(n/2 + 1/2))/2 - 1

trials <- vector()
for (i in 1:length(n)) {
	trials[i] <- eulersPhi((n[i] - 1)/2)
}
#trials <- eulersPhi((n - 1)/2)

#result <- t.test(dat$aveRL, true, mu=0, paired=TRUE)

tTL <- (sqrt(trials)*(aveTL - trueTL))/sqrt(varTL)
tCL <- (sqrt(trials)*(aveCL - trueCL))/sqrt(varCL)
tRL <- (sqrt(trials)*(aveRL - trueRL))/sqrt(varRL)
#any(is.na(true))
#warnings()
pTL <- 2*pt(-abs(tTL), df=trials-1)
pCL <- 2*pt(-abs(tCL), df=trials-1)
pRL <- 2*pt(-abs(tRL), df=trials-1)

rejectRL <- vector()
count <- 0
for(i in 1:length(pRL)){
	if (pRL[i] > 0.05) {
		count <- count + 1
		rejectRL[i] <- "NO"
	} else {
		rejectRL[i] <- "YES"
	}
}

ratio <- count / length(n)
#ratio

mean(tRL)
sd(tRL)
shapiro.test(tRL)
ad.test(tRL)

shapiro.test(tTL)
ad.test(tTL)
shapiro.test(tCL)
ad.test(tCL)

qqnorm(tTL, main="Normal Q-Q Plot Tail Length")
qqline(tRL, col="red")
plot(density(tTL))

qqnorm(tCL, main="Normal Q-Q Plot Cycle Length")
qqline(tCL, col="red")
plot(density(tCL))

qqnorm(tRL, main="Normal Q-Q Plot Rho Length")
qqline(tRL, col="red")
plot(density(tRL))

df = data.frame(n, tRL, pRL, rejectRL)

write.table(df, file=paste(dir, "discLogTestValues.csv", sep=""), append=FALSE, quote=FALSE, sep="\t\t", eol="\n", row.names = FALSE)
#warnings()