[R] many zeroes in rgamma ... what's going on?

Prof Brian Ripley ripley at stats.ox.ac.uk
Mon Nov 19 09:02:30 CET 2007


On Sun, 18 Nov 2007, Gregory Gentlemen wrote:

> Dear Dr. Dalgaard,
>
> Thank you for your insight! In fact, I did read the example 
> documentation, however, it pretty much told me the same thing that my 
> little simulation did: there is ALOT of point mass at zero.
>
> Is there any fix to this problem? Seeing that rgamma won't work 
> accurately, if I wanted to plot a density of an inverse gamma 
> distribution with small scale and shape parameters, how would I do so?

You haven't understood the issue snown in the example.  This is not about
'won't work accurately', but 'can't work accurately':  half the mass is on 
numbers which cannot be represented in your computer.

and Vincent Goulet wrote

> Package actuar has the {d,p,q,r}invgamma() functions (and quite a few
> others), if this can be of any help to you.

But it cannot, because the reciprocals cannot be represented either (and 
actually the issue is a little worse because there is no gradual 
overflow).  And indeed that is what happens if you try the suggestion.


> Peter Dalgaard <p.dalgaard at biostat.ku.dk> wrote: Gregory Gentlemen wrote:
>> Hello fellow R users,
>>
>> I wanted to view the density on the standard deviation scale of a 
>> gamma(0.001, 0.001) prior for the precision. I did this as seen in the 
>> code below and found that for some reason rgamma is giving many values 
>> equal to zero, which is strange since a gamma distribution is 
>> continuous. What is going on here?
>>
>> Thanks for any help in advance.
>> Greg
>>
> That sort of shape parameter gives a distribution with most of its mass
> squashed against the y axis, so random numbers underflow to zero. But
> why did you not read the Example section of help(rgamma)? The effect is
> clearly indicated there.
>
>>
>>> x1 <- rgamma(10000, shape=0.001, scale=0.001)
>>> sd1 <- 1/sqrt(x1)
>>> truehist(sd1, xlim=c(0, 1.5))
>>>
>> Error in truehist(sd1, xlim = c(0, 1.5)) :
>>         'nbins' must result in a positive integer
>>
>>> summary(sd1)
>>>
>>       Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
>>  2.266e+01  9.311e+66 3.250e+153        Inf        Inf        Inf


-- 
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



More information about the R-help mailing list