R : Copyright 1999, The R Development Core Team
Version 0.64.0 Patched (unreleased snapshot) (May 7, 1999)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type	"?license" or "?licence" for distribution details.

R is a collaborative project with many contributors.
Type	"?contributors" for a list.

Type	"demo()" for some demos, "help()" for on-line help, or
    	"help.start()" for a HTML browser interface to help.
Type	"q()" to quit R.

> ####	d|ensity
> ####	p|robability (cumulative)
> ####	q|uantile
> ####	r|andom number generation
> ####
> ####	Functions for  ``d/p/q/r''
> 
> source(paste(getenv("SRCDIR"),"all.equal.R", sep="/"))
> 
> 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        
[1] TRUE
> 
> ##--- 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")
+ }
n= 33 : ......................
n= 32 : ......................
n= 33 : ......................
n= 44 : ......................
n= 39 : ......................
n= 44 : ......................
n= 41 : ......................
n= 46 : ......................
n= 30 : ......................
n= 43 : ......................
n= 44 : ......................
n= 33 : ......................
n= 37 : ......................
n= 39 : ......................
n= 37 : ......................
n= 50 : ......................
>