[Rd] Inaccurate qgamma() (PR#11030)

[Duncan Murdoch]

On 24/03/2008 1:35 PM, murdoch at stats.uwo.ca wrote:
> I haven't looked inside to see what is causing this, but there's a big 
> discontinuity in qgamma:
> 
>   curve(qgamma(x, shape=19), from=1e-10, to=2e-10)
> 
> This appears in both R-patched and R-devel.

With debugging turned on, the inaccurate value prints this:

 > qgamma(1.2e-10, shape=19)
qgamma(p=1.2e-10, alpha=     19, scale=      1, l.t.= 1, log_p= 0):  nu 
 > .32: Wilson-Hilferty; x = -6.33328
         ==> ch =    5.03581: Ph.II iter; ch=5.03581, p2=8.8438e-11
      it=2,  ch=-0.635427, p2=4.94066e-324
      it=3,  ch=1.2e-10, p2=4.94066e-324
      it=4,  ch=nan, p2=4.94066e-324

  it=1: p=1.2e-10, x = 2.51791, p.=3.1562e-11; p1:=D{p}=-8.8438e-11
no Newton step done since delta{p} >= last delta
[1] 2.517907

Things are fine if I use log.p=TRUE:

 > qgamma(log(1.2e-10), shape=19, log.p=TRUE)
qgamma(p=-22.8435, alpha=     19, scale=      1, l.t.= 1, log_p= 1):  nu 
 > .32: Wilson-Hilferty; x = -6.33328
         ==> ch =    5.03581: Ph.II iter; ch=5.03581, p2=8.8438e-11
      it=2,  ch=-0.635427, p2=4.94066e-324
      it=3,  ch=1.2e-10, p2=4.94066e-324
      it=4,  ch=nan, p2=4.94066e-324

  it=1: p=-22.8435, x = 2.51791, p.=-24.1791; p1:=D{p}=-1.33554
          it=2,  d{p}=-0.0581595
          it=3,  d{p}=-0.00011856
          it=4,  d{p}=-4.9439e-10
          it=5,  d{p}=1.42247e-15
[1] 2.729837

Maybe we should switch to this scale when the first try fails?

Duncan Murdoch