



#################################################################
#I.  Get Data
#################################################################

#loads data with marginals estimated
load("allDatTrainTest.RData")
library(evd)


#choose the marginal transformation used
empirical <- F
d <- 5
fCols <- c("f1", "f2", "f3", "f4", "f5")
uCols <- c("u1", "u2", "u3", "u4", "u5")
if(empirical){
	fCols <- c("fe1", "fe2", "fe3", "fe4", "fe5")
	uCols <- c("emp1", "emp2", "emp3", "emp4", "emp5")
}  
origCols <- c("alx", "arl", "mc", "rt", "ts")
predLoc <- 2
origDataCol <- origCols[predLoc]
frDataCol <- fCols[predLoc]
uDataCol <- uCols[predLoc]
rm(empirical)

#info about the marginal for location 2
scale = 7.56
shape = 0.073
rate =  0.0700467
#thold = 59.44
thold = 57.96577
#marginalPar <- list(scale = scale, shape = shape, rate = rate, thold = thold)
tholdInv <- qgev(1 - rate, 1, 1, 1)


#gets the g, gInv, and gPrime functions needed to transform 
#marginals from orig scale to fr scale and back
overlap <- 4
origData <- dat[, origDataCol]
uData <- dat[, uDataCol]
sl <- sort.list(origData)
sortOrigData <- origData[sl]
sortUData <- uData[sl]
below <- sortOrigData < thold + overlap

#fits log-smoothing spline below thold
gSpline <- smooth.spline(sortOrigData[below], 
				sortUData[below], nknots = 50, spar = 0.3)
xHold <- seq(0, max(sortOrigData[below]), 0.02)
uHold <- predict(gSpline, xHold)$y



#defines forward g and inverse g and forward gPrime functions
gFun <- function(x)
{
	z <- numeric(length(x))
	below <- x <= thold

	#below thold, use spline
	if(length(x[below]) > 0)
	{
		u <- approx(xHold, uHold, x[below])$y
		s <- seq(1, length(x[below]))		
		positive <- s[u > 0 & !is.na(u)]
		negative <- s[u <= 0 & !is.na(u)]		
		z[positive] <- qgev(u[positive],1,1,1)
		z[negative] <- 0
	}

	#above thold, use gpd
	if(length(x[!below]) > 0)
	{
		zAbove <- z[!below]
		u <- 1 - rate*(1 - pgpd(x[!below], thold, scale, shape))
		okayErs <- u < 1
		zAbove[okayErs] <- qgev(u[okayErs], 1, 1, 1)
		zAbove[!okayErs] <- Inf
		z[!below] <- zAbove
	}
	
	return(z)
}


	
gInvFun <- function(z)
{
	x <- numeric(length(z))
	below <- z <= tholdInv

	#below thold, use gSplineInv
	if(length(z[below] > 0))
	{
		u <- pgev(z[below],1,1,1)
		x[below] <- approx(uHold, xHold, u)$y
	}

	#above thold, use gpd inv
	if(length(z[!below] > 0))
	{
		intoFun <- 1 - 1/rate*(1 - pgev(z[!below],1,1,1))
		xAbove <- x[!below]
		xAbove[intoFun < 1] <- qgpd(intoFun[intoFun < 1], thold,
				 scale, shape)
		xAbove[intoFun == 1] <- Inf
		x[!below] <- xAbove
	}
		
	return(x)
}
	


#see notes May 19, 2010	
gDerFun <- function(x)
{
	gPrime <- numeric(length(x))
	below <- x <= thold

	#below thold, use spline
	if(length(x[below]) > 0)
	{
		u <- approx(xHold, uHold, x[below])$y
		s <- seq(1, length(x[below]))		
		positive <- s[u > 0]
		negative <- s[u <= 0]		
		gPrime[negative] <- 0
		gPrime[positive] <- (log(u[positive]))^(-2)*
					1/u[positive]*
					predict(gSpline, x[positive], deriv = 1)$y
	}
		
	#above thold use gpd/Frechet
	if(length(x[!below]) > 0)
	{
		ins <- 1 + shape/scale*(x[!below]-thold)
		gPrime[!below] <- (-log(1 - rate*ins^(-1/shape)))^(-2)*
					  (1-rate*ins^(-1/shape))^(-1)*
					  rate/scale*ins^(-1/shape-1)
	}
								  
	return(gPrime)
}





rm(origData, uData, sortOrigData, sortUData, sl, below, overlap)


##################
#  training set
##################
#ulimately, need angular observations to fit the ang msr model
#if NOT empirical, loads the parametric Frechet margin data
fObsTrain <- train[, fCols]
rowSumTrain <- apply(fObsTrain, 1, sum)
slTrain <- sort.list(rowSumTrain, decreasing = T)
trainSort <- train[slTrain,]
fObsTrainSort <- fObsTrain[slTrain,]
rowSumTrainSort <- rowSumTrain[slTrain]
angTrainSort <- fObsTrainSort/rowSumTrainSort
colnames(angTrainSort) <- c("w1", "w2", "w3", "w4", "w5")
rm(fObsTrain, fObsTrainSort, rowSumTrain, rowSumTrainSort, slTrain)
#keeps only angular measures (sorted) and sorted train set w/ all
#numExc <- 300  #num which are used to fit model
numExc <- 210
#210 if at .93 quantile
#rowSumTrainSort[numExc]
#81.8


##################
#  test set
##################
#ultimately need testSort
fObsTest <- test[, fCols]
#sorts using *only* observed components
fObsTest <- fObsTest[, -predLoc]
rowSumTest <- apply(fObsTest, 1, sum)
slTest <- sort.list(rowSumTest, decreasing = T)
testSort <- test[slTest,]
rm(fObsTest, rowSumTest, slTest)
#numPred <- 150
numPred <- 105
#105 if .93 quantile





#################################################################
#II.  Fit Angular Measure
#################################################################

##################
#  Pairwise Beta
##################
#store parms in pbPar in order alpha, beta(vec)
source("pbDensity.R")
source("pbSpecMsr.R")
source("diriDensity.R")
source("maxLikelihood.R")

init <- rep(2,11)
mlPb <- maxLikelihood(angTrainSort[1:numExc,], model = "pb", 
		init = init)
pbPar <- mlPb$par
pbPar
#numExc = 300
# [1] 0.4314172 7.7498761 0.5321871 0.7036177 0.5824262 1.3867521 1.3397459
# [8] 1.1528388 0.7492636 1.9020939 0.6441370

#numExc = 210
#[1] 0.3723437 6.1124280 0.5125025 0.6426672 0.5596758 0.9580730 #0.9829687 1.0110072 0.7607537 1.6442564 0.5561418

mlPb$value
#numExc = 300
#[1] -803.5344
#numExc = 210
#[1] -592.4468

#a different starting value, yielded same results
#init <- c(0.5, 4, 1, 1, 1, 3, 3, 1, 1, 2, 1)

##look at likelihood contributions
#sl <- sort.list(mlPb$likeContrib)
#keep <- 1:50
#cbind(mlPb$likeContrib[sl[keep]], round(trainSort[sl[keep], c(fCols)], 2))
##low lhood when 1 & 3 high??
#keep <- 250:300
#cbind(mlPb$likeContrib[sl[keep]], round(trainSort[sl[keep], c(fCols)], 2))
##highest lhoods when one location dominates
rm(init, mlPb, sl, keep)




######
# mixed Dirichlet model of Boldi Davison
#results stored in this file--decided against including this model
load("testOptP5K4n210.RData")
mdPar <- testOpt[7, -c(1,2)]
source("mixDiriDensity.R")
source("diriStdDensity.R")


##################
#  Kriging
##################
#estimates (mean, cov) based on all data in train set
#prediction just for the top ones in test set
fitSet <- trainSort[,origCols]
predSet <- testSort[,origCols]

#constants (may need to change for fitSet, predSet)
#number of exceedances for thresholds
#predNExc <- 150
predNExc <- 105
predExc <- predSet[1:predNExc,] 

#estimate mean and covariance from *all* data in fitSet
meanEst <- apply(fitSet, 2, mean)
covEst <- cov(fitSet)

#perform *simple* kriging for prediction exceedances only,
#obtain point estimates
KCov <- covEst[-predLoc, -predLoc]
cCov <- covEst[predLoc, -predLoc]
krigWeights <- solve(KCov) %*% cCov

#gotta worry about mean when making predictions
meanMtx <- matrix(rep(meanEst, dim(predExc)[1]), byrow = T,
			ncol = dim(predExc)[2])
obsResids <- predExc[, -predLoc] - meanMtx[, -predLoc]
krigEsts <- meanEst[predLoc] + t(krigWeights) %*% t(obsResids)

#get prediction variance for conditional distribution
predVar <- covEst[predLoc, predLoc] - t(cCov) %*% solve(KCov) %*% cCov
predVar
#41.42


##################
#  Indicator Kriging
##################
source("marysMonotoneSpline.R")

levelJump <- .25
levelVec <- seq(10,105, levelJump)
numLevels <- length(levelVec)
krigWeightsLevelMtx <- matrix(nrow = length(levelVec), ncol = d-1)

#creates a matrix of kriging weights for each level in levelVec
for(i in seq(1, length(levelVec)))
{
	level <- levelVec[i]
	iMsmts <- fitSet < level
	covLevel <- cov(iMsmts)
	KCovLevel <- matrix(nrow = d, ncol = d)
	KCovLevel[1:(d-1), 1:(d-1)] <- covLevel[-predLoc, -predLoc]
	KCovLevel[d, 1:(d-1)] <- 1
	KCovLevel[1:(d-1), d] <- 1
	KCovLevel[d,d] <- 0
	cCovLevel <- c(covLevel[predLoc, -predLoc], 1)
	temp <- solve(KCovLevel) %*% cCovLevel
	krigWeightsLevelMtx[i,] <- temp[1:(d-1)]
}

#steps thru the predSet getting the indicators for each level
#in levelVec, then multiplying by the krigWeightsMtx and
#thus getting cdf estimate at each level
levelMtx <- matrix(rep(levelVec, d-1), ncol = d-1)
indCdfMtx <- matrix(nrow = length(levelVec), ncol = numPred)

#i <- 2
for(i in seq(1, numPred))
{
	predMtx <- matrix(rep(predSet[i, -predLoc], length(levelVec)), 
					ncol = d-1, byrow = T) 
	iMtx <- 1*(predMtx < levelMtx)
	predLevel <- apply(krigWeightsLevelMtx * iMtx, 1, sum)
	indCdfMtx[, i] <- predLevel
}	



#loops thru creating a matrix of cdf ests and pdf ests
indKrig <- list()
indKrig$xNotMtx <- matrix(rep(levelVec, numPred), ncol = numPred)
indKrig$xNotDensValMtx <- matrix(nrow = numLevels, ncol = numPred)
indKrig$xNotCDF <- matrix(nrow = numLevels, ncol = numPred)
indKrig$xNotConditionalExpVec <- numeric(numPred)

for(i in seq(numPred))
{
	out <- mspl(levelVec, indCdfMtx[,i], kt = 15, knots = F)
	FHat <- out$yhat
	FHat[FHat < 0] <- 0
	FHat[FHat > 1] <- 1
	derHat <- (FHat[-c(1,2)] - FHat[-c(numLevels-1, numLevels)])/
				(2*levelJump)
	derHat <- c(derHat[1], derHat, derHat[numLevels-2])
	indKrig$xNotCDF[,i] <- FHat
	indKrig$xNotDensValMtx[,i] <- derHat
	#trapezoidal rule to get conditional expectation
	integrand <- levelVec * derHat
	trapezoidWghts <- c(1, rep(2, numLevels-2), 1)
	condExp <- sum(integrand*trapezoidWghts)*
			(levelVec[numLevels]-levelVec[1])/(2*numLevels)
	indKrig$xNotConditionalExpVec[i] <- condExp
}
	
	
	
	
	
	
################################################################
#III.  Prediction
################################################################
source("unequalSimpson.R")
source("condDensApprox.R")


############
##  Individual replications
############
##marginal distribution
datPredLoc <- dat[, 1+predLoc]
dens <- density(datPredLoc, bw = 2.5)
x <- dens$x
y1 <- dens$y
y2 <- rate * 1/scale*(1 + shape*(x-thold)/scale)^(-1/shape-1)
y2[y2 > max(y1)] <- max(y1)
blend <- 1/pi*atan(1*(x-thold))+1/2
margY <- (1- blend)*y1 + blend*y2
margX <- x
rm(dens, x, y1, y2)


#m indexes replication
m <- 2
z <- testSort[m, fCols]
pbOut <- condDensApprox(5, z, model = "pairwiseBeta",
			par = pbPar, predLoc = predLoc)
mdOut <- condDensApprox(5, z, model = "mixedDirichlet",
			par = mdPar, predLoc = predLoc)
x <- seq(0, max(pbOut$xNot[is.finite(pbOut$xNot)]), l = 1000)
y <- dnorm(x, krigEsts[m], sqrt(predVar))
yLim <- max(c(pbOut$xDensVal[is.finite(pbOut$xDensVal)],
			  mdOut$xDensVal[is.finite(mdOut$xDensVal)], 
			  y))
plot(pbOut$xNot, pbOut$xDensVal, type = 'l', xlab = "x", 
		ylab = "Conditional Density", ylim = c(0, yLim))
lines(mdOut$xNot, mdOut$xDensVal, col = 3)
lines(x,y, col = 4)
points(testSort[m, (1 + predLoc)], 0, pch = 19, col = 2)
lines(margX, margY, col = 2)
print(pbOut$xConditionalExp)
print(mdOut$xConditionalExp)

m <- 8
z <- testSort[m, fCols]
pbOut <- condDensApprox(5, z, model = "pairwiseBeta",
			par = pbPar, predLoc = predLoc)
mdOut <- condDensApprox(5, z, model = "mixedDirichlet",
			par = mdPar, predLoc = predLoc)
x <- seq(0, max(pbOut$xNot[is.finite(pbOut$xNot)]), l = 1000)
y <- dnorm(x, krigEsts[m], sqrt(predVar))
yLim <- max(c(pbOut$xDensVal[is.finite(pbOut$xDensVal)],
			  mdOut$xDensVal[is.finite(mdOut$xDensVal)], 
			  y))
plot(pbOut$xNot, pbOut$xDensVal, type = 'l', xlab = "x", 
		ylab = "Conditional Density", ylim = c(0, yLim))
lines(mdOut$xNot, mdOut$xDensVal, col = 3)
lines(x,y, col = 4)
points(testSort[m, (1 + predLoc)], 0, pch = 19, col = 2)
lines(margX, margY, col = 2)
print(pbOut$xConditionalExp)
print(mdOut$xConditionalExp)

m <- 40
z <- testSort[m, fCols]
pbOut <- condDensApprox(5, z, model = "pairwiseBeta",
			par = pbPar, predLoc = predLoc)
mdOut <- condDensApprox(5, z, model = "mixedDirichlet",
			par = mdPar, predLoc = predLoc)
x <- seq(0, max(pbOut$xNot[is.finite(pbOut$xNot)]), l = 1000)
y <- dnorm(x, krigEsts[m], sqrt(predVar))
yLim <- max(c(pbOut$xDensVal[is.finite(pbOut$xDensVal)],
			  mdOut$xDensVal[is.finite(mdOut$xDensVal)], 
			  y))
plot(pbOut$xNot, pbOut$xDensVal, type = 'l', xlab = "x", 
		ylab = "Conditional Density", ylim = c(0, yLim))
lines(mdOut$xNot, mdOut$xDensVal, col = 3)
lines(x,y, col = 4)
points(testSort[m, (1 + predLoc)], 0, pch = 19, col = 2)
lines(margX, margY, col = 2)
print(pbOut$xConditionalExp)
print(mdOut$xConditionalExp)


rm(m, z, pbOut, mdOut, x, y)


###########
#  Looping thru top N replications
###########

#loops thru all the exceedances 
#and tabulates results for a qq plot and rank histogram
#done on Frechet scale
source("condDensLoop.R")
pbLoopOut <- condDensLoop(p = 5, 
					zMtx <- testSort[1:numPred,fCols], 
					model = "pairwiseBeta", par = pbPar, 
					predLoc = predLoc)
mdLoopOut <- condDensLoop(p = 5, 
					zMtx <- testSort[1:numPred,fCols], 
					model = "mixedDirichlet", par = mdPar, 
					predLoc = predLoc)


#################################################################
#  RESULTS
#################################################################
actualVal <- testSort[1:numPred, origDataCol]

##rank histogram and qq plot
source("cdfInverse.R")
pbCDI <- cdfInverse(pbLoopOut, T)
mdCDI <- cdfInverse(mdLoopOut, T)

resids <- actualVal - krigEsts
krigCDI <- pnorm(resids, sd = sqrt(predVar))

indKrigCDI <- numeric(numPred)
for(i in seq(1, numPred))
{
	indKrigCDI[i] <- approx(x = indKrig$xNotMtx[,i], 
		y = indKrig$xNotCDF[,i], actualVal[i])$y
}
#band-aid to fix number 12 which comes in way high
indKrigCDI[is.na(indKrigCDI)] <- 1

qtile <- qnorm(seq(0, 1, 0.1), sd = sqrt(predVar))
binOfEach <- cut(resids, qtile, labels = F)
count <- tabulate(binOfEach)

hist(pbCDI, freq = F, breaks = seq(0, 1, .1),
	xlab = "quantile", ylab = "emp. probability",
	main = "", ylim = c(0,2))
#hist(mdCDI, freq = F, breaks = seq(0, 1, .1),
#	xlab = "quantile", ylab = "emp. probability",
#	main = "", ylim = c(0,2))
abline(h = 1, lty = 3)
points(seq(.05, .95, .1), count/numPred/.1, type = 'h', col = 2)

numPts <- dim(pbLoopOut$cumTrapsMtx)[1]
#numReps <- dim(pbLoopOut$cumTrapsMtx)[2]
empirical <- seq(1, numPred)/(numPred+1)
plot(empirical, sort(pbCDI), ylab = "model", 
	xlim = c(0,1), ylim = c(0,1))
#plot(empirical, sort(mdCDI), ylab = "model", 
#	xlim = c(0,1), ylim = c(0,1))
abline(0,1, lty = 3)
points(empirical, sort(krigCDI), col = 2)


#comparing point estimates
#original scale:  actual vs approx condtional expectation
xLim <- range(c(pbLoopOut$xConditionalExpVec, 
			mdLoopOut$xConditionalExpVec, krigEsts))
plot(pbLoopOut$xConditionalExpVec, actualVal, type = 'n', 
			xlim = xLim)
text(pbLoopOut$xConditionalExpVec, actualVal, as.character(seq(1:numPred)))
text(mdLoopOut$xConditionalExpVec, actualVal, as.character(seq(1:numPred)), col = 3)
text(krigEsts, actualVal, col = 4, as.character(seq(1:numPred)))
abline(0,1)
#similar clouds

ptEst <- pbLoopOut$xConditionalExpVec
#ptEst <- mdLoopOut$xConditionalExpVec
#ptEst <- krigEsts
error <- as.vector(round((ptEst - actualVal), 2))
plot(error, type = 'n')
text(seq(1, length(error)), error, as.character(seq(1, length(error))))


############
# RMSE
############
RMSE <- sqrt(mean( (predExc[, predLoc] - krigEsts)^2 ))
RMSE
#[1] 11.00630

RMSE <- sqrt(mean( 
		(predExc[, predLoc] - pbLoopOut$xConditionalExpVec)^2 ))
RMSE
#[1] 12.52352

RMSE <- sqrt(mean( 
		(predExc[, predLoc] - mdLoopOut$xConditionalExpVec)^2 ))
RMSE
#[1] 12.61635

RMSE <- sqrt(mean(
		(predExc[, predLoc] - indKrig$xNotConditionalExpVec)^2 ))
RMSE 
#[1] 11.25895


###########
##  Coverage (high quantile)
###########
source("estimateHighQuantile.R")

qtile <- 0.95

pbQuantile <-estimateHighQuantile(pbLoopOut$zNotMtx, 
	pbLoopOut$cumTrapsMtx,pbLoopOut$zActualMtx[predLoc,], qtile)
1 - mean(pbQuantile)
#0.93

kQuantile <- resids > qnorm(qtile, 0, sd = sqrt(predVar))
1 - mean(kQuantile)
#0.83

mdQuantile <-estimateHighQuantile(mdLoopOut$zNotMtx, 
	mdLoopOut$cumTrapsMtx,mdLoopOut$zActualMtx[predLoc,], qtile)
1 - mean(mdQuantile)
#0.89

#indicator kriging
ikQuantile <- estimateHighQuantile(indKrig$xNotMtx, indKrig$xNotCDF, actualVal, qtile)
1-mean(ikQuantile)
#0.86

#below is actual, pb, krig, md, indKrig
#.90
#.88
#.81
#.84
#.83

#.99
#.97
#.92
#.95
#.90

#.75
#.70
#.65
#.72
#.73

#.50
#.44
#.54
#.52
#.49




###########
##  Likelihood score
###########
likelihoodTestScore <- numeric(numPred)

for(i in seq(1, numPred))
{
	xNot <- pbLoopOut$xNotMtx[,i]
	xNot <- xNot[!is.na(xNot)]
	s <- seq(1, length(xNot))
	lessThan <- xNot < testSort[i, origDataCol]
	maxIndex <- max(s[lessThan])
	likelihoodTestScore[i] <- pbLoopOut$xDensValMtx[maxIndex, i]
}

sum(log(likelihoodTestScore))
# [1] -413.0309

#kriging
krigLhoodTestScore <- dnorm(testSort[1:numPred, origDataCol], mean = krigEsts, sd = sqrt(predVar))
sum(log(krigLhoodTestScore))
# [1] -445.5251


###########
##  CRPS
###########
crpsVec <- matrix(ncol = 3, nrow = numPred)

for(i in seq(1, numPred))
{
	len <- sum(!is.na(pbLoopOut$xNotMtx[,i]))
	F_0 <- as.numeric(pbLoopOut$xNotMtx[1:len,i] > actualVal[i])
	integrand <- (pbLoopOut$cumTraps[1:len, i] - F_0)^2
	#plot(pbLoopOut$xNotMtx[1:len,i], pbLoopOut$cumTraps[1:len,i])
	#points(pbLoopOut$xNot[1:len,i], F_0, col = 2)
#	plot(pbLoopOut$xNotMtx[1:len,i], integrand)
	crpsVec[i,1] <- sum((integrand[1:(len-2)]+integrand[2:(len-1)])/2*diff(pbLoopOut$xNotMtx[1:(len-1),i]))
	cdfKrig <-  pnorm(pbLoopOut$xNotMtx[1:len,i], mean = krigEsts[i], sd = sqrt(predVar))
	integrand2 <- (cdfKrig - F_0)^2
	crpsVec[i,2] <- sum((integrand2[1:(len-2)]+integrand2[2:(len-1)])/2*diff(pbLoopOut$xNotMtx[1:(len-1),i]))
	F_00 <- as.numeric(indKrig$xNotMtx[,i] > actualVal[i])
	integrand3 <- (indKrig$xNotCDF[,i] - F_00)^2
	crpsVec[i,3] <- sum(integrand3*trapezoidWghts)*(levelVec[numLevels]-levelVec[1])/(2*numLevels)
}

apply(crpsVec, 2, sum)
#[1] 716.7346 668.0703 652.1014

#plotting crps
plot(actualVal, crpsVec[,1], ylim = c(0, max(crpsVec)))
points(actualVal, crpsVec[,2], col = 2)
points(actualVal, crpsVec[,3], col = 3)

#checking this surprising result
crpsVec[1:20,]

plotCrps <- function(i)
{
quartz(height = 6, width = 10)
	par(mfrow = c(1,2))
len <- sum(!is.na(pbLoopOut$xNotMtx[,1]))-1
plot(pbLoopOut$xNotMtx[,i], pbLoopOut$cumTrapsMtx[,i], type = 'l')
F_0 <- pbLoopOut$xNotMtx[1:len, i] > actualVal[i]
lines(pbLoopOut$xNotMtx[1:len,i], F_0, col = 2)
cdfKrig <-  pnorm(pbLoopOut$xNotMtx[1:len,i], mean = krigEsts[i], sd = sqrt(predVar))
lines(pbLoopOut$xNotMtx[1:len,i], cdfKrig, col = 3)
lines(indKrig$xNotMtx[,i], indKrig$xNotCDF[,i], col = 4)

m <- i
xMax <- 130
yAng <- pbLoopOut$xDensVal[, m]
xKrig <- seq(0, xMax, l = 1000)
yKrig <- dnorm(x, krigEsts[m], sqrt(predVar))
yIndKrig <- indKrig$xNotDensValMtx[,m]
#yMax <- max(c(yAng, yKrig, yIndKrig), na.rm = T)
yMax <- 0.08

#ang Msr

plot(pbLoopOut$xNot[, m], pbLoopOut$xDensVal[, m], type = 'l',
		xlab = "", 
		ylab = "",
		ylim = c(0, yMax),
		xlim = c(0, 130))
points(testSort[m, (1 + predLoc)], 0, pch = 19, col = 2)

#simple kriging
x <- seq(0, xMax)
y <- dnorm(x, krigEsts[m], sqrt(predVar))
lines(x,y, col = 3)

#indicator kriging
lines(indKrig$xNotMtx[,m], indKrig$xNotDensValMtx[,m], col = 4)


zChar <- as.character(round(testSort[m, origCols]))
#Note:  Arlington listed as location 5
mtext(
	paste("(", zChar[1], ", ", zChar[3], ", ", zChar[4], ", ",
	zChar[5], ", ", zChar[2], ")", sep = ""),
	1, line = 3, cex = 1.2)
}






#########
# QVS
#########
qtile <- .80

pbQVS <- numeric(numPred)
krigQVS <- numeric(numPred)
indKrigQVS <- numeric(numPred)

for(i in seq(1, numPred))
{
	print(i)
	qtPB <- approx(pbLoopOut$cumTrapsMtx[,i], 
				pbLoopOut$xNotMtx[,i], xout = qtile)$y
	qtKrig <- qnorm(qtile, krigEsts[i], sqrt(predVar))
	#for indicator kriging, we have to trim off the ends
	s <- seq(1, length(indKrig$xNotCDF[,i]))
	keep <- indKrig$xNotCDF[,i] != 1 & indKrig$xNotCDF[,i] !=0
	s <- s[keep]
	s <- c(s[1]-1, s, s[length(s)]+1)
	indKrigCDF <- indKrig$xNotCDF[s,i]
	indKrigX <- indKrig$xNotMtx[s,i]
	s2 <- seq(1, length(indKrigCDF))
	s2 <- s2[!duplicated(indKrigCDF)]
	indKrigCDF <- indKrigCDF[s2]
	indKrigX <- indKrigX[s2]
	qtIndKrig <- approx(indKrigCDF, 
				indKrigX, xout = qtile)$y

	diffPB <- actualVal[i] - qtPB
	diffKrig <- actualVal[i] - qtKrig
	diffIndKrig <- actualVal[i] - qtIndKrig
	
	if(diffPB < 0)
	{ pbQVS[i] <- (qtile - 1)*diffPB }
	else
	{ pbQVS[i] <- qtile*diffPB }
	
	if(diffKrig < 0)
	{ krigQVS[i] <- (qtile - 1)*diffKrig }
	else
	{ krigQVS[i] <- qtile*diffKrig }
	
	if(diffIndKrig < 0)
	{ indKrigQVS[i] <- (qtile - 1)*diffIndKrig }
	else
	{ indKrigQVS[i] <- qtile*diffIndKrig }
}
	
qtile
sum(pbQVS)
sum(krigQVS)
sum(indKrigQVS)
				
#.99
###> sum(pbQVS)
###[1] 40.97406
###> sum(krigQVS)
###[1] 65.8088
###> sum(indKrigQVS)
###[1] 67.80456

#.95
###> sum(pbQVS)
###[1] 134.7734
###> sum(krigQVS)
###[1] 170.0366
###> sum(indKrigQVS)
###[1] 153.4065

#.90
###> sum(pbQVS)
###[1] 225.6814
###> sum(krigQVS)
###[1] 246.2555
###> sum(indKrigQVS)
###[1] 238.6319

#.85
###> sum(pbQVS)
###[1] 295.8181
###> sum(krigQVS)
###[1] 300.2415
###> sum(indKrigQVS)
###[1] 303.9192

#[1] 0.75
# # # > sum(pbQVS)
# # # [1] 398.9665
# # # > sum(krigQVS)
# # # [1] 378.2651
# # # > sum(indKrigQVS)
# # # [1] 377.1899

###> qtile
###[1] 0.5
###> sum(pbQVS)
###[1] 502.5117
###> sum(krigQVS)
###[1] 444.8451
###> sum(indKrigQVS)
###[1] 452.6777

###> qtile
###[1] 0.1
###> sum(pbQVS)
###[1] 235.0169
###> sum(krigQVS)
###[1] 237.3748
###> sum(indKrigQVS)
###[1] 212.5893

###> qtile
###[1] 0.05
###> sum(pbQVS)
###[1] 159.7788
###> sum(krigQVS)
###[1] 160.9128
###> sum(indKrigQVS)
###[1] 137.3183



