[R] how can generate from trunceted gamma distribution in R ?

[Peter Ruckdeschel]

> Duncan Murdoch wrote:
>> On 11/12/2009 7:12 AM, khazaei at ceremade.dauphine.fr wrote:
>>> Hi, all
>>> How can generate a sample from truncated inverse  gamma distribution
>>> in R?
>>
>> Using the inverse CDF method or rejection sampling are possible,
>> depending on what your truncation is like.  If your truncation forces
>> the observations far out into the tails, you need to be careful about
>> rounding and underflow when using the the inverse CDF method.
>>
>> Duncan Murdoch
>>
> 
> I think perusal of this paper might be a good idea:
> 
> Sampling Truncated Normal, Beta, and Gamma Densities
> Paul Damien and Stephen G. Walker
> Journal of Computational and Graphical Statistics, Vol. 10, No. 2 (Jun.,
> 2001), pp. 206-215
> 
> Remembering that the inverse gamma is the inverse of a gamma, you may be
> able to get a truncated inverse gamma from a truncated (at the other
> end) gamma. Alternatively, the methodology outlined in the paper most
> likely can be modified for the inverse gamma.
> 
> David Scott
> 

While David Scott's reply definitely gives a more problem specific
solution than we could offer, you might also want to look at package
"distr" on CRAN where a general truncation operator for distributions is
provided --- see ?Truncate (after installing/attaching package distr).

The inverse Gamma so far is not implemented to distr as an S4-class
(you could easily do this yourself, though). But, as David Scott
mentioned you can produce it by something along the lines

require(distr)
G0 <- Gammad(scale = 2.3, shape = 1.4) ## generates a Gamma distribution
G <- 1/G0 ## the corresponding inverse Gamma
d(G)(2) ### density of G at 2
p(G)(4) ## cdf of G at 4 ...

## example for Truncated G
TG <- Truncate(G, lower=0, upper=0.9)
 ## the lower=0  is somehow redundant in this case, will see if this
 ## can be set automatically in a next release..

q(TG)(0.99) ## upper 1% quantile

## and some functionals:
require(distrEx)
E(TG)
mad(TG)
sd(TG)

## I am not claiming that this code gives extremely accurate results,
## but for higher accuracy, you could easily overload operator "/" for
## operands "numeric", "Gammad" (and if you like, correspondingly,
## write methods for E())

Note that our package even takes into account that you might
want to use log-scales if you are interested in the tails,
so it takes up Duncan Murdoch's comment in some sense,
automatically.


Best, Peter