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