The R Project SVN R

####    d|ensity
####    p|robability (cumulative)
####    q|uantile
####    r|andom number generation
####
####    Functions for  ``d/p/q/r''

if(!interactive()) .Random.seed <- c(0,rep(7654, 3))

###--- Discrete Distributions: Simple Consistency Checks  pZZ = cumsum(dZZ)

## Currently, just Wilcoxon  [should do this for all !]
is.sym <- TRUE
for(n in rpois(5, lam=6))
    for(m in rpois(15, lam=8))
    {
        x <- -1:(n*m + 1)
        fx <- dwilcox(x, n, m)
        Fx <- pwilcox(x, n, m)
        is.sym <- is.sym & all(fx == dwilcox(x, m, n))
        eq <- all.equal(Fx, cumsum(fx), tol= 1e-14)
        if(!is.logical(eq) || !eq) print(eq)
    }
is.sym

##--- Cumulative Poisson '==' Cumulative Chi^2 :
##--- Abramowitz & Stegun, p.941 :  26.4.21 (26.4.2)
n1 <- 20; n2 <- 16
for(lambda in rexp(n1))
    for(k in rpois(n2, lambda)) {
    tst <- all.equal(1 - pchisq(2*lambda, 2*(k+1)),
             sum(dpois(0:k, lambda=lambda)))
    if(!(is.logical(tst) && tst))
        cat("lambda=", format(lambda),".  k =",k, " --> tst=", tst,"\n")
    }

##--- Cumulative Binomial '==' Cumulative F :
##--- Abramowitz & Stegun, p.945-6;  26.5.24  AND  26.5.28 :
n0 <- 50; n1 <- 16; n2 <- 20; n3 <- 8
for(n in rbinom(n1, size = 2*n0, p = .4)) {
    cat("n=",n,": ")
    for(p in c(0,1,rbeta(n2, 2,4))) {
    cat(".")
    for(k in rbinom(n3, size = n,  prob = runif(1))) {
        ## For X ~ Bin(n,p), compute 1 - P[X > k] = P[X <= k]  in two ways:
        tst <- all.equal(if(k==n || p==0) 1 else
                 pf((k+1)/(n-k)*(1-p)/p, df1=2*(n-k), df2=2*(k+1)),
                 sum(dbinom(0:k, size = n, prob = p)))
        if(!(is.logical(tst) && tst))
        cat("n=", n,"; p =",format(p),".  k =",k, " --> tst=",tst,"\n")
    }
    }
    cat("\n")
}

##---  Gamma (incl. chi^2) Density :
x <- round(rgamma(100, shape = 2),2)
for(sh in round(rlnorm(30),2)) {
    Ga <- gamma(sh)
    for(sig in round(rlnorm(30),2)) {
        tst <- all.equal((d1 <- dgamma(  x,   shape = sh, scale = sig)),
                         (d2 <- dgamma(x/sig, shape = sh, scale = 1) / sig),
                         tol = 1e-15)
        if(!(is.logical(tst) && tst))
            cat("ERROR: dgamma() doesn't scale:",tst,"\n",
                "  x =", formatC(x),"\n  shape,scale=",formatC(c(sh, sig)),"\n")
        tst <- all.equal(d1, (d3 <- 1/(Ga * sig^sh) * x^(sh-1) * exp(-x/sig)),
                         tol= 1e-14)
        if(!(is.logical(tst) && tst))
            cat("NOT Equal:",tst,"\n x =", formatC(x),
                "\n  shape,scale=",formatC(c(sh, sig)),"\n")
    }
}