  The timing of radio pulsars in binary systems provides a superb testing
ground of general relativity. Here we propose a Bayesian approach to carry out
these tests, and a relevant efficient numerical implementation, that has
several conceptual and practical advantages with respect to traditional methods
based on least-square-fits that have been used so far: (i) it accounts for the
actual structure of the likelihood function - and it is not predicated on the
Laplace approximation which is implicitly built in least-square fits that can
potentially bias the inference - (ii) it provides the ratio of the evidences of
any two models under consideration as the statistical quantity to compare
different theories, and (iii) it allows us to put joint constraints from the
monitoring of multiple systems, that can be expressed in terms of ratio of
evidences or probability intervals of global (thus not system-dependent)
parameters of the theory, if any exists. Our proposed approach optimally
exploits the progress in timing of radio pulsars and the increase in the number
of observed systems. We demonstrate the power of this framework using simulated
data sets that are representative of current observations.
