Modeling epidemic spread is critical for informing policy decisions aimed at mitigation. Accordingly, in this work we present a new data-driven method based on Gaussian process regression (GPR) to model epidemic spread through the difference on the logarithmic scale of the infected cases. We bound the variance of the predictions made by GPR, which quantifies the impact of epidemic data on the proposed model. Next, we derive a high-probability error bound on the prediction error in terms of the distance between the training points and a testing point, the posterior variance, and the level of change in the spreading process, and we assess how the characteristics of the epidemic spread and infection data influence this error bound. We present examples that use GPR to model and predict epidemic spread by using real-world infection data gathered in the UK during the COVID-19 epidemic. These examples illustrate that, under typical conditions, the prediction for the next twenty days has 94.29% of the noisy data located within the 95% confidence interval, validating these predictions. We further compare the modeling and prediction results with other methods, such as polynomial regression, k-nearest neighbors (KNN) regression, and neural networks, to demonstrate the benefits of leveraging GPR in disease spread modeling.