phypers {DPQ}

The Four (4) Symmetric 'phyper()' Calls

Description

Compute the four (4) symmetric phyper() calls which mathematically would be identical but in practice typically slightly differ numerically.

Usage

phypers(m, n, k, q = .suppHyper(m, n, k), tol = sqrt(.Machine$double.eps))

Arguments

m

the number of white balls in the urn.

n

the number of black balls in the urn.

k

the number of balls drawn from the urn, hence must be in 0,1,\dots, m+n.

q

vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls. By default all “non-trivial” abscissa values i.e., for which the mathematical value is strictly inside (0,1).

tol

a non-negative number, the tolerance for the all.equal() checks.

Value

a list with components

q

Description of 'comp1'

phyp

a numeric matrix of 4 columns with the 4 different calls to phyper() which are theoretically equivalent because of mathematical symmetry.

Author(s)

Martin Maechler

References

Johnson et al

See Also

R's phyper. In package DPQmpfr, phyperQ() uses (package gmp based) exact rational arithmetic, summing up dhyperQ(), terms computed by chooseZ(), exact (long integer) arithmetic binomial coefficients.

Examples

## The function is defined as
function(m,n,k, q = .suppHyper(m,n,k), tol = sqrt(.Machine$double.eps)) {
    N <- m+n
    pm <- cbind(ph = phyper(q,     m,  n , k), # 1 = orig.
                p2 = phyper(q,     k, N-k, m), # swap m <-> k (keep N = m+n)
                ## "lower.tail = FALSE"  <==>  1 - p..(..)
                Ip2= phyper(m-1-q, N-k, k, m, lower.tail=FALSE),
                Ip1= phyper(k-1-q, n,   m, k, lower.tail=FALSE))

    ## check that all are (approximately) the same :
    stopifnot(all.equal(pm[,1], pm[,2], tolerance=tol),
              all.equal(pm[,2], pm[,3], tolerance=tol),
              all.equal(pm[,3], pm[,4], tolerance=tol))
    list(q = q, phyp = pm)
}


str(phs <- phypers(20, 47, 31))
with(phs, cbind(q, phyp))
with(phs,
     matplot(q, phyp, type = "b"), main = "phypers(20, 47, 31)")

## differences:
with(phs, phyp[,-1] - phyp[,1])
## *relative*
relE <- with(phs, { phM <- rowMeans(phyp); 1 - phyp/phM })
print.table(cbind(q = phs$q, relE / .Machine$double.eps), zero.print = ".")