##
## RNG tests using DKW inequality for rate of convergence
##
## P(sup F_n-F >t)<2exp(-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 nt^2 this agrees with the
## large-sample approximation to the Kolmogorov-Smirnov statistic.
## 


superror<-function(rfoo,pfoo,sample.size,...){
    x<-rfoo(sample.size,...)
    tx<-table(x)
    xi<-as.numeric(names(tx))
    f<-pfoo(xi,...)
    fhat<-cumsum(tx)/sample.size
    max(abs(fhat-f))
}

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=FALSE,stop.on.failure=TRUE){
    rfoo<-eval(as.name(paste("r",stub,sep="")))
    pfoo<-eval(as.name(paste("p",stub,sep="")))
    s<-superror(rfoo,pfoo,sample.size,...)
    if (print.result | print.detail){
        printargs<-substitute(list(...))
        printargs[[1]]<-as.name(stub)
        cat(deparse(printargs))
        if (print.detail)
            cat("\nsupremum error = ",signif(s,2)," with p-value=",min(1,round(pdkwbound(sample.size,s),4)),"\n")
    }
    rval<-(s<qdkwbound(sample.size,pthreshold))
    if (print.result)
        cat(c(" FAILED\n"," PASSED\n",)[rval+1])
    if (stop.on.failure & !rval)
        stop("dkwtest failed")
    rval
  }


dkwtest("binom",size=1,prob=0.2)
dkwtest("binom",size=100,prob=0.2)
dkwtest("binom",size=1,prob=0.2)
dkwtest("binom",size=1,prob=0.8)
dkwtest("binom",size=100,prob=0.8)
dkwtest("binom",size=100,prob=0.2)
dkwtest("binom",size=1,prob=0.2)

dkwtest("pois",lambda=9.5)
dkwtest("pois",lambda=0.95)
dkwtest("pois",lambda=95)
dkwtest("pois",lambda=0.95)

dkwtest("nbinom",size=1,prob=0.2)
dkwtest("nbinom",size=100,prob=0.2)
dkwtest("nbinom",size=1,prob=0.2)
dkwtest("nbinom",size=1,prob=0.8)
dkwtest("nbinom",size=100,prob=0.8)
dkwtest("nbinom",size=100,prob=0.2)
dkwtest("nbinom",size=1,prob=0.2)

dkwtest("norm")
dkwtest("norm",mean=5,sd=3)

dkwtest("gamma",shape=0.1)
dkwtest("gamma",shape=10)
dkwtest("gamma",shape=0.1)
dkwtest("gamma",shape=10)

dkwtest("hyper",m=40,n=30,k=20)
dkwtest("hyper",m=4,n=3,k=2)
dkwtest("hyper",m=40,n=30,k=20)
dkwtest("hyper",m=4,n=3,k=2)
dkwtest("hyper",m=5,n=3,k=2)

dkwtest("signrank",n=1)
dkwtest("signrank",n=10)
dkwtest("signrank",n=1)
dkwtest("signrank",n=30)

dkwtest("wilcox",m=40,n=30)
dkwtest("wilcox",m=4,n=3)
dkwtest("wilcox",m=40,n=30)
dkwtest("wilcox",m=4,n=3)
dkwtest("wilcox",m=5,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=2,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)