R/hyperbCvMTest.R

In HyperbolicDist: The Hyperbolic Distribution

### Cramer-von Mises goodness of fit for the hyperbolic distribution
hyperbCvMTest <- function(x, Theta, conf.level = 0.95, ...){
  if (!missing(conf.level) && (length(conf.level) != 1 ||
                               !is.finite(conf.level) ||
                               conf.level < 0 || conf.level > 1))
      stop("conf.level must be a single number between 0 and 1")
  DNAME <- deparse(substitute(x))
  NX <- length(x)
  if (NX < 2)
      stop("not enough x observations")
  METHOD <- "Cramer-von Mises test of hyperbolic distribution"
  zvals <- phyperb(sort(x),Theta)
  STATISTIC <- sum((zvals - ((2*(1:length(x)) - 1)/(2*length(x))))^2) +
           1/(12*length(x))
  xi <- 1/sqrt(1 + Theta[2])
  chi <- Theta[1]*xi/sqrt(1 + Theta[1]^2)
  PARAMETER <- c(xi,chi)
  names(STATISTIC) <- "Wsq"
  names(PARAMETER) <- c("xi","chi")
  names(METHOD) <- "Cramer-von Mises test of hyperbolic distribution"
  pValResult <- hyperbCvMTestPValue(xi, chi, STATISTIC)
  RVAL <- list(statistic = STATISTIC, method = METHOD, data.name = DNAME,
               parameter = PARAMETER, p.value = pValResult$pValue,
               warn = pValResult$warn )
  class(RVAL) <- "hyperbCvMTest"
  print(RVAL, ...)
} ## End hyperbCvMTest()

### Calculate P-Value of Cramer-von Mises test of the hyperbolic distribution
hyperbCvMTestPValue <- function(xi, chi, Wsq, digits = 3){
  xiList <- c(0.99, 0.95, seq( 0.90, 0.1, by = -0.1))
  alphaList <- c(0.25, 0.1, 0.05, 0.025, 0.01)
  chiList <- seq(0, 0.8, by = 0.2)
  exactChi <- FALSE
  exactXi <- FALSE
  warn <- c(FALSE, FALSE)
  data(hyperbWSqTable, envir = environment())
  wsqTable <- get("hyperbWSqTable", envir = environment())
  if(abs(chi) > xi ) stop ("Chi must be less than or equal to Xi")
  tol <- .Machine$double.eps
  if(chi > 0.8){
    chi <- 0.79999
    warn[2] <- TRUE
  }
  for (j in 1:4){
    if(chiList[j] < abs(chi) && (chiList [j + 1] - abs(chi))> tol){
      chiLo <- chiList[j]
      chiUp <- chiList[j + 1]
      wsqTable <- wsqTable[,j] + ((chi - chiList[j])/
                         (chiList[j+1] - chiList[j]))*
                           (wsqTable[,j + 1] - wsqTable[,j])
    }
  }
  for (j in 1:5){
    if(identical(all.equal(chiList[j], chi), TRUE)){
      exactChi <- TRUE
      wsqTable <- wsqTable[,j]
    }
  }
  # correct xi values
  if(xi < 0.9 && abs(chi) > 0.6) xiList[4] <- xiList[4] + 0.01
  if(xi < 0.7 && abs(chi) > 0.4) xiList[6] <- xiList[6] + 0.01
  if(xi < 0.5 && abs(chi) > 0.2) xiList[8] <- xiList[8] + 0.01
  if(xi < 0.3 && abs(chi) > 0.0) xiList[10] <- xiList[10] + 0.01
  for (i in 1:11){
    #print(i)
    if(xi > 0.99){
      xi <- 0.99
      warn[1] <- TRUE
    }
    if(xi < 0.1){
      xi <- 0.1
      warn[1] <- TRUE
    }
    if(identical(all.equal(xiList[i], xi), TRUE)){
      exactXi <- TRUE
      wsqTable <- wsqTable[(5*i - 4):(5*i)]
    }
    if(xi < xiList[i] && xi - xiList[i+1] > tol){
      xiUp <- xiList[i]
      xiLo <- xiList[i+1]
      wsqTable <- wsqTable[(5*i - 4):(5*i + 5)]
    }
  }

  if(exactXi == FALSE) {
    wsqTable <- wsqTable[1:5] +
                ((xi - xiLo)/(xiUp - xiLo))*(wsqTable[6:10] - wsqTable[1:5])
  }

  if(wsqTable[1] > Wsq) pValue <- "> 0.25"
  if(wsqTable[1] == Wsq) pValue <- "0.25"
  if(wsqTable[5] < Wsq) pValue <- "< 0.01"
  if(wsqTable[5] == Wsq) pValue <- "0.01"

  for(i in 1:4){
    if(Wsq >= wsqTable[i] && Wsq < wsqTable[i + 1]){
      pValue <- alphaList[i] + ((Wsq-wsqTable[i])/
                    (wsqTable[i + 1] - wsqTable[i]))*
                      (alphaList[i + 1] - alphaList[i])
    }
  }

  if(is.numeric(pValue)) pValue <- round(pValue, digits)
  pValResult <- list(pValue = pValue, warn = sum(warn))
  pValResult
}## End hyperbCvMTestPValue()

### Print results of Cramer-von Mises goodness-of-fit test
### for the hyperbolic distribution
print.hyperbCvMTest <- function (x, prefix = "\t", ...){
  cat("\n")
  writeLines(strwrap(x$method, prefix = "\t"))
  cat("\n")
  cat("data: ", x$data.name, "\n")
  out <- character()
  if (!is.null(x$statistic))
    out <- c(out, paste(names(x$statistic), "=",
             format(round(x$statistic, 4))))
  if (!is.null(x$parameter))
    out <- c(out, paste(names(x$parameter), "=",
             format(round(x$parameter, 3))))
  if (!is.null(x$p.value)) {
    fp <- as.character(x$p.value)
    out <- c(out, paste("p-value",
             if (substr(fp, 1, 1) == "<" ||substr(fp, 1, 1) == ">"){
               fp
             }else{
               paste("=", fp)
             }))
  }
  if(x$warn != 0){
    warning("Estimated parameters are outside the table.\np-value may be incorrect")
  }
  writeLines(strwrap(paste(out, collapse = ", ")))
  cat("\n")
  invisible(x)
} ## End of print.hyperbCvMTest()