The R Project SVN R

R : Copyright 1999, The R Development Core Team
Version 0.65.0 Under development (unstable) (May 28, 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''
>
> 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 : ......................
>
> ##---  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")
+     }
+ }
>