Bugzilla – Bug 12324

qgamma inaccuracy

Last modified: 2010-07-14 11:07:14 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.

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

NOTES: Needs new algorithm for really extreme case: warns for now

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

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.