R/jointExceedanceCurves.R

In harrysouthworth/texmex: Statistical Modelling of Extreme Values

#' Joint exceedance curves
#'
#' Calculate bivariate joint exceedance curves
#'
#' Calculates pairs of points (x,y) for which the point exceedance
#' probability P(X>x and Y>y) is constant.  This is available only in
#' two dimensions: for higher dimensional data, the bivariate margin
#' will be used and other variables ignored. Takes as input either a
#' two column matrix of observations, output from \code{mexMonteCarlo}
#' (in which case samples from all fitted models are used to calculate curves)
#' or output from a call to the \code{predict} method for an object of
#' class \code{mex} (in which case just the single fitted model is used
#' for estimation, with the importance sample generated in the call
#' to \code{predict} being used to calculate the joint exceedance curve).
#'
#' @return Returns an object of class \code{jointExcCurve}. This is a
#' list of length two, one for each variable for which the curve is
#' calculated.  Each item of the list is a vector of coordinate values
#' for the variable in question.  Attributes include \code{names} and
#' the exceedance probability used to calculate the curve \code{ExceedanceProb}.
#'
#' The curve is calculated by finding pairs of points (x,y) for which the empirical
#' probability P(X>x, Y>y) of both variables exceeding their corresponding value is equal
#' to the specified \code{ExceedanceProb}.  Note that when this is calculated
#' for an object of class \code{predict.mex} (returned by a call to the \code{predict}
#' method for an object of class \code{mex}) then the exceedance probability is
#' interpreted as the UNCONDITIONAL exceedance probability of the importance sample,
#' ie the probability of sampled values occurring from the original modelled joint
#' distribution, and NOT the conditional distribution used to generate the
#' importance sample.
#'
#' Estimated curve can be added to a ggplot of the data (and/or importance sample) by using the function \code{geom_jointExcCurve}, see examples below.
#'
#' @examples
#' \donttest{
#'  # for data frame of raw data
#'  Sigma <- matrix(c(1, .5, .5, 1), ncol=2)
#'  m1 <- rmvnorm(5000,sigma=Sigma)
#'  m1 <- as.data.frame(m1)
#'  j1 <- JointExceedanceCurve(m1,0.01)
#'  j2 <- JointExceedanceCurve(m1,0.005)
#'  j3 <- JointExceedanceCurve(m1,0.001)
#'  ggplot(m1,aes(V1,V2)) + geom_point(colour="dark blue",alpha=0.5) +
#'  geom_jointExcCurve(j1,colour="orange") +
#'  geom_jointExcCurve(j2,colour="orange") +
#'  geom_jointExcCurve(j3,colour="orange")
#'
#'  # using importance sample generated by call to predict for object of class mex
#'  m <- mex(winter,mqu=0.7,dqu=0.7,which="NO")
#'  m2 <- predict(m,nsim=5000,pqu=0.999)
#'  g <- ggplot(m2,plot.=FALSE)
#'  j4 <- JointExceedanceCurve(m2,0.0005,which=c("NO","NO2"))
#'  j5 <- JointExceedanceCurve(m2,0.0002,which=c("NO","NO2"))
#'  j6 <- JointExceedanceCurve(m2,0.0001,which=c("NO","NO2"))
#'  g[[2]] +
#'      geom_jointExcCurve(j4,aes(NO,NO2),col="orange") +
#'      geom_jointExcCurve(j5,aes(NO,NO2),col="orange") +
#'      geom_jointExcCurve(j6,aes(NO,NO2),col="orange")
#'
#' # for augmented dataset, generated by MC sampling from collection of fitted H+T models
#'  m <- mexAll(winter,mqu=0.7,dqu=rep(0.7,5))
#'  m3 <- mexMonteCarlo(nSample=5000,mexList=m)
#'  j7 <- JointExceedanceCurve(m3,0.05,which=c("NO","NO2"))
#'  j8 <- JointExceedanceCurve(m3,0.02,which=c("NO","NO2"))
#'  j9 <- JointExceedanceCurve(m3,0.01,which=c("NO","NO2"))
#'  ggplot(as.data.frame(m3$MCsample[,c("NO","NO2")]),aes(NO,NO2)) +
#'      geom_point(col="light blue",alpha=0.5) +
#'      geom_jointExcCurve(j7,col="orange") +
#'      geom_jointExcCurve(j8,col="orange") +
#'      geom_jointExcCurve(j9,col="orange")
#' }
#' @aliases geom_jointExcCurve JointExceedanceCurve
#'
#' @usage JointExceedanceCurve(Sample, ExceedanceProb,...)
#' \method{print}{jointExcCurve}(x,...)
#'
#' @export JointExceedanceCurve
#'
JointExceedanceCurve <- function(Sample, ExceedanceProb,...) {
    theCall <- match.call()
    UseMethod("JointExceedanceCurve", Sample)
}

#' @rdname JointExceedanceCurve
#' @export
JointExceedanceCurve.default <- function(Sample, ExceedanceProb,n=50,x=NULL,...) {
    s <- calcJointExceedanceCurve(Sample,ExceedanceProb,n,x)
    names(s) <- names(Sample)
    attr(s,"ExceedanceProb") <- ExceedanceProb
    s
}

#' @rdname JointExceedanceCurve
#' @export
#' @method JointExceedanceCurve mexMC
#' @param which Vector length two identifying which margins to use for joint exceedance curve estimation. Can be integer vector, giving column numbers of original data matrix, or character vector identifying variables by name (these must match column names in original data).
#' @param ... Further aguments to be passed to methods
JointExceedanceCurve.mexMC <- function(Sample, ExceedanceProb,n=50,x=NULL,which=1:2,...) {
    S <- Sample$MCsample[,which]
    Sample <- as.matrix(S)
    Sample <- Sample[!is.na(Sample[,1]),]
    Sample <- Sample[!is.na(Sample[,2]),]

    s <- calcJointExceedanceCurve(Sample,ExceedanceProb,n,x)
    names(s) <- names(S)
    attr(s,"ExceedanceProb") <- ExceedanceProb
    s
}

#' @rdname JointExceedanceCurve
#' @export
#' @method JointExceedanceCurve predict.mex
JointExceedanceCurve.predict.mex <- function(Sample, ExceedanceProb,n=50,x=NULL,which=1:2,...) {
    CondExceedanceProb <- 1-Sample$pqu
    if(is.numeric(which)){ # since column order of predict order sample is not original data column order
        d <- dim(Sample$data$simulated)[2]
        w <- Sample$which
        which <- order(c(w,c(1:d)[-w]))[which]
    }
    S <- Sample$data$simulated[,which]
    Sample <- as.matrix(S)
    if(ExceedanceProb > CondExceedanceProb) stop("ExceedanceProb must be less than the probability of exceeding the threshold used for importance sampling in the call to predict")

    s <- calcJointExceedanceCurve(Sample,ExceedanceProb/CondExceedanceProb,n,x)
    names(s) <- names(S)
    attr(s,"ExceedanceProb") <- ExceedanceProb
    s
}

#' @rdname JointExceedanceCurve
#' @export calcJointExceedanceCurve
#' @param Sample Monte Carlo (or other) sample from which to calculate joint exceedance curve
#' @param ExceedanceProb Takes values between 0 and 1, constant value of joint exceedance probability for which the curve will be calculated
#' @param n If \code{x=NULL} then this is HALF the number of points at which the curve will be estimated (ie the curve is calculated at 2n locations)
#' @param x If specified by the user, the values of in the first dimension of \code{Sample} at which to calculate the curve. Defaults to \code{NULL} otherwise should be a numeric vector within the range of the first dimension of \code{Sample}.
calcJointExceedanceCurve  <- function(Sample,ExceedanceProb,n=50,x=NULL) {
    # mx, my are marginal upper limits
    # px, py are plotting points
    mx <- quantile(Sample[,1],1-ExceedanceProb)
    my <- quantile(Sample[,2],1-ExceedanceProb)
    if(is.null(x)){
        px <- seq(min(Sample[,1]),mx,length=n)
        py <- seq(min(Sample[,2]),my,length=n)
    } else {
        px <- x
    }

    f <- function(z,x,y) {
        sapply(z,function(z){
            g <- function(w) mean(x > z & y > w) - ExceedanceProb
            if(g(range(y)[1]) * g(range(y)[2]) < 0) {
                out <- uniroot(g,lower=range(y)[1],upper=range(y)[2])$root
            } else {
                out <- min(y[x>z])
            }
            out
        }
        )
    }

    #calculate curve values at plotting points
    cx <- f(px,Sample[,1],Sample[,2])
    if(is.null(x)){
        cy <- f(py,Sample[,2],Sample[,1])
        res <- data.frame(x=c(px,cy),y=c(cx,py))
        # sort results
        res <- res[order(res[,1]),]
        res[,2] <- sort(res[,2],decreasing=TRUE) # gets rid of errors introduced by tied values which are sorted on in previous line.  Shape has to be convex so must be decreasing.
        row.names(res) <- NULL
    } else {
        res <- data.frame(x=px,y=cx)
    }
    oldClass(res) <- "jointExcCurve"
    res
}

#' @rdname JointExceedanceCurve
#' @export
print.jointExcCurve <- function(x, ...){
    P <- attributes(x)$ExceedanceProb
    cat("\n Estimated curve with constant joint exceedance probability equal to",P,"\n")
    res <- as.data.frame(cbind(x[[1]],x[[2]]))
    colnames(res) <- attributes(x)$names
    print(res)
    invisible(x)
}

#' @rdname JointExceedanceCurve
#' @export geom_jointExcCurve
geom_jointExcCurve <- function(x,...){
    dat <- as.data.frame(cbind(x[[1]],x[[2]]))
    colnames(dat) <- attributes(x)$names
    geom_line(data=dat,...)
}