We present the electron-phonon averaged via Gaussian process regression (EPA-GPR) method, in which the electron-phonon coupling matrix is represented as a function of two energies and is in turn modeled as a Gaussian process. The EPA-GPR method can be used as an efficient method to estimate thermoelectric properties of materials for fast-screening applications, comparable to the original electron-phonon averaged (EPA) method and the electron-phonon averaged via moving-least-squares (EPA-MLS) method. The EPA-GPR method does not require specification of any open parameter, unlike the other EPA-related methods, since all the hyperparameters in the model can be unambiguously estimated within the type II maximum likelihood (ML-II) approximation. Thus, the EPA-GPR method is a parameter-free estimation method. Additionally, the concept of Gaussian processes in the EPA-GPR method allows us to quantify the uncertainty in estimated properties of thermoelectric materials. One can randomly realize the electron-phonon coupling coefficients from the identified Gaussian process, and those realized samples can be further analyzed in the solution process of the semiclassical Boltzmann transport equation for charge carriers. The results of the semiclassical Boltzmann transport equation provide the statistical properties of the thermoelectric properties of interest. The means, standard deviations, histograms, and confidence intervals of the Seebeck coefficient, the electrical conductivity, and the power factor can be constructed and analyzed. Furthermore the proposed EPA-GPR method is applied to a p-type half-Heusler compound, i.e., HfCoSb, as a case example, the results of which clearly present the advantages of the method.