R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> library(DPQ)
> 
> ### Some "manual" tests for ../R/pnbeta-approx2.R
> ###                         ~~~~~~~~~~~~~~~~~~~~~
> 
> ## Some values from  Table 2 and 3 (p.152) of the paper:
> ## T.1
> pnbetaAppr2(0.868,  10, 20, ncp=150)# 0.937569774351
[1] 0.9375698
> ## T.2
> pnbetaAppr2(0.644,  10,  5, ncp= 20)# 0.0497384
[1] 0.04973844
> pnbetaAppr2(0.80,   30, 30, ncp=140)# 0.7816822
[1] 0.7816822
> 
> ## more extreme cases
> x <- c(0.9, 0.95, 0.96, 0.97, 0.99)
> pnbetaAppr2(x, 10, 20, 1000)# 7.953858e-08 8.276896e-02 3.709090e-01 8.216320e-01 9.999989e-01
[1] 7.953858e-08 8.276896e-02 3.709090e-01 8.216320e-01 9.999989e-01
> pbeta       (x, 10, 20, 1000)# 6.959528e-08 8.289800e-02 3.709609e-01 8.214000e-01 9.999992e-01
[1] 6.959528e-08 8.289800e-02 3.709609e-01 8.214000e-01 9.999992e-01
> ## in the tail, they differ very considerably; in the center (P in (1e-7, 1-1e-7), things are ok:
> pnbetaAppr2(x, 10, 30, 5000)
[1] 2.917795e-67 6.793172e-25 4.510392e-17 7.178493e-10 7.961821e-01
> pbeta       (x, 10, 30, 5000)
[1] 1.358084e-75 1.603009e-33 1.816631e-26 2.185984e-10 7.960078e-01
> 
> ## really large lambda: now "traditional" pbeta() seems hopeless
> x <- 1 - 2^(-5:-15)
> px <- cbind(appr2 = pnbetaAppr2(x, 10, 30, 1e6, log.p=TRUE),
+             R     = pbeta      (x, 10, 30, 1e6, log.p=TRUE))
> cbind(x, px)
              x         appr2             R
 [1,] 0.9687500 -6.824379e+03          -Inf
 [2,] 0.9843750 -3.956463e+03          -Inf
 [3,] 0.9921875 -2.251496e+03          -Inf
 [4,] 0.9960938 -1.243353e+03          -Inf
 [5,] 0.9980469 -6.563719e+02          -Inf
 [6,] 0.9990234 -3.235695e+02 -4.070011e+02
 [7,] 0.9995117 -1.428812e+02 -1.851310e+02
 [8,] 0.9997559 -5.185496e+01 -8.417972e+01
 [9,] 0.9998779 -1.236178e+01 -1.244352e+01
[10,] 0.9999390 -8.248203e-01 -8.247854e-01
[11,] 0.9999695 -5.689747e-04 -5.485346e-04
> # either '-Inf'; or full warnings about precision
> ## now (TOMS algo): much better and basically correct (but still -Inf !!)
> matplot(x, px, type ="b")
> 
> proc.time()
   user  system elapsed 
   0.18    0.06    0.23