library(evd)

#simulates 1 mv logistic rv, three-dimensional
#alpha is dependence parameter of mv logistic, 
#dep increases as alpha goes to zero
alpha <- .3
n <- 1
z <- rmvevd(n = n, dep = alpha, model = "log", d = 3, 
			mar = c(1,1,1))
print(z)

source("h.R")
source("compareCondDens.R")
source("numericalIntegration.R") #note function is called "unequalSimpson"

out <- compareCondDens(z, dep = alpha, plot = T)


#small one, bad approx
set.seed(38)
zSmall <- rmvevd(n = n, dep = alpha, model = "log", d = 3, 
			mar = c(1,1,1))
print(zSmall)
#[1] 0.2250443 0.2428516 0.3146199
outSmall <- compareCondDens(zSmall, dep = alpha, plot = T)


#moderate, pretty good approx
set.seed(233)
zMod <- rmvevd(n = n, dep = alpha, model = "log", d = 3, 
			mar = c(1,1,1))
print(zMod)
outMod <- compareCondDens(zMod, dep = alpha, plot = T)


#big one, very good approx
set.seed(269)
zBig <- rmvevd(n = n, dep = alpha, model = "log", d = 3, 
			mar = c(1,1,1))
print(zBig)
outBig <- compareCondDens(zBig, dep = alpha, plot = T)








#repeated simulation
set.seed(38)
numSims <- 5000
loopOut <- matrix(nrow = numSims, ncol = 5)
colnames(loopOut) <- c("z1", "z2", "z3", "z1 + z2", "p")

for(i in seq(1, numSims))
{
	print(i)
	z <- rmvevd(n = n, dep = alpha, model = "log", d = 3, 
			mar = c(1,1,1))
	out <- compareCondDens(z, dep = alpha, plot = F)
	s <- seq(1, dim(out)[1])
	greater <- out[,1] > z[3]
	line <- min(s[greater])
	x1 <- out[(line-1), 1]
	x2 <- out[line, 1]
	y1 <- out[(line-1), 3]
	y2 <- out[line, 3]
	loopOut[i, 1:3] <- z
	loopOut[i, 4] <- z[1] + z[2]
	loopOut[i, "p"] <- y1 + (y2 - y2)/(x2 - x1)*(z[3] - x1)
}	

sl <- sort.list(loopOut[, 4], decreasing = T)
sortLoopOut <- loopOut[sl, ]

hist(sortLoopOut[1:1000, 5], freq = F, main = "", xlab = "")
abline(h = 1, lty = 2)
varFreq <- 1000*.9*.1
varProp <- varFreq/100^2
sdProp <- sqrt(varProp)
abline(h = 1 + 1.96*sdProp, lty = 3)
abline(h = 1 - 1.96*sdProp, lty = 3)
