ARS1D: Derivative-free Adaptive Rejection Sampler for 1-D functions AUTHOR Tom Stepleton VERSION This is version 0.5, the initial release, dated 26 April 2008. See the CHANGES section below for more details on this and prior releases. ABOUT This is an implementation of the derivative-free Adaptive Rejection Sampler algorithm. It generates samples efficiently from 1-D log-concave density functions. "Derivative-free" means that it's not necessary to supply the derivative of the density function to the sampler. For more information about the Adaptive Rejection Sampler algorithm, refer to these slides: http://www.stat.duke.edu/~cnk/Links/slides.pdf A more formal reference is W.R. Gilks (1992), "Derivative-free Adaptive Rejection Sampling for Gibbs Sampling," /Bayesian Statistics 4/, (eds. J. Bernardo, J. Berger, A.P. Dawid, and A.F.M. Smith.) Oxford University Press, 641-649. The sampler can sample functions with bounded support (e.g. the Beta distribution PDF), semi-infinite support (e.g. the Gamma distribution PDF), or infinite support (e.g. normal distributions). The sampler is written in templated C++. A MATLAB MEX interface is supplied. The code is well commented and features Doxygen annotation. The sampler takes pains to avoid numerical stability issues. Internal representations are dynamically rescaled and otherwise adjusted to prevent overflow conditions. If necessary, the code using the sampler need only provide the logarithms of density function values instead of the values themselves. Density functions need not sum to 1. The sampler has been compiled and run on a MacOS X 10.5 computer with a 2.33 GHz Intel processor. Each round of sampling in the MATLAB test script test_ars1d.m generates 10,000,000 samples of the furnished density function in around 3 seconds, which includes setup and teardown of supporting structures for the MATLAB interface and a few hundred callbacks to MATLAB to evaluate the (interpreted) density function. The sampler is new code and has not been thoroughly vetted for accuracy or stability. While the author has taken pains to ensure its correctness, you use the code at your own risk! BUILDING THE MATLAB MEX FUNCTION The sampler uses the Boost suite of C++ libraries for random number generation. Boost is available for free at http://www.boost.org/. You will need the Boost library With Boost, hopefully all it should take to compile this code is editing the Makefile to correctly identify the paths to the Boost library and the mex executable. Documentation for the MATLAB interface to the sampler exists in the "help text" header of the ars1d.m function file. USE IN C++ PROJECTS The sampler is entirely contained within ars1d.h. As templated code, there is no library to build; simply include the header file and furnish a Boost library random number generator. The MEX function source code mex_ars1d.cc should serve as a useful example of how to incorporate the sampler into your project. Doxygen documentation for the ARS1D::Sampler class and affiliated classes has been included in this source code distribution (see doc/html/index.html). CONCERNS ABOUT THE BOOST LIBRARY uniform_real CLASS As of 26 April 2008, the Boost library uniform_real class documentation includes the following disclaimer: Note: The current implementation is buggy, because it may not fill all of the mantissa with random bits. The uniform_real class transforms a Boost random number generator's native output (usually uniformly-distributed unsigned integers) into uniformly- distributed real numbers within a specified interval. Its approach for doing so is common to many programs. Assuming without loss of generality that the desired interval is [0,a), and that the sampler samples natural numbers in {0, 1, ..., n-1} uniformly, pseudocode for the present Boost uniform_real method is: integer r = sample_a_natural_number() real s = a * (integer_to_real(r) / integer_to_real(n)) return s Whether this method is acceptable is up to you, the user. The sampler will sample numbers in the intervals [0,1) and [0,c+s), where c is the integral of your density function over its support and s is a positive number. For the first few (< 25) samples s may be much larger than c, but subsequently s diminishes rapidly toward 0. This occurs as the piecewise exponential upper bound approximation to the density function is refined; indeed, s is the difference between the upper bound's integral and the density function's integral. CHANGES Version Notes ------- ----- 0.5 Initial release