Duncan Murdoch wrote:
On 11/12/2009 7:12 AM, khaz...@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
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
