  We introduce a method for performing a robust Bayesian analysis of
non-Gaussianity present in pulsar timing data, simultaneously with the pulsar
timing model, and additional stochastic parameters such as those describing red
spin noise and dispersion measure variations. The parameters used to define the
presence of non-Gaussianity are zero for Gaussian processes, giving a simple
method of defining the strength of non-Gaussian behaviour. We use simulations
to show that assuming Gaussian statistics when the noise in the data is drawn
from a non-Gaussian distribution can significantly increase the uncertainties
associated with the pulsar timing model parameters. We then apply the method to
the publicly available 15 year Parkes Pulsar Timing Array data release 1
dataset for the binary pulsar J0437$-$4715. In this analysis we present a
significant detection of non-Gaussianity in the uncorrelated non-thermal noise,
but we find that it does not yet impact the timing model or stochastic
parameter estimates significantly compared to analysis performed assuming
Gaussian statistics. The methods presented are, however, shown to be of
immediate practical use for current European Pulsar Timing Array (EPTA) and
International Pulsar Timing Array (IPTA) datasets.
