Nothing
### 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()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.