[Rd] rgamma gives zeros. (PR#9184)

[ripley at stats.ox.ac.uk]

I think you fail to understand the floating point arithmetic involved. The 
smallest floating point number that can be represented to full precision 
is

> .Machine$double.xmin
[1] 2.225074e-308

and

> pgamma(.Machine$double.xmin, 0.001, rate=0.01)
[1] 0.4904533

Your distribution is so extreme that in computer representation it is 
largely discrete.

Now, on your platform there are denormalized doubles, so

> table(rgamma(1e6, 0.001, rate=0.01) < .Machine$double.xmin)

 FALSE   TRUE
509073 490927
> table(rgamma(1e6, 0.001, rate=0.01) == 0)

 FALSE   TRUE
525187 474813

> x <- rgamma(1e6, 0.001, rate=0.01)
> sort(x[x> 0]) [1:10]
 [1] 4.940656e-322 4.940656e-322 4.940656e-322 4.940656e-322 4.940656e-322
 [6] 4.940656e-322 4.940656e-322 4.940656e-322 4.940656e-322 4.940656e-322

is pretty much what I expected.

On Tue, 29 Aug 2006, pxi at stat.cmu.edu wrote:

> Full_Name: Peiyi Xi
> Version: R 2.2.0

You are specifically asked NOT to report on obselete versions of R.

> OS: Windows XP Professional
> Submission from: (NULL) (128.2.3.141)
> 
> 
> When I use rgamma(n, shape, rate ) to generate gamma samples, it
> gives zeros when both shape and rate are very small.
> 
> But we know that if x follows a gamma distribution, x should be positive.
> 
> e.g.
> > temp=rgamma(10, 0.001, rate=0.01)
> > temp
>  [1]  2.438078e-33 5.101136e-130  1.760830e-54 2.724731e-166  0.000000e+00
>  [6] 4.461151e-146  1.332914e-55 2.336396e-277  0.000000e+00  0.000000e+00
> > temp[5]
> [1] 0
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
> 

-- 
Brian D. Ripley,                  ripley at 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