>>>>> "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).