Bug 12324 - qgamma inaccuracy
qgamma inaccuracy
Status: CLOSED FIXED
Product: R
Classification: Unclassified
Component: Accuracy
old
All All
: P5 normal
Assigned To: Jitterbug compatibility account
Depends on:
Blocks:
  Show dependency treegraph
 
Reported: 2008-08-06 01:59 UTC by Jitterbug compatibility account
Modified: 2010-07-14 11:07 UTC (History)
1 user (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Jitterbug compatibility account 2008-08-06 01:59:48 UTC
From: skylab.gupta@gmail.com
Full_Name: 
Version: 2.7.1 (2008-06-23)
OS: windows vista
Submission from: (NULL) (216.82.144.137)


Hello,

I have been working with various probability distributions in R, and it seems
the gamma distribution is inaccurate for some inputs.

For example, qgamma(1e-100, 5e-101, lower.tail=FALSE) gives: 1.0. However, it
seems this is incorrect; I think the correct answer should be
0.082372029620717283. When I check these numbers using pgamma, I get:

pgamma(1,5e-101, lower.tail=FALSE) = 9.1969860292859463e-102
and
pgamma(0.082372029620717283,5e-101, lower.tail=FALSE) =
1.0000000000000166e-100.

Similarly, for example:
qgamma(1e-100,0.005,lower.tail=FALSE) = 109.36757177007101
pgamma(109.36757177007101, 0.005, lower.tail=FALSE) = 1.4787306506694758e-52.

This looks completely wrong. The correct value, I think, should be
219.59373661415756. In fact,
pgamma(219.59373661415756, 0.005, lower.tail=FALSE) = 9.9999999999999558e-101.

In fact, when I do the following in R, the results are completely wrong,

x<-c(5e-1,5e-2,5e-3,5e-4,5e-5,5e-6,5e-7,5e-8,5e-9,5e-10)
z1 <-qgamma(1e-100,x,lower.tail=FALSE)
y<-pgamma(z1,x,lower.tail=FALSE)

The value of y that I get should be close to 1e-100, but they are not:

> y
 [1] 1.000000e-100  1.871683e-51  1.478731e-52  1.444034e-53  1.440606e-54
 [6]  1.440264e-55  1.440230e-56  1.440226e-57  1.440226e-58  1.440226e-59

The correct values of z1 should be:
z1true <- c(226.97154111939946, 222.15218724493326, 219.59373661415756,
217.27485383840451, 214.98015408183574, 212.68797118872064, 210.39614286838227,
208.10445550564617, 205.81289009100664, 203.52144711679352)

With these values of z1true, we get the expected values:
y<-pgamma(z1true,x,lower.tail=FALSE)
> y
 [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100


I am using the precompiled binary version of R, under Windows Vista.
-----------
> version
    _                           
platform       i386-pc-mingw32             
arch           i386                        
os             mingw32                     
system         i386, mingw32               
status                                     
major          2                           
minor          7.1                         
year           2008                        
month          06                          
day            23                          
svn rev        45970                       
language       R                           
version.string R version 2.7.1 (2008-06-23)
------------

So, it seems qgamma is inaccurate for small probability values in the upper
tail, when the shape parameter is also small.

Comment 1 Jitterbug compatibility account 2008-08-13 00:50:50 UTC
From: Prof Brian Ripley <ripley@stats.ox.ac.uk>
This is a really extreme usage.  AFAICS the code works well enough down to 
shape=1e-10 or so, e.g.

> qgamma(1e-10, 5e-11, lower.tail=FALSE)
[1] 0.08237203

I would be interested to know what substantive problem you were trying to 
solve here that required such values.

I am pretty sure that a completely different algorithm will be required. 
For completeness we may write that in due course, but for now (R 2.7.2) I 
suggest just issuing a warning for miniscule 'shape'.

On Thu, 7 Aug 2008, skylab.gupta@gmail.com wrote:

> Full_Name:
> Version: 2.7.1 (2008-06-23)
> OS: windows vista
> Submission from: (NULL) (216.82.144.137)
>
>
> Hello,
>
> I have been working with various probability distributions in R, and it seems
> the gamma distribution is inaccurate for some inputs.
>
> For example, qgamma(1e-100, 5e-101, lower.tail=FALSE) gives: 1.0. However, it
> seems this is incorrect; I think the correct answer should be
> 0.082372029620717283. When I check these numbers using pgamma, I get:
>
> pgamma(1,5e-101, lower.tail=FALSE) = 9.1969860292859463e-102
> and
> pgamma(0.082372029620717283,5e-101, lower.tail=FALSE) =
> 1.0000000000000166e-100.
>
> Similarly, for example:
> qgamma(1e-100,0.005,lower.tail=FALSE) = 109.36757177007101
> pgamma(109.36757177007101, 0.005, lower.tail=FALSE) = 1.4787306506694758e-52.
>
> This looks completely wrong. The correct value, I think, should be
> 219.59373661415756. In fact,
> pgamma(219.59373661415756, 0.005, lower.tail=FALSE) = 9.9999999999999558e-101.
>
> In fact, when I do the following in R, the results are completely wrong,
>
> x<-c(5e-1,5e-2,5e-3,5e-4,5e-5,5e-6,5e-7,5e-8,5e-9,5e-10)
> z1 <-qgamma(1e-100,x,lower.tail=FALSE)
> y<-pgamma(z1,x,lower.tail=FALSE)
>
> The value of y that I get should be close to 1e-100, but they are not:
>
>> y
> [1] 1.000000e-100  1.871683e-51  1.478731e-52  1.444034e-53  1.440606e-54
> [6]  1.440264e-55  1.440230e-56  1.440226e-57  1.440226e-58  1.440226e-59
>
> The correct values of z1 should be:
> z1true <- c(226.97154111939946, 222.15218724493326, 219.59373661415756,
> 217.27485383840451, 214.98015408183574, 212.68797118872064, 210.39614286838227,
> 208.10445550564617, 205.81289009100664, 203.52144711679352)
>
> With these values of z1true, we get the expected values:
> y<-pgamma(z1true,x,lower.tail=FALSE)
>> y
> [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
>
>
> I am using the precompiled binary version of R, under Windows Vista.
> -----------
>> version
>    _
> platform       i386-pc-mingw32
> arch           i386
> os             mingw32
> system         i386, mingw32
> status
> major          2
> minor          7.1
> year           2008
> month          06
> day            23
> svn rev        45970
> language       R
> version.string R version 2.7.1 (2008-06-23)
> ------------
>
> So, it seems qgamma is inaccurate for small probability values in the upper
> tail, when the shape parameter is also small.
>
> ______________________________________________
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

Comment 2 Jitterbug compatibility account 2008-12-12 22:17:00 UTC
NOTES:
 Needs new algorithm for really extreme case: warns for now
Comment 3 Jitterbug compatibility account 2008-12-12 23:20:27 UTC
Audit (from Jitterbug):
Sat Aug  9 08:06:45 2008	ripley	moved from incoming to Accuracy
Fri Dec 12 17:20:27 2008	ripley	changed notes
Comment 4 Martin Maechler 2010-07-14 11:07:14 UTC
These examples have been fixed by

------------------------------------------------------------------------
r50393 | maechler | 2009-11-11 14:58:22 +0100 (Wed, 11 Nov 2009) | 1 line

pgamma() and qgamma() upper-tail improvements, partly related to PR#12324
------------------------------------------------------------------------

hence in R versions 2.11.0 and up.