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