tests/hypergeom.R

In DPQmpfr: DPQ (Density, Probability, Quantile) Distribution Computations using MPFR

library(DPQ)
library(DPQmpfr)

## Imported from 'DPQ', support of the hypergeometric distrib. as function of its parameters:
## .suppHyper <- DPQ::.suppHyper
.suppHyper <- function(m,n,k) max(0, k-n) : min(k, m)

## case where phyper(*, log=TRUE) fails :
phyper(11, 15, 0, 12) # 0 -- correct
phyper(11, 15, 0, 12, log=TRUE) # NaN  || correct value:  log(0) = -Inf
(phyp5.0.12 <- cumsum(dhyperQ(0:12, m=15,n=0,k=12)))

require(Rmpfr)
.N <- asNumeric # short form for here

m <- 1000; n <- 1001; k <- 1500
dhyp.1k <- dhyperQ(.suppHyper(m,n,k), m,n,k)
phyp.1k <- cumsum(dhyp.1k)
## Rel.Error: number of correct digits
.N(-log10(abs(1 - phyper(.suppHyper(m,n,k), m,n,k) / .bigq2mpfr(phyp.1k, 99))))

## even a bit larger
m <- 9999; n <- 500; k <- 10100
range(x <- .suppHyper(m,n,k)) # 9600 9999
dQ <- dhyperQ(x, m,n,k) # these add exactly to 1 :
stopifnot(sum(dQ) == 1)


dM <- .bigq2mpfr(dQ, precB=128)
tail( log(dM) ) # -1395.7599 ... -1443.7188576

1 - phyper(x, m,n,k, log.p=TRUE) / log(cumsum(dM)) # from 0 --> 1  ???
## -->
## this is *bad*: cumsum(dM) suffers from cancellation still!
cbind(x, phyper = phyper(x, m,n,k), cum.dQ = asNumeric(cumsum(dM)),
      phyp.log  = phyper(x, m,n,k, log.p=TRUE), log.cumQ = asNumeric(log(cumsum(dM))))

## this shows the cancellation:
cbind(x,
      cumD = .N(1 - cumsum(dM)),
      revCS= .N(rev(cumsum(rev(dM)))))[ 1:99, ]
## Using bigrational ("Q") is still exact:   ---- well there's problem somewhere !
cbind(x,
      cumD = .N(1 - cumsum(    dQ)),
      revCS= .N(rev(cumsum(rev(dQ)))))[ 1:99, ]

## seems good
I_CSd <- .bigq2mpfr(1 - cumsum(dQ), precB=1024)
roundMpfr(I_CSd, precBits = 8)
.N(log(I_CSd)) # should be "exact"

## rel.err :
relE.up.ln <- .N(1- phyper(x, m,n,k, lower.tail = FALSE, log.p = TRUE) / log(I_CSd))
plot(x, relE.up.ln, type = "o", cex = 1/2)
abline(h = c(-1,1)*.Machine$double.eps, col = "thistle", lty=2)
plot(x, log(I_CSd), type="l")
## rel.err of phyper(*, log=TRUE) : --> is exploding!
1- phyper(x, m,n,k, log.p=TRUE) / log(cumsum(dM))