  Gaussian Processes (GP) have become popular machine-learning methods for
kernel-based learning on datasets with complicated covariance structures. In
this paper, we present a novel extension to the GP framework using a
contaminated normal likelihood function to better account for heteroscedastic
variance and outlier noise. We propose a scalable inference algorithm based on
the Sparse Variational Gaussian Process (SVGP) method for fitting sparse
Gaussian process regression models with contaminated normal noise on large
datasets. We examine an application to geomagnetic ground perturbations, where
the state-of-the-art prediction model is based on neural networks. We show that
our approach yields shorter prediction intervals for similar coverage and
accuracy when compared to an artificial dense neural network baseline.
