Longhai Li, Department of Mathematics and Statistics, University of Saskatchewan
Adaptive
Rejection Sampling (ARS) is an efficient method for sampling from a
target distribution whose log PDF is concave. The description of
ARS by the original inventor (Wally Gilks) is as follows:
ARS works by constructing an envelope function of the log of the target density, which is then used in rejection sampling. Whenever a point is rejected by ARS, the envelope is updated to correspond more closely to the true log density, thereby reducing the chance of rejecting subsequent points. Fewer ARS rejection steps implies fewer point-evaluations of the log density. Such evaluations are typically very expensive computationally in applications of Gibbs sampling. |
This
method is very useful and conceptually simple. The original authors
provide a FORTRAN routine. However, it seemed to me that there wasn't a
C routine available. Therefore I posted this C function, which can be
called directly in R with .C function.
version of 2010-12-02: ars.c
To use the main function sample_ars, you need only to include the file ars.c in your C code. See the following examples for demonstration.
Truncated Normal Distribution
C sampling function, R wrapper function and testing code, Testing plots
Posterior of Hyperparameter of Inverse-Gamma (Inverse-Chisquare) Distribution
C sampling function, R wrapper function and testing code, Testing plots
Truncated Beta
C sampling function, R wrapper function and testing code, Testing plots
The target distribution can be badly unnormalized as long as its log PDF values are not underflow. This is convenient when the target distribution is a univariate conditional distribution, whose value is known only to be some multiple of the value of the joint distribution of all parameters, which is typical in Gibbs sampling. This is realized by computing only the log of target PDF and its derivative, and also only computing the log of the integral under each piece of exponential linear function.
The bounds of target distribution can be defined implicitly in the function eval_logf, not necessarily specified in lb and ub. This is realized by stepping out from the initial tangent point to both directions until reaching a point with log PDF equal to -INFINITY (ie reaching a point outside the bounds). The bounds will be further refined as new point is drawn. This is convenient when the bounds cannot be found in closed form but a point can be checked with a condition whether still inside bounds, for example an interval along one axis within an ellips.
Only one initial tangent point is required to start the program. More tangent points will be found by stepping out from the first point until an upper hull with finite integral is found.
