No need for integrate(), uniroot(), etc Exploit the structure of the distribution using something like quantile <- scale * ( 1/qbeta(p, shape1=p, shape2=q) - 1 )^(-1/a) where the parameterization is as in C. Kleiber and S. Kotz: Statistical Size Distributions in Economics and Actuarial Sciences, Wiley 2003 (you seem to know the book) or contact me off-list, I have code for this. Best, CK
____________________________________________ Dr. Christian Kleiber FB Statistik Universitaet Dortmund Vogelpothsweg 78 44221 Dortmund Germany phone. ++49-(0)231-755 5419 fax. ++49-(0)231-755 5284 e-mail: kleiber at statistik.uni-dortmund.de