[Rd] Buglet in qbeta?

[Martin Maechler]

>>>>> "JL" == Josef Leydold <leydold at statmath.wu.ac.at>
>>>>>     on Wed, 7 Oct 2009 14:48:26 +0200 writes:

    JL> Hi,
    JL> I sometimes play around with extreme parameters for distributions and
    JL> found that qbeta is not always monotone as the following example shows.
    JL> I don't know whether this is serious enough to submit a bug report (as
    JL> this example is near to the limitations of floating point arithmetic).

well, it's not near enough the limits to *not* warrant a bug
report!

I "know" (I have saved a collection of ..) about several
accuracy problems of pbeta() / qbeta(), most cases indeed
"borderline" (at the limits of FP arithmetic) but this one may
well be a new one.

A bit more succint than the dump of your numbers is a simple
graphic

  p <- (1:100)/100; qp <- qbeta(p, 0.01,5)
  plot(p,qp, type="o", log = "y")

  ## or, already suggesting a fix to the bug:

  plot(p,qp, type="o", log = "xy")

which is definitely a bug...  though not a terrible one.

@ R-core: please leave me a bit of time before touching qbeta.c,
  	  as I already have a modified local version of that.

Thanks for the report, Josef!
--
Martin Maechler, ETH Zurich

    >> x <- qbeta((0:100)/100,0.01,5)
    >> x
    JL> [1]  0.000000e+00 1.253990e-201 1.589622e-171 6.462785e-154 2.015085e-141
    JL> [6] 9.892240e-132 8.192553e-124 4.056003e-117 2.554424e-111 3.330774e-106
    JL> [11] 1.253990e-101  1.728076e-97  1.038529e-93  3.109063e-90 5.141594e-87
    JL> [16]  5.098238e-84  3.238117e-81  1.390549e-78  4.222258e-76 9.411402e-74
    JL> [21]  1.589622e-71  2.090373e-69  2.190596e-67  1.866714e-65 1.316493e-63
    JL> [26]  7.803602e-62  3.941205e-60  1.716606e-58  6.517745e-57 2.178181e-55
    JL> [31]  6.462785e-54  1.715788e-52  4.104801e-51  8.906113e-50 1.762731e-48
    JL> [36]  3.199622e-47  5.352348e-46  8.288322e-45  1.193037e-43 1.602341e-42
    JL> [41]  2.015085e-41  2.380564e-40  2.649862e-39  2.787018e-38 2.776910e-37
    JL> [46]  2.627517e-36  2.366341e-35  2.032732e-34  1.668853e-33 1.311905e-32
    JL> [51]  9.892240e-32  2.220446e-15  2.331468e-15  2.553513e-15 1.110223e-16
    JL> [56]  3.330669e-16  9.992007e-16  8.881784e-16  4.440892e-16 7.771561e-16
    JL> [61]  1.554312e-15  8.881784e-16  9.992007e-16  2.109424e-15 2.331468e-15
    JL> [66]  2.553513e-15  2.553513e-15  1.110223e-16  6.661338e-16 1.221245e-15
    JL> [71]  1.443290e-15  2.220446e-16  6.661338e-16  2.664535e-15 1.054712e-14
    JL> [76]  4.019007e-14  1.512124e-13  5.589973e-13  2.031153e-12 7.261081e-12
    JL> [81]  2.554423e-11  8.847001e-11  3.017724e-10  1.014153e-09 3.359099e-09
    JL> [86]  1.096950e-08  3.532966e-08  1.122586e-07  3.520157e-07 1.089674e-06
    JL> [91]  3.330818e-06  1.005659e-05  3.000050e-05  8.845884e-05 2.579427e-04
    JL> [96]  7.446202e-04  2.133444e-03  6.108393e-03  1.783085e-02 5.699554e-02
    JL> [101]  1.000000e+00
    >> order(x)
    JL> [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
    JL> [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
    JL> [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  55  68  72
    JL> [55]  56  59  69  73  60  58  62  57  63  70  71  61  64  52  53  65  54  66
    JL> [73]  67  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
    JL> [91]  91  92  93  94  95  96  97  98  99 100 101
    >> pbeta(x,0.01,5)
    JL> [1] 0.0000000 0.0100000 0.0200000 0.0300000 0.0400000 0.0500000 0.0600000
    JL> [8] 0.0700000 0.0800000 0.0900000 0.1000000 0.1100000 0.1200000 0.1300000
    JL> [15] 0.1400000 0.1500000 0.1600000 0.1700000 0.1800000 0.1900000 0.2000000
    JL> [22] 0.2100000 0.2200000 0.2300000 0.2400000 0.2500000 0.2600000 0.2700000
    JL> [29] 0.2800000 0.2900000 0.3000000 0.3100000 0.3200000 0.3300000 0.3400000
    JL> [36] 0.3500000 0.3600000 0.3700000 0.3800000 0.3900000 0.4000000 0.4100000
    JL> [43] 0.4200000 0.4300000 0.4400000 0.4500000 0.4600000 0.4700000 0.4800000
    JL> [50] 0.4900000 0.5000000 0.7285871 0.7289426 0.7296061 0.7070842 0.7148952
    JL> [57] 0.7227924 0.7219416 0.7169548 0.7209782 0.7259930 0.7219416 0.7227924
    JL> [64] 0.7282134 0.7289426 0.7296061 0.7296061 0.7070842 0.7198677 0.7242443
    JL> [71] 0.7254552 0.7120024 0.7198677 0.7299167 0.7400284 0.7499948 0.7599988
    JL> [78] 0.7700008 0.7799998 0.7900000 0.8000000 0.8100000 0.8200000 0.8300000
    JL> [85] 0.8400000 0.8500000 0.8600000 0.8700000 0.8800000 0.8900000 0.9000000
    JL> [92] 0.9100000 0.9200000 0.9300000 0.9400000 0.9500000 0.9600000 0.9700000
    JL> [99] 0.9800000 0.9900000 1.0000000
    >> version
    JL> _                                                                
    JL> platform       x86_64-unknown-linux-gnu                                         
    JL> arch           x86_64                                                           
    JL> os             linux-gnu                                                        
    JL> system         x86_64, linux-gnu                                                
    JL> status         Under development (unstable)                                     
    JL> major          2                                                                
    JL> minor          11.0                                                             
    JL> year           2009                                                             
    JL> month          10                                                               
    JL> day            07                                                               
    JL> svn rev        49963                                                            
    JL> language       R                                                                
    JL> version.string R version 2.11.0 Under development (unstable) (2009-10-07
    JL> r49963)

    JL> p.s. there are similar results for R-2.9.2 in Windows (with
    JL> different round-off errors).