Gaussian process regression is widely applied in computational science and engineering for surrogate modeling owning to its kernel-based and probabilistic nature. In this work, we propose a Bayesian approach that integrates the variability of input data into the Gaussian process regression for function and partial differential equation approximation. Leveraging two types of observables -- noise-corrupted outputs with certain inputs and those with prior-distribution-defined uncertain inputs, a posterior distribution of uncertain inputs is estimated via Bayesian inference. Thereafter, such quantified uncertainties of inputs are incorporated into Gaussian process predictions by means of marginalization. The setting of two types of data aligned with common scenarios of constructing surrogate models for the solutions of partial differential equations, where the data of boundary conditions and initial conditions are typically known while the data of solution may involve uncertainties due to the measurement or stochasticity. The effectiveness of the proposed method is demonstrated through several numerical examples including multiple one-dimensional functions, the heat equation and Allen-Cahn equation. A consistently good performance of generalization is observed, and a substantial reduction in the predictive uncertainties is achieved by the Bayesian inference of uncertain inputs.