## ## RNG tests using DKW inequality for rate of convergence ## ## P(sup | F_n - F | > t) < 2 exp(-2nt^2) ## ## The 2 in front of exp() was derived by Massart. It is the best possible ## constant valid uniformly in t,n,F. For large n*t^2 this agrees with the ## large-sample approximation to the Kolmogorov-Smirnov statistic. ## pdkwbound <- function(n,t) 2*exp(-2*n*t*t) qdkwbound <- function(n,p) sqrt(log(p/2)/(-2*n)) dkwtest <- function(stub = "norm", ..., sample.size = 10000, pthreshold = 0.001, print.result = TRUE, print.detail = TRUE, stop.on.failure = TRUE) { rfoo <- eval(as.name(paste("r", stub, sep=""))) pfoo <- eval(as.name(paste("p", stub, sep=""))) x <- rfoo(sample.size,...) tx <- table(signif(x, 12)) # such that xi will be sort(unique(x)) xi <- as.numeric(names(tx)) f <- pfoo(xi,...) fhat <- cumsum(tx)/sample.size s <- max(abs(fhat-f)) if (print.result || print.detail) { printargs <- substitute(list(...)) printargs[[1]] <- as.name(stub) cat(deparse(printargs)) } rval <- (s < qdkwbound(sample.size,pthreshold)) if (print.result) cat(c(" FAILED\n"," PASSED")[rval+1]) if (print.detail || !rval) cat ("\n", " first =", round(x[4],1), " last =", round(x[sample.size],1), " median =", round(median(x),4), " mean =", round(mean(x),4), " sd =", round(sd(x),4), "\n", " supremum error = ", round(signif(s,2),4), " with p-value =", min(1,round(pdkwbound(sample.size,s),4))) cat("\n\n") if (stop.on.failure && !rval) stop("dkwtest failed") rval } .proctime00 <- proc.time() # start timing RNGkind("default","default") set.seed(123456) dkwtest("binom",size = 1,prob = 0.2) dkwtest("binom",size = 2,prob = 0.2) dkwtest("binom",size = 100,prob = 0.2) dkwtest("binom",size = 1e4,prob = 0.2) dkwtest("binom",size = 1,prob = 0.8) dkwtest("binom",size = 100,prob = 0.8) dkwtest("binom",size = 100,prob = 0.999) dkwtest("pois",lambda = 0.095) dkwtest("pois",lambda = 0.95) dkwtest("pois",lambda = 9.5) dkwtest("pois",lambda = 95) dkwtest("nbinom",size = 1,prob = 0.2) dkwtest("nbinom",size = 2,prob = 0.2) dkwtest("nbinom",size = 100,prob = 0.2) dkwtest("nbinom",size = 1e4,prob = 0.2) dkwtest("nbinom",size = 1,prob = 0.8) dkwtest("nbinom",size = 100,prob = 0.8) dkwtest("nbinom",size = 100,prob = 0.999) dkwtest("norm") dkwtest("norm",mean = 5,sd = 3) dkwtest("gamma",shape = 0.1) dkwtest("gamma",shape = 0.2) dkwtest("gamma",shape = 10) dkwtest("gamma",shape = 20) dkwtest("hyper",m = 40,n = 30,k = 20) dkwtest("hyper",m = 40,n = 3,k = 20) dkwtest("hyper",m = 6,n = 3,k = 2) dkwtest("hyper",m = 5,n = 3,k = 2) dkwtest("hyper",m = 4,n = 3,k = 2) dkwtest("signrank",n = 1) dkwtest("signrank",n = 2) dkwtest("signrank",n = 10) dkwtest("signrank",n = 30) dkwtest("wilcox",m = 40,n = 30) dkwtest("wilcox",m = 40,n = 10) dkwtest("wilcox",m = 6,n = 3) dkwtest("wilcox",m = 5,n = 3) dkwtest("wilcox",m = 4,n = 3) dkwtest("chisq",df = 1) dkwtest("chisq",df = 10) dkwtest("logis") dkwtest("logis",location = 4,scale = 2) dkwtest("t",df = 1) dkwtest("t",df = 10) dkwtest("t",df = 40) dkwtest("beta",shape1 = 1, shape2 = 1) dkwtest("beta",shape1 = 2, shape2 = 1) dkwtest("beta",shape1 = 1, shape2 = 2) dkwtest("beta",shape1 = 2, shape2 = 2) dkwtest("beta",shape1 = .2,shape2 = .2) dkwtest("cauchy") dkwtest("cauchy",location = 4,scale = 2) dkwtest("f",df1 = 1,df2 = 1) dkwtest("f",df1 = 1,df2 = 10) dkwtest("f",df1 = 10,df2 = 10) dkwtest("f",df1 = 30,df2 = 3) dkwtest("weibull",shape = 1) dkwtest("weibull",shape = 4,scale = 4) ## regression test for PR#7314 dkwtest("hyper", m=60, n=100, k=50) dkwtest("hyper", m=6, n=10, k=5) dkwtest("hyper", m=600, n=1000, k=500) ## regression test for non-central t bug dkwtest("t", df=20, ncp=3) ## regression test for non-central F bug dkwtest("f", df1=10, df2=2, ncp=3) cat('Time elapsed: ', proc.time() - .proctime00,'\n')