R Under development (unstable) (2023-11-02 r85465 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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 more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #### Used to be part of /u/maechler/R/MM/NUMERICS/hyper-dist.R -- till Jan.24 2020 > > library(DPQ) > > #### First version committed is the code "as in Apr 1999": > #### file | 11782 | Apr 22 1999 | hyper-dist.R > > options(rErr.eps = 1e-30) > rErr <- function(approx, true, eps = getOption("rErr.eps", 1e-30)) + { + ifelse(Mod(true) >= eps, + 1 - approx / true, # relative error + true - approx) # absolute error (e.g. when true=0) + } > > if(!dev.interactive(orNone=TRUE)) pdf("hyper-dist-ex.pdf") > > > ### ----------- ----------- > > rErr.phypDiv <- function(q, m,n,k) { + ph <- phyper(q, m=m, n=n, k=k) + cbind( + rERR.Ibeta= rErr(phyperIbeta (q, m=m, n=n, k=k) , ph), + rERR.as151= rErr(phyperApprAS152 (q, m=m, n=n, k=k) , ph), + rERR.1mol = rErr(phyper1molenaar (q, m=m, n=n, k=k) , ph), + rERR.2mol = rErr(phyper2molenaar (q, m=m, n=n, k=k) , ph), + rERR.Peizer=rErr(phyperPeizer (q, m=m, n=n, k=k) , ph)) + } > > ## Plotting of the rel.error: > p.rErr.phypDiv <- function(m,n,k, q = .suppHyper(m,n,k), abslog = FALSE, logx= "", type = "b") + { + rE <- rErr.phypDiv(m=m, n=n, k=k, q=q) + nc <- ncol(rE) + cc <- adjustcolor(1:nc, 0.5) + ppch <- as.character(1:nc) + cl <- sys.call() + cl[[1]] <- quote(phyper) + matplot(q, if(abslog) abs(rE) else rE, type = type, lwd=2, col=cc, pch=ppch, + ylab = if(abslog) quote(abs(rE)) else quote(rE), + log = paste0(logx, if(abslog) "y" else ""), + main = sprintf("%s", deparse1(cl)))# + mtext(paste(if(abslog) "|rel.Err|" else "rel.Err.", + "of diverse phyper() pbinom* approximations")) + if(!abslog) abline(h=0, col=adjustcolor("gray", 0.5), lty=2) + legend("topright", legend = paste(1:ncol(rE), colnames(rE), sep=": "), + lwd=2, col=cc, lty=1:4, pch=ppch, bty="n") + invisible(rE) + } > > > k <- c(10*(1:9),100*(1:9),1000*(1:9)) > k1 <- k > options(digits=4) > > rErr12 <- p.rErr.phypDiv(q=k, 2*k,2*k,2*k) > p.rErr.phypDiv(q=k, 2*k,2*k,2*k, abslog=TRUE) > > k <- c(10*(7:9),100*(1:9),1000*(1:9), 1e4*(1:9)) > k2 <- k > > p.rErr.phypDiv(q=k, 2*k,2*k,2*k, abslog=TRUE, logx="x") > > k <- c(10*(1:9),100*(1:3)) # small k only (afterwards --> all phyper*() == 1 !)(revert) > ## Here, the Ibeta fails (NaN) ; *also* Molenaar's (???? contradiction to below) > (rErr1215 <- p.rErr.phypDiv(q=k, 1.2*k, 2*k, 1.5*k)) rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peizer [1,] NaN 1.287e-02 3.718e-05 6.630e-05 -2.392e-06 [2,] NaN 1.805e-03 3.092e-07 4.254e-07 -3.749e-09 [3,] NaN 2.646e-04 1.994e-09 2.535e-09 -9.461e-12 [4,] NaN 4.002e-05 1.192e-11 1.460e-11 -2.864e-14 [5,] NaN 6.175e-06 6.894e-14 8.271e-14 0.000e+00 [6,] NaN 9.664e-07 4.441e-16 5.551e-16 1.110e-16 [7,] NaN 1.528e-07 0.000e+00 0.000e+00 0.000e+00 [8,] NaN 2.435e-08 0.000e+00 0.000e+00 0.000e+00 [9,] NaN 3.905e-09 0.000e+00 0.000e+00 0.000e+00 [10,] NaN 6.293e-10 0.000e+00 0.000e+00 0.000e+00 [11,] NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [12,] NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 Warning message: In pbeta(1 - xi, lam - x + cc, x - cc + 1) : NaNs produced > p.rErr.phypDiv(q=k, 1.2*k, 2*k, 1.5*k, abslog=TRUE, logx="x") Warning messages: 1: In pbeta(1 - xi, lam - x + cc, x - cc + 1) : NaNs produced 2: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 21 y values <= 0 omitted from logarithmic plot > > (rErr1215 <- p.rErr.phypDiv(q=k, 1.2*k, 2*k, 1.5*k)) rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peizer [1,] NaN 1.287e-02 3.718e-05 6.630e-05 -2.392e-06 [2,] NaN 1.805e-03 3.092e-07 4.254e-07 -3.749e-09 [3,] NaN 2.646e-04 1.994e-09 2.535e-09 -9.461e-12 [4,] NaN 4.002e-05 1.192e-11 1.460e-11 -2.864e-14 [5,] NaN 6.175e-06 6.894e-14 8.271e-14 0.000e+00 [6,] NaN 9.664e-07 4.441e-16 5.551e-16 1.110e-16 [7,] NaN 1.528e-07 0.000e+00 0.000e+00 0.000e+00 [8,] NaN 2.435e-08 0.000e+00 0.000e+00 0.000e+00 [9,] NaN 3.905e-09 0.000e+00 0.000e+00 0.000e+00 [10,] NaN 6.293e-10 0.000e+00 0.000e+00 0.000e+00 [11,] NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [12,] NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 Warning message: In pbeta(1 - xi, lam - x + cc, x - cc + 1) : NaNs produced > p.rErr.phypDiv(q=k, 1.2*k, 2*k, 1.5*k, abslog=TRUE, logx="x") Warning messages: 1: In pbeta(1 - xi, lam - x + cc, x - cc + 1) : NaNs produced 2: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 21 y values <= 0 omitted from logarithmic plot > > > k <- c(10*(7:9), t(outer(10^(2:5), 1:9))) > x <- round(.8*k); ph.k2k <- phyper(x,1.6*k, 2*k,1.8*k) > (rErr1618 <- p.rErr.phypDiv(q=x, 1.6*k, 2*k, 1.8*k, logx="x")) ## AS151 is really unusable here ! rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peizer [1,] 2.711e-04 0.8682 -2.045e-04 -2.761e-05 -7.714e-07 [2,] 2.235e-04 0.8947 -1.685e-04 -2.271e-05 -5.565e-07 [3,] 1.884e-04 0.9155 -1.419e-04 -1.911e-05 -4.170e-07 [4,] 1.616e-04 0.9320 -1.217e-04 -1.637e-05 -3.220e-07 [5,] 5.864e-05 0.9914 -4.407e-05 -5.902e-06 -5.847e-08 [6,] 3.229e-05 0.9988 -2.425e-05 -3.243e-06 -2.147e-08 [7,] 2.111e-05 0.9998 -1.585e-05 -2.118e-06 -1.053e-08 [8,] 1.518e-05 1.0000 -1.139e-05 -1.522e-06 -6.054e-09 [9,] 1.158e-05 1.0000 -8.694e-06 -1.161e-06 -3.850e-09 [10,] 9.217e-06 1.0000 -6.917e-06 -9.234e-07 -2.628e-09 [11,] 7.560e-06 1.0000 -5.673e-06 -7.572e-07 -1.887e-09 [12,] 6.347e-06 1.0000 -4.762e-06 -6.356e-07 -1.411e-09 [13,] 5.427e-06 1.0000 -4.072e-06 -5.434e-07 -1.085e-09 [14,] 1.934e-06 1.0000 -1.451e-06 -1.935e-07 -1.944e-10 [15,] 1.056e-06 1.0000 -7.924e-07 -1.057e-07 -7.027e-11 [16,] 6.876e-07 1.0000 -5.158e-07 -6.878e-08 -3.014e-11 [17,] 4.927e-07 1.0000 -3.696e-07 -4.928e-08 3.679e-12 [18,] 3.752e-07 1.0000 -2.814e-07 -3.753e-08 -2.016e-11 [19,] 2.980e-07 1.0000 -2.235e-07 -2.981e-08 1.187e-12 [20,] 2.441e-07 1.0000 -1.831e-07 -2.441e-08 -1.793e-11 [21,] 2.047e-07 1.0000 -1.535e-07 -2.047e-08 6.614e-11 [22,] 1.748e-07 1.0000 -1.311e-07 -1.748e-08 -1.525e-10 [23,] 6.196e-08 1.0000 -4.647e-08 -6.197e-09 4.220e-10 [24,] 3.377e-08 1.0000 -2.532e-08 -3.377e-09 -4.321e-10 [25,] 2.195e-08 1.0000 -1.646e-08 -2.195e-09 -5.324e-10 [26,] 1.571e-08 1.0000 -1.178e-08 -1.571e-09 -9.303e-10 [27,] 1.195e-08 1.0000 -8.966e-09 -1.196e-09 9.785e-10 [28,] 9.489e-09 1.0000 -7.117e-09 -9.490e-10 -1.233e-09 [29,] 7.769e-09 1.0000 -5.826e-09 -7.768e-10 -2.975e-14 [30,] 6.512e-09 1.0000 -4.884e-09 -6.512e-10 1.799e-09 [31,] 5.560e-09 1.0000 -4.170e-09 -5.561e-10 6.661e-15 [32,] 1.967e-09 1.0000 -1.476e-09 -1.968e-10 1.491e-09 [33,] 1.071e-09 1.0000 -8.035e-10 -1.071e-10 -1.096e-08 [34,] 6.960e-10 1.0000 -5.220e-10 -6.961e-11 2.731e-14 [35,] 4.982e-10 1.0000 -3.736e-10 -4.982e-11 2.949e-08 [36,] 3.789e-10 1.0000 -2.843e-10 -3.785e-11 3.877e-08 [37,] 3.007e-10 1.0000 -2.256e-10 -3.001e-11 -4.886e-08 [38,] 2.462e-10 1.0000 -1.846e-10 -2.463e-11 7.361e-14 [39,] 2.062e-10 1.0000 -1.547e-10 -2.063e-11 4.818e-14 > p.rErr.phypDiv(q=x, 1.6*k, 2*k, 1.8*k, abslog=TRUE, logx="x") > ## interestingly, Peizer suffers from some erratic behavior for large q > ## --> the L has 4 log(.) terms, each of which has argument ~= 1 <--> could use log1p(.) or ??? > > k <- k1 # the old set > par(mfrow=c(2,1)) > for(Reps in c(0,1)) { + options(rErr.eps = Reps) ## 0: always RELATIVE error; 1: always ABSOLUTE + x <- round(.6*k); ph.k2k <- phyper(x,1.6*k, 2*k,1.8*k) + ## ^^^ + print(ph.mat <- + cbind(k=k, phyper= ph.k2k, + rERR.Ibeta= rErr(phyperIbeta (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.as151= rErr(phyperApprAS152 (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.1mol = rErr(phyper1molenaar (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.2mol = rErr(phyper2molenaar (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.Peiz = rErr(phyperPeizer (x,1.6*k,2*k,1.8*k) , ph.k2k)) + ) + ## The two molenaar's ``break down''; Peizer remains decent! + Err.phyp <- abs(ph.mat[,-(1:2)]) + matplot(k, Err.phyp, + ylim= if((ee <- .Options$rErr.eps) > 0) c(1e-200,1), + main="|relE.{phyper(x = 0.6k, 1.6k, 2k, 1.8k)}|", + type='o', log='xy') + } k phyper rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peiz [1,] 10 1.573e-01 -0.023303 0.6444 2.162e-02 4.521e-04 6.654e-04 [2,] 20 4.810e-02 -0.003428 0.8714 1.667e-02 -6.849e-03 4.014e-04 [3,] 30 1.633e-02 0.020289 0.9530 9.682e-03 -1.450e-02 2.899e-04 [4,] 40 5.805e-03 0.044219 0.9827 2.331e-03 -2.220e-02 2.303e-04 [5,] 50 2.118e-03 0.067790 0.9936 -5.154e-03 -2.994e-02 1.935e-04 [6,] 60 7.863e-04 0.090861 0.9977 -1.272e-02 -3.773e-02 1.686e-04 [7,] 70 2.954e-04 0.113394 0.9991 -2.034e-02 -4.556e-02 1.505e-04 [8,] 80 1.120e-04 0.135381 0.9997 -2.802e-02 -5.344e-02 1.369e-04 [9,] 90 4.274e-05 0.156828 0.9999 -3.576e-02 -6.136e-02 1.263e-04 [10,] 100 1.640e-05 0.177745 1.0000 -4.355e-02 -6.935e-02 1.178e-04 [11,] 200 1.340e-09 0.360344 1.0000 -1.246e-01 -1.524e-01 7.908e-05 [12,] 300 1.253e-13 0.502335 1.0000 -2.118e-01 -2.417e-01 6.610e-05 [13,] 400 1.239e-17 0.612786 1.0000 -3.057e-01 -3.379e-01 5.959e-05 [14,] 500 1.264e-21 0.698716 1.0000 -4.069e-01 -4.416e-01 5.568e-05 [15,] 600 1.315e-25 0.765574 1.0000 -5.159e-01 -5.533e-01 5.308e-05 [16,] 700 1.387e-29 0.817594 1.0000 -6.334e-01 -6.736e-01 5.122e-05 [17,] 800 1.478e-33 0.858069 1.0000 -7.599e-01 -8.033e-01 4.982e-05 [18,] 900 1.588e-37 0.889563 1.0000 -8.962e-01 -9.430e-01 4.873e-05 [19,] 1000 1.716e-41 0.914068 1.0000 -1.043e+00 -1.093e+00 4.786e-05 [20,] 2000 4.451e-81 0.993008 1.0000 -3.308e+00 -3.414e+00 4.395e-05 [21,] 3000 1.332e-120 0.999431 1.0000 -8.085e+00 -8.309e+00 4.265e-05 [22,] 4000 4.225e-160 0.999954 1.0000 -1.816e+01 -1.863e+01 4.199e-05 [23,] 5000 1.384e-199 0.999996 1.0000 -3.940e+01 -4.039e+01 4.160e-05 [24,] 6000 4.628e-239 1.000000 1.0000 -8.418e+01 -8.628e+01 4.134e-05 [25,] 7000 1.569e-278 1.000000 1.0000 -1.786e+02 -1.830e+02 4.115e-05 [26,] 8000 5.376e-318 1.000000 1.0000 1.000e+00 1.000e+00 1.000e+00 [27,] 9000 0.000e+00 NaN NaN NaN NaN NaN k phyper rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peiz [1,] 10 1.573e-01 -3.665e-03 1.013e-01 3.400e-03 7.110e-05 1.046e-04 [2,] 20 4.810e-02 -1.649e-04 4.191e-02 8.017e-04 -3.294e-04 1.931e-05 [3,] 30 1.633e-02 3.313e-04 1.556e-02 1.581e-04 -2.367e-04 4.734e-06 [4,] 40 5.805e-03 2.567e-04 5.705e-03 1.353e-05 -1.289e-04 1.337e-06 [5,] 50 2.118e-03 1.436e-04 2.105e-03 -1.092e-05 -6.343e-05 4.099e-07 [6,] 60 7.863e-04 7.144e-05 7.844e-04 -9.998e-06 -2.966e-05 1.325e-07 [7,] 70 2.954e-04 3.349e-05 2.951e-04 -6.008e-06 -1.346e-05 4.447e-08 [8,] 80 1.120e-04 1.516e-05 1.119e-04 -3.137e-06 -5.983e-06 1.533e-08 [9,] 90 4.274e-05 6.702e-06 4.273e-05 -1.528e-06 -2.622e-06 5.398e-09 [10,] 100 1.640e-05 2.915e-06 1.640e-05 -7.143e-07 -1.137e-06 1.932e-09 [11,] 200 1.340e-09 4.830e-10 1.340e-09 -1.670e-10 -2.042e-10 1.060e-13 [12,] 300 1.253e-13 6.293e-14 1.253e-13 -2.653e-14 -3.028e-14 8.280e-18 [13,] 400 1.239e-17 7.591e-18 1.239e-17 -3.787e-18 -4.186e-18 7.382e-22 [14,] 500 1.264e-21 8.829e-22 1.264e-21 -5.142e-22 -5.580e-22 7.036e-26 [15,] 600 1.315e-25 1.007e-25 1.315e-25 -6.784e-26 -7.275e-26 6.979e-30 [16,] 700 1.387e-29 1.134e-29 1.387e-29 -8.786e-30 -9.345e-30 7.105e-34 [17,] 800 1.478e-33 1.269e-33 1.478e-33 -1.123e-33 -1.188e-33 7.365e-38 [18,] 900 1.588e-37 1.412e-37 1.588e-37 -1.423e-37 -1.497e-37 7.738e-42 [19,] 1000 1.716e-41 1.568e-41 1.716e-41 -1.790e-41 -1.876e-41 8.212e-46 [20,] 2000 4.451e-81 4.420e-81 4.451e-81 -1.473e-80 -1.520e-80 1.956e-85 [21,] 3000 1.332e-120 1.331e-120 1.332e-120 -1.077e-119 -1.106e-119 5.679e-125 [22,] 4000 4.225e-160 4.225e-160 4.225e-160 -7.671e-159 -7.870e-159 1.774e-164 [23,] 5000 1.384e-199 1.384e-199 1.384e-199 -5.453e-198 -5.590e-198 5.758e-204 [24,] 6000 4.628e-239 4.628e-239 4.628e-239 -3.896e-237 -3.993e-237 1.913e-243 [25,] 7000 1.569e-278 1.569e-278 1.569e-278 -2.803e-276 -2.872e-276 6.458e-283 [26,] 8000 5.376e-318 5.376e-318 5.376e-318 5.376e-318 5.376e-318 5.376e-318 [27,] 9000 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 5 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > > > hp1 <- list(n = 10, n1 = 12, n2 = 2) > > for(h.nr in 1:1) { + attach(hplist <- get(paste("hp",h.nr, sep='')) ) # defining n, n1, n2 + cat("\n"); print(unlist(hplist)) + N <- n1 + n2 + x <- 0:max(0,min(n1, n2 - n)) + d <- dhyper(x, n1, n2, n) + names(d) <- as.character(x) + print(d) + } n n1 n2 10 12 2 0 0 > > > ### Binomial Approximation(s) ================================================= > > ph.m <- cbind(ph = phyper (0:5, 5,15, 7), + pM1 = phyperBinMolenaar.1(0:5, 5,15, 7), + pM2 = phyperBinMolenaar.2(0:5, 5,15, 7), + pM3 = phyperBinMolenaar.3(0:5, 5,15, 7), + pM4 = phyperBinMolenaar.4(0:5, 5,15, 7)) > > cc <- adjustcolor(1:(1+ncol(ph.m)), 0.5); ppch <- as.character(0:4) > matplot(0:5, ph.m, type = "b", lwd=2, col=cc, pch=ppch, + main = "all 4 phyper() binomial approximations via Molenaar's") > legend("right", legend = c("phyper", paste("ph.appr.Mol", 1:4)), + lwd=2, col=cc, lty=1:5, pch=ppch, bty="n") > > > ## relative errors : > ph.err <- 1 - ph.m[,-1] / ph.m[,"ph"] > > matplot(0:5, ph.err, type = "b", lwd=2, col=cc, pch=ppch, + main = "rel.Err. of the 4 phyper() binom.Molenaar's approximations") > legend("topright", legend = paste("ph.appr.Mol", 1:4), + lwd=2, col=cc[-1], lty=2:5, pch=ppch[-1], bty="n") > ## quite interesting: if we take the *mean* of approx. 1 and 2 -- should become good? > > ## have rErr() above {but does not work with *matrix* ?} > > rErr.Mol.bin <- function(m,n,k, q = .suppHyper(m,n,k), lower.tail=TRUE, log.p=FALSE) { + ph <- phyper (q, m=m, n=n, k=k, lower.tail=lower.tail, log.p=log.p) + phM <- phyperAllBinM(m=m, n=n, k=k, q=q, lower.tail=lower.tail, log.p=log.p) + apply(phM, 2, rErr, true = ph) + } > ## Plotting of the rel.error: > p.rErr.Mol.bin <- function(m,n,k, q = .suppHyper(m,n,k), abslog = FALSE, + lower.tail=TRUE, log.p=FALSE) + { + rE <- rErr.Mol.bin(m=m, n=n, k=k, q=q, lower.tail=lower.tail, log.p=log.p) + cc <- adjustcolor(1:4, 0.5) + ppch <- as.character(1:4) + matplot(q, if(abslog) abs(rE) else rE, type = "b", lwd=2, col=cc, pch=ppch, + log = if(abslog) "y" else "", + main = sprintf("phyper(*, m = %d, n = %d, k = %d)", m,n,k)) + mtext(paste(if(abslog) "|rel.Err|" else "rel.Err.", + "of the 4 phyper() binom.Molenaar's approximations")) + if(!abslog) abline(h=0, col=adjustcolor("gray", 0.5), lty=2) + legend("topright", legend = paste("phyp.Mol", 1:4), + lwd=2, col=cc, lty=1:4, pch=ppch, bty="n") + invisible(rE) + } > > p.rErr.Mol.bin (5,15, 7) # same plot, as the "manual" one above > p.rErr.Mol.bin (70,100, 7)# here, 1 and 4 are good > p.rErr.Mol.bin (70,100, 20)# 4 (and 1) ((maybe their *mean* !) > p.rErr.Mol.bin (70,100, 20, abslog=TRUE)# 4 (and 1) Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 2 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > p.rErr.Mol.bin (70,100, 50)# 4 > p.rErr.Mol.bin (70,100, 50, abslog=TRUE) -> rE.71c50 Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 28 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 8 y values <= 0 omitted from logarithmic plot > (p.rErr.Mol.bin (70,100, 70,abslog=TRUE) -> rE.71c70)# 1 & 2; but *none* is good for small q pM1 pM2 pM3 pM4 [1,] -5.586e-25 -5.586e-25 -1.629e-24 -1.629e-24 [2,] -6.027e-23 -6.027e-23 -2.182e-22 -2.182e-22 [3,] -2.767e-21 -2.767e-21 -1.378e-20 -1.378e-20 [4,] -6.336e-20 -6.336e-20 -5.483e-19 -5.483e-19 [5,] -3.145e-19 -3.145e-19 -1.546e-17 -1.546e-17 [6,] 2.806e-17 2.806e-17 -3.301e-16 -3.301e-16 [7,] 1.113e-15 1.113e-15 -5.562e-15 -5.562e-15 [8,] 2.466e-14 2.466e-14 -7.621e-14 -7.621e-14 [9,] 3.934e-13 3.934e-13 -8.678e-13 -8.678e-13 [10,] 4.898e-12 4.898e-12 -8.355e-12 -8.355e-12 [11,] 4.953e-11 4.953e-11 -6.897e-11 -6.897e-11 [12,] 4.166e-10 4.166e-10 -4.936e-10 -4.936e-10 [13,] 2.964e-09 2.964e-09 -3.092e-09 -3.092e-09 [14,] 1.805e-08 1.805e-08 -1.708e-08 -1.708e-08 [15,] 9.500e-08 9.500e-08 -8.372e-08 -8.372e-08 [16,] 4.354e-07 4.354e-07 -3.657e-07 -3.657e-07 [17,] 1.748e-06 1.748e-06 -1.428e-06 -1.428e-06 [18,] 6.182e-06 6.182e-06 -5.002e-06 -5.002e-06 [19,] 1.934e-05 1.934e-05 -1.573e-05 -1.573e-05 [20,] 5.376e-05 5.376e-05 -4.443e-05 -4.443e-05 [21,] 1.332e-04 1.332e-04 -1.129e-04 -1.129e-04 [22,] 2.951e-04 2.951e-04 -2.580e-04 -2.580e-04 [23,] 5.866e-04 5.866e-04 -5.305e-04 -5.305e-04 [24,] 1.049e-03 1.049e-03 -9.812e-04 -9.812e-04 [25,] 1.691e-03 1.691e-03 -1.632e-03 -1.632e-03 [26,] 2.467e-03 2.467e-03 -2.442e-03 -2.442e-03 [27,] 3.262e-03 3.262e-03 -3.284e-03 -3.284e-03 [28,] 3.920e-03 3.920e-03 -3.967e-03 -3.967e-03 [29,] 4.292e-03 4.292e-03 -4.302e-03 -4.302e-03 [30,] 4.289e-03 4.289e-03 -4.183e-03 -4.183e-03 [31,] 3.918e-03 3.918e-03 -3.641e-03 -3.641e-03 [32,] 3.275e-03 3.275e-03 -2.830e-03 -2.830e-03 [33,] 2.507e-03 2.507e-03 -1.958e-03 -1.958e-03 [34,] 1.757e-03 1.757e-03 -1.199e-03 -1.199e-03 [35,] 1.127e-03 1.127e-03 -6.457e-04 -6.457e-04 [36,] 6.603e-04 6.603e-04 -3.016e-04 -3.016e-04 [37,] 3.531e-04 3.531e-04 -1.193e-04 -1.193e-04 [38,] 1.720e-04 1.720e-04 -3.785e-05 -3.785e-05 [39,] 7.615e-05 7.615e-05 -8.131e-06 -8.131e-06 [40,] 3.058e-05 3.058e-05 -4.888e-09 -4.888e-09 [41,] 1.111e-05 1.111e-05 1.093e-06 1.093e-06 [42,] 3.647e-06 3.647e-06 6.860e-07 6.860e-07 [43,] 1.078e-06 1.078e-06 2.899e-07 2.899e-07 [44,] 2.867e-07 2.867e-07 9.772e-08 9.772e-08 [45,] 6.836e-08 6.836e-08 2.766e-08 2.766e-08 [46,] 1.459e-08 1.459e-08 6.721e-09 6.721e-09 [47,] 2.781e-09 2.781e-09 1.417e-09 1.417e-09 [48,] 4.722e-10 4.722e-10 2.607e-10 2.607e-10 [49,] 7.121e-11 7.121e-11 4.195e-11 4.195e-11 [50,] 9.515e-12 9.515e-12 5.909e-12 5.909e-12 [51,] 1.123e-12 1.123e-12 7.280e-13 7.280e-13 [52,] 1.168e-13 1.168e-13 7.838e-14 7.838e-14 [53,] 1.066e-14 1.066e-14 7.327e-15 7.327e-15 [54,] 7.772e-16 7.772e-16 5.551e-16 5.551e-16 [55,] 1.110e-16 1.110e-16 0.000e+00 0.000e+00 [56,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [57,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [58,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [59,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [60,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [61,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [62,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [63,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [64,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [65,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [66,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [67,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [68,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [69,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [70,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [71,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 66 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 16 y values <= 0 omitted from logarithmic plot > par(new=TRUE) > plot(0:70, phyper(0:70, 70,100, 70), type = "l", col="blue", ann=FALSE,axes=FALSE) > > ## rel error on log-scale ... really nonsense > (p.rErr.Mol.bin (70,100, 70, log.p=TRUE))# 1 & 2; but *none* is good for small q pM1 pM2 pM3 pM4 [1,] 2.166e-03 2.166e-03 5.722e-03 5.722e-03 [2,] 1.652e-03 1.652e-03 5.439e-03 5.439e-03 [3,] 1.128e-03 1.128e-03 5.131e-03 5.131e-03 [4,] 6.085e-04 6.085e-04 4.827e-03 4.827e-03 [5,] 1.003e-04 1.003e-04 4.538e-03 4.538e-03 [6,] -3.933e-04 -3.933e-04 4.272e-03 4.272e-03 [7,] -8.706e-04 -8.706e-04 4.033e-03 4.033e-03 [8,] -1.331e-03 -1.331e-03 3.824e-03 3.824e-03 [9,] -1.773e-03 -1.773e-03 3.649e-03 3.649e-03 [10,] -2.198e-03 -2.198e-03 3.510e-03 3.510e-03 [11,] -2.606e-03 -2.606e-03 3.409e-03 3.409e-03 [12,] -2.998e-03 -2.998e-03 3.348e-03 3.348e-03 [13,] -3.376e-03 -3.376e-03 3.330e-03 3.330e-03 [14,] -3.740e-03 -3.740e-03 3.357e-03 3.357e-03 [15,] -4.094e-03 -4.094e-03 3.433e-03 3.433e-03 [16,] -4.439e-03 -4.439e-03 3.559e-03 3.559e-03 [17,] -4.779e-03 -4.779e-03 3.739e-03 3.739e-03 [18,] -5.118e-03 -5.118e-03 3.977e-03 3.977e-03 [19,] -5.462e-03 -5.462e-03 4.278e-03 4.278e-03 [20,] -5.817e-03 -5.817e-03 4.646e-03 4.646e-03 [21,] -6.192e-03 -6.192e-03 5.087e-03 5.087e-03 [22,] -6.596e-03 -6.596e-03 5.606e-03 5.606e-03 [23,] -7.044e-03 -7.044e-03 6.212e-03 6.212e-03 [24,] -7.552e-03 -7.552e-03 6.910e-03 6.910e-03 [25,] -8.144e-03 -8.144e-03 7.709e-03 7.709e-03 [26,] -8.847e-03 -8.847e-03 8.614e-03 8.614e-03 [27,] -9.702e-03 -9.702e-03 9.631e-03 9.631e-03 [28,] -1.076e-02 -1.076e-02 1.076e-02 1.076e-02 [29,] 9.377e-03 9.377e-03 -9.312e-03 -9.312e-03 [30,] 7.354e-03 7.354e-03 -7.121e-03 -7.121e-03 [31,] 5.592e-03 5.592e-03 -5.169e-03 -5.169e-03 [32,] 4.094e-03 4.094e-03 -3.524e-03 -3.524e-03 [33,] 2.860e-03 2.860e-03 -2.228e-03 -2.228e-03 [34,] 1.890e-03 1.890e-03 -1.288e-03 -1.288e-03 [35,] 1.170e-03 1.170e-03 -6.698e-04 -6.698e-04 [36,] 6.721e-04 6.721e-04 -3.068e-04 -3.068e-04 [37,] 3.559e-04 3.559e-04 -1.202e-04 -1.202e-04 [38,] 1.725e-04 1.725e-04 -3.796e-05 -3.796e-05 [39,] 7.623e-05 7.623e-05 -8.140e-06 -8.140e-06 [40,] 3.059e-05 3.059e-05 -4.890e-09 -4.890e-09 [41,] 1.111e-05 1.111e-05 1.093e-06 1.093e-06 [42,] 3.647e-06 3.647e-06 6.860e-07 6.860e-07 [43,] 1.078e-06 1.078e-06 2.899e-07 2.899e-07 [44,] 2.867e-07 2.867e-07 9.772e-08 9.772e-08 [45,] 6.836e-08 6.836e-08 2.766e-08 2.766e-08 [46,] 1.459e-08 1.459e-08 6.721e-09 6.721e-09 [47,] 2.781e-09 2.781e-09 1.417e-09 1.417e-09 [48,] 4.722e-10 4.722e-10 2.607e-10 2.607e-10 [49,] 7.121e-11 7.121e-11 4.195e-11 4.195e-11 [50,] 9.515e-12 9.515e-12 5.909e-12 5.909e-12 [51,] 1.123e-12 1.123e-12 7.280e-13 7.280e-13 [52,] 1.167e-13 1.167e-13 7.835e-14 7.835e-14 [53,] 1.064e-14 1.064e-14 7.350e-15 7.350e-15 [54,] 8.472e-16 8.472e-16 5.993e-16 5.993e-16 [55,] 5.873e-17 5.873e-17 4.232e-17 4.232e-17 [56,] 3.526e-18 3.526e-18 2.578e-18 2.578e-18 [57,] 1.824e-19 1.824e-19 1.348e-19 1.348e-19 [58,] 8.081e-21 8.081e-21 6.017e-21 6.017e-21 [59,] 3.047e-22 3.047e-22 2.279e-22 2.279e-22 [60,] 9.706e-24 9.706e-24 7.273e-24 7.273e-24 [61,] 2.589e-25 2.589e-25 1.939e-25 1.939e-25 [62,] 5.727e-27 5.727e-27 4.273e-27 4.273e-27 [63,] 1.037e-28 1.037e-28 7.695e-29 7.695e-29 [64,] 1.513e-30 1.513e-30 1.115e-30 1.115e-30 [65,] 1.745e-32 1.745e-32 1.273e-32 1.273e-32 [66,] 1.544e-34 1.544e-34 1.115e-34 1.115e-34 [67,] 1.008e-36 1.008e-36 7.193e-37 7.193e-37 [68,] 4.545e-39 4.545e-39 3.203e-39 3.203e-39 [69,] 1.258e-41 1.258e-41 8.749e-42 8.749e-42 [70,] 1.601e-44 1.601e-44 1.098e-44 1.098e-44 [71,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 > ## however, *again* the *mean* of 1|2 with 3|4 seems fine > > > ## Total Variation error -- what Kuensch(1998) used for binom() -- hyper() approx. > TVerr <- function(approx, true) { + d <- true - approx + max(sum( d[d > 0]), + sum(- d[d < 0])) + } > supErr <- function(approx, true) max(abs(true - approx)) > > TVerr.bin <- function(m,n,k, q = .suppHyper(m,n,k), lower.tail=TRUE, log.p=FALSE) { + ph <- phyper (q, m=m, n=n, k=k, lower.tail=lower.tail, log.p=log.p) + phM <- phyperAllBin(m=m, n=n, k=k, q=q, lower.tail=lower.tail, log.p=log.p) + apply(phM, 2, TVerr, true = ph) + } > supErr.binM <- function(m,n,k, q = .suppHyper(m,n,k), lower.tail=TRUE, log.p=FALSE) { + ph <- phyper (q, m=m, n=n, k=k, lower.tail=lower.tail, log.p=log.p) + phM <- phyperAllBinM(m=m, n=n, k=k, q=q, lower.tail=lower.tail, log.p=log.p) + apply(phM, 2, supErr, true = ph) + } > > (ph.TV <- apply(ph.m[,-1], 2, TVerr, true = ph.m[,1])) pM1 pM2 pM3 pM4 0.02608 0.01202 0.01087 0.02116 > ## in all three cases, Molenaar is clearly better than "simple" > TVerr.bin(5,15,7) p1 p2 p3 p4 pM1 pM2 pM3 pM4 0.08959 0.05560 0.05560 0.08959 0.02608 0.01202 0.01087 0.02116 > TVerr.bin(1000,2000,7) # 1 & 4 are very accurate p1 p2 p3 p4 pM1 pM2 pM3 pM4 5.125e-04 1.147e-01 1.147e-01 5.125e-04 9.517e-07 2.790e-02 1.693e-02 8.649e-07 > TVerr.bin(1000,2000, 1500)# not very good p1 p2 p3 p4 pM1 pM2 pM3 pM4 2.13342 1.15743 1.15743 2.13342 0.06290 0.02311 0.02311 0.05688 > > TVerr.bin(100000, 200000, 150000) p1 p2 p3 p4 pM1 pM2 pM3 pM4 21.33334 11.57508 11.57508 21.33334 0.27892 0.09201 0.09201 0.27369 > ## but TV error really *grows* with the support of the distribution > supErr.binM(100000, 200000, 150000) pM1 pM2 pM3 pM4 0.0012773 0.0004212 0.0004212 0.0012586 > supErr.binM(1000000, 2000000, 150000) pM1 pM2 pM3 pM4 2.446e-05 4.171e-04 4.019e-04 2.434e-05 > supErr.binM(1e6, 1e9, 150000) pM1 pM2 pM3 pM4 3.277e-10 1.456e-08 1.247e-08 2.808e-10 > round(-log10(supErr.binM(1e6, 1e9, 1e7)), 2) #=> correct #{digits} pM1 pM2 pM3 pM4 5.98 7.98 7.99 5.99 > ## pM1 pM2 pM3 pM4 > ## 5.98 7.98 7.99 5.99 -- m = 1e6 < k = 1e7 ==> approx. with size = m are better > system.time(sE <- supErr.binM(1e7, 1e9, 1e7)) ## 20 sec. elapsed! ==> very slow !! user system elapsed 29.37 2.35 31.73 > round(-log10(sE), 2) pM1 pM2 pM3 pM4 6.00 6.00 6.01 6.01 > ## pM1 pM2 pM3 pM4 > ## 6.00 6.00 6.01 6.01 --- m = k ==> all are the same > > system.time(sE <- supErr.binM(1e7, 1e9, 100))# fast user system elapsed 0 0 0 > round(-log10(sE), 2) pM1 pM2 pM3 pM4 15.35 5.22 5.21 14.65 > ## pM1 pM2 pM3 pM4 > ## 15.35 5.22 5.21 14.65 k << m ==> 1 & 4 are much better {why '1' better than '4' ?} > > > > > ### Here, we currently only explore normal approx: > > k <- c(10*(1:9),100*(1:9),1000*(1:9)) > > options(digits=4) > ph.k2k <- phyper(k,2*k,2*k,2*k) # rel.err ~ 10^-7 > cbind(k=k, phyper= ph.k2k, + rERR.Ibeta= rErr(phyperIbeta (k,2*k,2*k,2*k) , ph.k2k), + rERR.as151= rErr(phyperApprAS152 (k,2*k,2*k,2*k) , ph.k2k), + rERR.1mol = rErr(phyper1molenaar (k,2*k,2*k,2*k) , ph.k2k), + rERR.2mol = rErr(phyper2molenaar (k,2*k,2*k,2*k) , ph.k2k), + rERR.Peizer=rErr(phyperPeizer (k,2*k,2*k,2*k) , ph.k2k)) k phyper rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peizer [1,] 10 0.6238 2.240e-03 1.240e-03 -1.799e-03 -2.708e-04 -3.725e-05 [2,] 20 0.5884 8.141e-04 4.511e-04 -6.433e-04 -9.445e-05 -6.803e-06 [3,] 30 0.5724 4.472e-04 2.480e-04 -3.515e-04 -5.115e-05 -2.496e-06 [4,] 40 0.5628 2.917e-04 1.618e-04 -2.287e-04 -3.313e-05 -1.222e-06 [5,] 50 0.5562 2.093e-04 1.162e-04 -1.638e-04 -2.367e-05 -7.019e-07 [6,] 60 0.5513 1.595e-04 8.853e-05 -1.247e-04 -1.799e-05 -4.460e-07 [7,] 70 0.5476 1.268e-04 7.036e-05 -9.903e-05 -1.426e-05 -3.038e-07 [8,] 80 0.5445 1.038e-04 5.765e-05 -8.109e-05 -1.167e-05 -2.178e-07 [9,] 90 0.5420 8.709e-05 4.835e-05 -6.798e-05 -9.772e-06 -1.624e-07 [10,] 100 0.5398 7.440e-05 4.131e-05 -5.805e-05 -8.340e-06 -1.249e-07 [11,] 200 0.5282 2.638e-05 1.465e-05 -2.055e-05 -2.944e-06 -2.215e-08 [12,] 300 0.5230 1.437e-05 7.982e-06 -1.119e-05 -1.601e-06 -8.047e-09 [13,] 400 0.5199 9.338e-06 5.187e-06 -7.268e-06 -1.040e-06 -3.922e-09 [14,] 500 0.5178 6.683e-06 3.712e-06 -5.201e-06 -7.439e-07 -2.245e-09 [15,] 600 0.5163 5.085e-06 2.825e-06 -3.957e-06 -5.658e-07 -1.423e-09 [16,] 700 0.5151 4.036e-06 2.242e-06 -3.140e-06 -4.490e-07 -9.688e-10 [17,] 800 0.5141 3.304e-06 1.835e-06 -2.570e-06 -3.675e-07 -6.959e-10 [18,] 900 0.5133 2.769e-06 1.538e-06 -2.154e-06 -3.079e-07 -5.193e-10 [19,] 1000 0.5126 2.364e-06 1.313e-06 -1.839e-06 -2.629e-07 -3.973e-10 [20,] 2000 0.5089 8.361e-07 4.645e-07 -6.504e-07 -9.294e-08 -7.817e-11 [21,] 3000 0.5073 4.551e-07 2.529e-07 -3.540e-07 -5.059e-08 -2.550e-11 [22,] 4000 0.5063 2.956e-07 1.642e-07 -2.300e-07 -3.286e-08 9.986e-12 [23,] 5000 0.5056 2.115e-07 1.175e-07 -1.645e-07 -2.351e-08 -7.108e-12 [24,] 6000 0.5052 1.609e-07 8.941e-08 -1.252e-07 -1.788e-08 -4.567e-11 [25,] 7000 0.5048 1.277e-07 7.095e-08 -9.934e-08 -1.419e-08 4.881e-11 [26,] 8000 0.5045 1.045e-07 5.807e-08 -8.131e-08 -1.162e-08 -6.557e-11 [27,] 9000 0.5042 8.760e-08 4.867e-08 -6.814e-08 -9.735e-09 7.399e-11 > > ## Here, the Ibeta fails (NaN) ; moleaar's are both very good : > ph.k2k <- phyper(k, 1.2*k, 2*k,1.5*k) > cbind(k=k, phyper= ph.k2k, + rERR.Ibeta= rErr(phyperIbeta (k,1.2*k,2*k,1.5*k) , ph.k2k), + rERR.as151= rErr(phyperApprAS152 (k,1.2*k,2*k,1.5*k) , ph.k2k), + rERR.1mol = rErr(phyper1molenaar (k,1.2*k,2*k,1.5*k) , ph.k2k), + rERR.2mol = rErr(phyper2molenaar (k,1.2*k,2*k,1.5*k) , ph.k2k), + rERR.Peizer=rErr(phyperPeizer (k,1.2*k,2*k,1.5*k) , ph.k2k)) k phyper rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peizer [1,] 10 0.9999 NaN 1.287e-02 3.718e-05 6.629e-05 -2.391e-06 [2,] 20 1.0000 NaN 1.805e-03 3.092e-07 4.254e-07 -3.749e-09 [3,] 30 1.0000 NaN 2.646e-04 1.994e-09 2.535e-09 -9.461e-12 [4,] 40 1.0000 NaN 4.002e-05 1.192e-11 1.460e-11 -2.853e-14 [5,] 50 1.0000 NaN 6.175e-06 6.894e-14 8.271e-14 0.000e+00 [6,] 60 1.0000 NaN 9.664e-07 4.441e-16 5.551e-16 1.110e-16 [7,] 70 1.0000 NaN 1.528e-07 0.000e+00 0.000e+00 0.000e+00 [8,] 80 1.0000 NaN 2.435e-08 0.000e+00 0.000e+00 0.000e+00 [9,] 90 1.0000 NaN 3.905e-09 0.000e+00 0.000e+00 0.000e+00 [10,] 100 1.0000 NaN 6.293e-10 0.000e+00 0.000e+00 0.000e+00 [11,] 200 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [12,] 300 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [13,] 400 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [14,] 500 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [15,] 600 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [16,] 700 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [17,] 800 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [18,] 900 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [19,] 1000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [20,] 2000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [21,] 3000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [22,] 4000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [23,] 5000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [24,] 6000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [25,] 7000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [26,] 8000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 [27,] 9000 1.0000 NaN 0.000e+00 0.000e+00 0.000e+00 0.000e+00 Warning message: In pbeta(1 - xi, lam - x + cc, x - cc + 1) : NaNs produced > > x <- round(.8*k); ph.k2k <- phyper(x,1.6*k, 2*k,1.8*k) > (ph.mat <- + cbind(k=k, phyper= ph.k2k, + rERR.Ibeta= rErr(phyperIbeta (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.as151= rErr(phyperApprAS152 (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.1mol = rErr(phyper1molenaar (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.2mol = rErr(phyper2molenaar (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.Peiz = rErr(phyperPeizer (x,1.6*k,2*k,1.8*k) , ph.k2k)) + ) k phyper rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peiz [1,] 10 0.6310 2.621e-03 0.2331 -2.048e-03 -2.954e-04 -5.166e-05 [2,] 20 0.5936 9.561e-04 0.3212 -7.320e-04 -1.017e-04 -9.477e-06 [3,] 30 0.5767 5.258e-04 0.3764 -3.998e-04 -5.484e-05 -3.482e-06 [4,] 40 0.5666 3.432e-04 0.4145 -2.601e-04 -3.543e-05 -1.707e-06 [5,] 50 0.5596 2.464e-04 0.4421 -1.863e-04 -2.527e-05 -9.805e-07 [6,] 60 0.5544 1.878e-04 0.4625 -1.418e-04 -1.918e-05 -6.231e-07 [7,] 70 0.5504 1.492e-04 0.4779 -1.126e-04 -1.520e-05 -4.246e-07 [8,] 80 0.5472 1.223e-04 0.4896 -9.218e-05 -1.243e-05 -3.045e-07 [9,] 90 0.5445 1.026e-04 0.4985 -7.727e-05 -1.040e-05 -2.271e-07 [10,] 100 0.5422 8.763e-05 0.5053 -6.599e-05 -8.876e-06 -1.746e-07 [11,] 200 0.5299 3.107e-05 0.5253 -2.335e-05 -3.128e-06 -3.098e-08 [12,] 300 0.5244 1.693e-05 0.5238 -1.272e-05 -1.700e-06 -1.126e-08 [13,] 400 0.5211 1.100e-05 0.5211 -8.261e-06 -1.104e-06 -5.487e-09 [14,] 500 0.5189 7.875e-06 0.5189 -5.911e-06 -7.896e-07 -3.142e-09 [15,] 600 0.5173 5.992e-06 0.5173 -4.497e-06 -6.005e-07 -1.991e-09 [16,] 700 0.5160 4.756e-06 0.5160 -3.569e-06 -4.765e-07 -1.356e-09 [17,] 800 0.5150 3.893e-06 0.5150 -2.921e-06 -3.899e-07 -9.719e-10 [18,] 900 0.5141 3.263e-06 0.5141 -2.448e-06 -3.268e-07 -7.254e-10 [19,] 1000 0.5134 2.786e-06 0.5134 -2.090e-06 -2.790e-07 -5.569e-10 [20,] 2000 0.5095 9.853e-07 0.5095 -7.391e-07 -9.859e-08 -9.905e-11 [21,] 3000 0.5077 5.364e-07 0.5077 -4.023e-07 -5.366e-08 -3.568e-11 [22,] 4000 0.5067 3.484e-07 0.5067 -2.613e-07 -3.485e-08 -1.527e-11 [23,] 5000 0.5060 2.493e-07 0.5060 -1.870e-07 -2.494e-08 1.862e-12 [24,] 6000 0.5055 1.897e-07 0.5055 -1.423e-07 -1.897e-08 -1.019e-11 [25,] 7000 0.5051 1.505e-07 0.5051 -1.129e-07 -1.505e-08 5.995e-13 [26,] 8000 0.5047 1.232e-07 0.5047 -9.240e-08 -1.232e-08 -9.049e-12 [27,] 9000 0.5045 1.032e-07 0.5045 -7.743e-08 -1.033e-08 3.336e-11 > matplot(k, abs(ph.mat[,-(1:2)]), type='o', log='xy') > > op <- if(require("sfsmisc")) mult.fig(2)$old.par else par(mfrow=c(2,1)) Loading required package: sfsmisc > for(Reps in c(0,1)) { + options(rErr.eps = Reps) ## 0: always RELATIVE error; 1: always ABSOLUTE + tit <- paste("phyper() approximations:", + if(Reps == 0) "relative" else "absolute", "error") + cat(tit,":\n") + x <- round(.6*k); ph.k2k <- phyper(x,1.6*k, 2*k,1.8*k) + print(ph.mat <- + cbind(k=k, phyper= ph.k2k, + rERR.Ibeta= rErr(phyperIbeta (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.as151= rErr(phyperApprAS152 (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.1mol = rErr(phyper1molenaar (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.2mol = rErr(phyper2molenaar (x,1.6*k,2*k,1.8*k) , ph.k2k), + rERR.Peiz = rErr(phyperPeizer (x,1.6*k,2*k,1.8*k) , ph.k2k), + rErr.binMol=rErr(phyperBinMolenaar.1(x,1.6*k,2*k,1.8*k) , ph.k2k)) + ) + ## The two molenaar's ``break down''; Peizer remains decent! + Err.phyp <- abs(ph.mat[,-(1:2)]) + matplot(k, Err.phyp, + ylim= if((ee <- .Options$rErr.eps) > 0) c(1e-200,1), + main = tit, type='o', log='xy', + lty=1:6, col=1:6, pch=as.character(1:6)) + if(Reps == 1) + legend("bottomleft", + c("Ibeta", "AS 152", "1_Molenaar", "2_Molenaar", "Peizer", + "binom_Mol."), + bty = "n", lty=1:6, col=1:6, pch=as.character(1:6)) + } phyper() approximations: relative error : k phyper rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peiz [1,] 10 1.573e-01 -0.023303 0.6444 2.162e-02 4.521e-04 6.654e-04 [2,] 20 4.810e-02 -0.003428 0.8714 1.667e-02 -6.849e-03 4.014e-04 [3,] 30 1.633e-02 0.020289 0.9530 9.682e-03 -1.450e-02 2.899e-04 [4,] 40 5.805e-03 0.044219 0.9827 2.331e-03 -2.220e-02 2.303e-04 [5,] 50 2.118e-03 0.067790 0.9936 -5.154e-03 -2.994e-02 1.935e-04 [6,] 60 7.863e-04 0.090861 0.9977 -1.272e-02 -3.773e-02 1.686e-04 [7,] 70 2.954e-04 0.113394 0.9991 -2.034e-02 -4.556e-02 1.505e-04 [8,] 80 1.120e-04 0.135381 0.9997 -2.802e-02 -5.344e-02 1.369e-04 [9,] 90 4.274e-05 0.156828 0.9999 -3.576e-02 -6.136e-02 1.263e-04 [10,] 100 1.640e-05 0.177745 1.0000 -4.355e-02 -6.935e-02 1.178e-04 [11,] 200 1.340e-09 0.360344 1.0000 -1.246e-01 -1.524e-01 7.908e-05 [12,] 300 1.253e-13 0.502335 1.0000 -2.118e-01 -2.417e-01 6.610e-05 [13,] 400 1.239e-17 0.612786 1.0000 -3.057e-01 -3.379e-01 5.959e-05 [14,] 500 1.264e-21 0.698716 1.0000 -4.069e-01 -4.416e-01 5.568e-05 [15,] 600 1.315e-25 0.765574 1.0000 -5.159e-01 -5.533e-01 5.308e-05 [16,] 700 1.387e-29 0.817594 1.0000 -6.334e-01 -6.736e-01 5.122e-05 [17,] 800 1.478e-33 0.858069 1.0000 -7.599e-01 -8.033e-01 4.982e-05 [18,] 900 1.588e-37 0.889563 1.0000 -8.962e-01 -9.430e-01 4.873e-05 [19,] 1000 1.716e-41 0.914068 1.0000 -1.043e+00 -1.093e+00 4.786e-05 [20,] 2000 4.451e-81 0.993008 1.0000 -3.308e+00 -3.414e+00 4.395e-05 [21,] 3000 1.332e-120 0.999431 1.0000 -8.085e+00 -8.309e+00 4.265e-05 [22,] 4000 4.225e-160 0.999954 1.0000 -1.816e+01 -1.863e+01 4.199e-05 [23,] 5000 1.384e-199 0.999996 1.0000 -3.940e+01 -4.039e+01 4.160e-05 [24,] 6000 4.628e-239 1.000000 1.0000 -8.418e+01 -8.628e+01 4.134e-05 [25,] 7000 1.569e-278 1.000000 1.0000 -1.786e+02 -1.830e+02 4.115e-05 [26,] 8000 5.376e-318 1.000000 1.0000 1.000e+00 1.000e+00 1.000e+00 [27,] 9000 0.000e+00 NaN NaN NaN NaN NaN rErr.binMol [1,] 3.958e-02 [2,] 2.892e-02 [3,] 1.952e-02 [4,] 1.056e-02 [5,] 1.767e-03 [6,] -6.954e-03 [7,] -1.566e-02 [8,] -2.437e-02 [9,] -3.312e-02 [10,] -4.190e-02 [11,] -1.331e-01 [12,] -2.318e-01 [13,] -3.389e-01 [14,] -4.552e-01 [15,] -5.815e-01 [16,] -7.189e-01 [17,] -8.681e-01 [18,] -1.030e+00 [19,] -1.206e+00 [20,] -4.072e+00 [21,] -1.066e+01 [22,] -2.579e+01 [23,] -6.058e+01 [24,] -1.405e+02 [25,] -3.243e+02 [26,] -7.467e+02 [27,] NaN phyper() approximations: absolute error : k phyper rERR.Ibeta rERR.as151 rERR.1mol rERR.2mol rERR.Peiz [1,] 10 1.573e-01 -3.665e-03 1.013e-01 3.400e-03 7.110e-05 1.046e-04 [2,] 20 4.810e-02 -1.649e-04 4.191e-02 8.017e-04 -3.294e-04 1.931e-05 [3,] 30 1.633e-02 3.313e-04 1.556e-02 1.581e-04 -2.367e-04 4.734e-06 [4,] 40 5.805e-03 2.567e-04 5.705e-03 1.353e-05 -1.289e-04 1.337e-06 [5,] 50 2.118e-03 1.436e-04 2.105e-03 -1.092e-05 -6.343e-05 4.099e-07 [6,] 60 7.863e-04 7.144e-05 7.844e-04 -9.998e-06 -2.966e-05 1.325e-07 [7,] 70 2.954e-04 3.349e-05 2.951e-04 -6.008e-06 -1.346e-05 4.447e-08 [8,] 80 1.120e-04 1.516e-05 1.119e-04 -3.137e-06 -5.983e-06 1.533e-08 [9,] 90 4.274e-05 6.702e-06 4.273e-05 -1.528e-06 -2.622e-06 5.398e-09 [10,] 100 1.640e-05 2.915e-06 1.640e-05 -7.143e-07 -1.137e-06 1.932e-09 [11,] 200 1.340e-09 4.830e-10 1.340e-09 -1.670e-10 -2.042e-10 1.060e-13 [12,] 300 1.253e-13 6.293e-14 1.253e-13 -2.653e-14 -3.028e-14 8.280e-18 [13,] 400 1.239e-17 7.591e-18 1.239e-17 -3.787e-18 -4.186e-18 7.382e-22 [14,] 500 1.264e-21 8.829e-22 1.264e-21 -5.142e-22 -5.580e-22 7.036e-26 [15,] 600 1.315e-25 1.007e-25 1.315e-25 -6.784e-26 -7.275e-26 6.979e-30 [16,] 700 1.387e-29 1.134e-29 1.387e-29 -8.786e-30 -9.345e-30 7.105e-34 [17,] 800 1.478e-33 1.269e-33 1.478e-33 -1.123e-33 -1.188e-33 7.365e-38 [18,] 900 1.588e-37 1.412e-37 1.588e-37 -1.423e-37 -1.497e-37 7.738e-42 [19,] 1000 1.716e-41 1.568e-41 1.716e-41 -1.790e-41 -1.876e-41 8.212e-46 [20,] 2000 4.451e-81 4.420e-81 4.451e-81 -1.473e-80 -1.520e-80 1.956e-85 [21,] 3000 1.332e-120 1.331e-120 1.332e-120 -1.077e-119 -1.106e-119 5.679e-125 [22,] 4000 4.225e-160 4.225e-160 4.225e-160 -7.671e-159 -7.870e-159 1.774e-164 [23,] 5000 1.384e-199 1.384e-199 1.384e-199 -5.453e-198 -5.590e-198 5.758e-204 [24,] 6000 4.628e-239 4.628e-239 4.628e-239 -3.896e-237 -3.993e-237 1.913e-243 [25,] 7000 1.569e-278 1.569e-278 1.569e-278 -2.803e-276 -2.872e-276 6.458e-283 [26,] 8000 5.376e-318 5.376e-318 5.376e-318 5.376e-318 5.376e-318 5.376e-318 [27,] 9000 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 rErr.binMol [1,] 6.225e-03 [2,] 1.391e-03 [3,] 3.188e-04 [4,] 6.128e-05 [5,] 3.743e-06 [6,] -5.467e-06 [7,] -4.625e-06 [8,] -2.729e-06 [9,] -1.415e-06 [10,] -6.873e-07 [11,] -1.784e-10 [12,] -2.904e-14 [13,] -4.198e-18 [14,] -5.752e-22 [15,] -7.647e-26 [16,] -9.972e-30 [17,] -1.283e-33 [18,] -1.636e-37 [19,] -2.070e-41 [20,] -1.812e-80 [21,] -1.419e-119 [22,] -1.090e-158 [23,] -8.385e-198 [24,] -6.504e-237 [25,] -5.089e-276 [26,] -4.014e-315 [27,] 0.000e+00 Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 6 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > par(op) > > hp1 <- list(n = 10, n1 = 12, n2 = 2) > > for(h.nr in 1:1) { + attach(hplist <- get(paste("hp",h.nr, sep='')) ) # defining n, n1, n2 + cat("\n"); print(unlist(hplist)) + N <- n1 + n2 + x <- 0:max(0,min(n1, n2 - n)) + d <- dhyper(x, n1, n2, n) + names(d) <- as.character(x) + print(d) + } The following objects are masked from hplist <- get(paste("hp", h.nr, sep = "")) (pos = 4): n, n1, n2 n n1 n2 10 12 2 0 0 > > ### ----------- ----------- lfastchoose() etc ----------------------------- > > > sapply(0:10, function(n) lfastchoose(n, c(0,n))) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 0 0 0 0 0 0 0 0 0 0 0 [2,] 0 0 0 0 0 0 0 0 0 0 0 > > ## System lchoose gives non-sense: > sapply(0:10, function(n) lchoose(n, c(0,n))) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 0 0 0 0 0 0 0 0 0 0 0 [2,] 0 0 0 0 0 0 0 0 0 0 0 > ## This is ok: > sapply(0:10, function(n) f05lchoose(n, c(0,n))) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 0 0 0 0 0 0 0 0 0 0 0 [2,] 0 0 0 0 0 0 0 0 0 0 0 > > > ###---------- some gamma testing: ------------- > > pi - gamma(1/2)^2 [1] -4.441e-16 > for(n in 1:20) cat(n,":",formatC(log(prod(1:n)) - lgamma(n+1)),"\n") 1 : 0 2 : 0 3 : 0 4 : 0 5 : 0 6 : 0 7 : 0 8 : 0 9 : 0 10 : -3.553e-15 11 : -3.553e-15 12 : 0 13 : 3.553e-15 14 : -3.553e-15 15 : -3.553e-15 16 : -3.553e-15 17 : 0 18 : 0 19 : 0 20 : 0 > for(n in 1:20) cat(n,":",formatC(prod(1:n) - gamma(n+1)),"\n") 1 : 0 2 : 0 3 : 0 4 : 0 5 : 0 6 : 0 7 : 0 8 : 0 9 : 0 10 : 0 11 : 0 12 : 0 13 : 0 14 : 0 15 : 0 16 : 0 17 : 0 18 : 0 19 : 0 20 : 0 > > ## in math/gamma.c : > p1 <- c(0.83333333333333101837e-1, + -.277777777735865004e-2, + 0.793650576493454e-3, + -.5951896861197e-3, + 0.83645878922e-3, + -.1633436431e-2) > > ### Bernoulli numbers Bern() : > > n <- 0:12 > Bn <- n; for(i in n) Bn[i+1] <- Bern(i); cbind(n, Bn) n Bn [1,] 0 1.00000 [2,] 1 0.50000 [3,] 2 0.16667 [4,] 3 0.00000 [5,] 4 -0.03333 [6,] 5 0.00000 [7,] 6 0.02381 [8,] 7 0.00000 [9,] 8 -0.03333 [10,] 9 0.00000 [11,] 10 0.07576 [12,] 11 0.00000 [13,] 12 -0.25311 > system.time(Bern(30)) user system elapsed 0 0 0 > system.time(Bern(30)) user system elapsed 0 0 0 > if(FALSE){ ## rm(.Bernoulli) + system.time(Bern(30)) + system.time(Bern(80)) + } > > ### lgammaAsymp() --- Asymptotic log gamma function : > > lgammaAsymp(3,0) [1] 0.6655 > for(n in 0:10) print(log(2) - lgammaAsymp(3,n)) [1] 0.02768 [1] -9.985e-05 [1] 3.029e-06 [1] -2.375e-07 [1] 3.468e-08 [1] -8.081e-09 [1] 2.743e-09 [1] -1.277e-09 [1] 7.819e-10 [1] -6.091e-10 [1] 5.889e-10 > for(n in 0:12) print((log(2*3*4*5) - lgammaAsymp(5+1,n)) / .Machine$double.eps) [1] 6.249e+13 [1] -5.747e+10 [1] 450434216 [1] -9221832 [1] 354336 [1] -21832 [1] 1972 [1] -236 [1] 44 [1] -4 [1] 8 [1] 4 [1] 4 > for(n in 0:12) print((log(720) - lgammaAsymp(6+1,n)) / .Machine$double.eps) [1] 5.358e+13 [1] -3.626e+10 [1] 209501224 [1] -3165248 [1] 89852 [1] -4092 [1] 276 [1] -24 [1] 4 [1] 4 [1] 4 [1] 4 [1] 4 > for(n in 0:12) print((log(720) - lgamma (6+1) ) / .Machine$double.eps) [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 [1] 0 > for(n in 0:12) print((log(5040) - lgammaAsymp(7+1,n)) / .Machine$double.eps) [1] 4.689e+13 [1] -2.433e+10 [1] 107827568 [1] -1250968 [1] 27296 [1] -952 [1] 56 [1] 0 [1] 8 [1] 8 [1] 8 [1] 8 [1] 8 > > > ## look ok: > curve(lgamma(x),-70,10, n= 1001) # ok > curve(lgamma(x),-70,10, n=20001) # slowness of x11 !! > > curve(abs(gamma(x)),-70,10, n= 201, log = 'y') Warning message: In gamma(x) : NaNs produced > curve(abs(gamma(x)),-70,70, n=10001, log = 'y') Warning message: In gamma(x) : NaNs produced > par(new=T) > curve(lgamma(x), -70,70, n= 1001, col = 'red') Warning messages: 1: In eval(expr, envir = ll, enclos = parent.frame()) : full precision may not have been achieved in 'lgamma' 2: In eval(expr, envir = ll, enclos = parent.frame()) : full precision may not have been achieved in 'lgamma' 3: In eval(expr, envir = ll, enclos = parent.frame()) : full precision may not have been achieved in 'lgamma' > > ## .., we should NOT use log - scale with negative values... [graph ok, now] > curve(gamma(x),- 40,10, n= 2001, log = 'y') Warning messages: 1: In gamma(x) : NaNs produced 2: In xy.coords(x, y, xlabel, ylabel, log) : 780 y values <= 0 omitted from logarithmic plot > curve(abs(gamma(x)),-40,10, n= 2001, log = 'y') Warning message: In gamma(x) : NaNs produced > ##- Warning: NAs produced in function "gamma" > > x <- seq(-40,10, length = 201) > plot(x, gamma(x), log = 'y', type='h') Warning messages: 1: In gamma(x) : NaNs produced 2: In xy.coords(x, y, xlabel, ylabel, log) : 60 y values <= 0 omitted from logarithmic plot > > > ###------------------- early dhyper() problem --------------------------------- > > ##- Date: 30 Apr 97 11:05:00 EST > ##- From: "Ennapadam Venkatraman" > ##- Subject: R-beta: dhyper bug?? > ##- To: "r-testers" > ##- Sender: owner-r-help@stat.math.ethz.ch > > ## I get the following incorrect answer ---- FIXED > > dhyper(34,410,312,49) [1] 0.02189 > # [1] 2.244973e-118 > ##- when the coorect value is > ##- > dhyper(34,410,312,49)# (from S-Plus) [1] 0.02189 > ##- [1] 0.0218911 > > ##================== R is now correct ! > stopifnot(all.equal(dhyper(34,410,312,49), + 0.021891095726, tol=1e-11)) > ##E.S. Venkatraman (venkat@biosta.mskcc.org) > > > ### phyperR() === R version of pre-2004 C version in /src/nmath/phyper.c : > > ## This takes long, currently (1999??) {18 sec, on sophie [Ultra 1]} -- now all fast! > k <- (1:100)*1000 > system.time(phyper.k <- phyper(k, 2*k, 2*k, 2*k)) user system elapsed 0 0 0 > > k <- 1000; phyperR(k, 2*k, 2*k, 2*k) - phyper.k[1] [1] -1.321e-13 > > ## Debug: ------------ > if(FALSE) # for now, when I source this + if(interactive()) + debug(phyperR) > phyperR(k, 2*k, 2*k, 2*k) [1] 0.5126 > ## before while(): > xb <- 2000 > xr <- 0 > ltrm <- -2768.21584356709 > NR <- 2000 > NB <- 0# new value > > > proc.time() user system elapsed 35.25 3.26 38.53