The seismic data is often inadequately sampled in the spatial dimension, which severely distorts the results of seismic imaging. Seismic data interpolation is an important approach to recover the missed data. We propose a numerical scheme based on Gaussian process regression to interpolate the inadequately sampled seismic data. The observed seismic traces are treated as the training data to estimate the unknown seismic traces through Gaussian process regression. The proposed scheme not only predicts the unknown traces, but also provides the uncertainty quantification of the predictions. We also compare the proposed scheme with the POCS interpolation method through synthetic and field data. The numerical results demonstrate that the inadequately sampled synthetic and field data are effectively reconstructed by the proposed scheme.