  Pulsar-timing datasets have been analyzed with great success using
probabilistic treatments based on Gaussian distributions, with applications
ranging from studies of neutron-star structure to tests of general relativity
and searches for nanosecond gravitational waves. As for other applications of
Gaussian distributions, outliers in timing measurements pose a significant
challenge to statistical inference, since they can bias the estimation of
timing and noise parameters, and affect reported parameter uncertainties. We
describe and demonstrate a practical end-to-end approach to perform Bayesian
inference of timing and noise parameters robustly in the presence of outliers,
and to identify these probabilistically. The method is fully consistent (i.e.,
outlier-ness probabilities vary in tune with the posterior distributions of the
timing and noise parameters), and it relies on the efficient sampling of the
hierarchical form of the pulsar-timing likelihood. Such sampling has recently
become possible with a "no-U-turn" Hamiltonian sampler coupled to a highly
customized reparametrization of the likelihood; this code is described
elsewhere, but it is already available online. We recommend our method as a
standard step in the preparation of pulsar-timing-array datasets: even if
statistical inference is not affected, follow-up studies of outlier candidates
can reveal unseen problems in radio observations and timing measurements;
furthermore, confidence in the results of gravitational-wave searches will only
benefit from stringent statistical evidence that datasets are clean and
outlier-free.
