In Bayesian statistical inference and computationally intensive frequentist inference, one is interested in obtaining samples from densities which may be many dimensional, multi-modal, scale-sensitive, and unnormalized. One of the simplest Monte Carlo methods is the rejection sampler due to von Neumann. Here we introduce an auto-validating version of the rejection sampler via interval analysis. In particular, we adaptively partition the domain into boxes, use interval analysis to obtain an upper bound on the density's shape in each box, and thus construct a rigorous envelope function and a corresponding density which is easily sampled from.

We show that such a sampler does provide us with independent samples from a large class of target densities (Elementary Lipschitz class, ie densities that can be represented as a finite composition of standard arithmetic operations, standard functions and simple functions with known singularities) in a guaranteed manner. The guarantees here are synonyms for computer-assisted proofs.

We illustrate the efficiency of the sampler by theory and by examples in up to 10 dimensions. Our sampler is immune to the `pathologies' of some infamous densities, including the multi-variate versions of highly spiky witch's hat and multi-modal Rosenbrock, that primarily arise from the assumption that real numbers are representable and manipulable in computers with finite memory. We also draw provably i.i.d. samples from the posterior distribution over piece-wise Euclidean spaces of small primate phylogenies and thereby solve an open problem in statistical molecular evolution.