Recent work on brain tumor growth modeling for glioblastoma using reaction-diffusion equations suggests that the diffusion coefficient and the proliferation rate can be related to clinically relevant information. However, estimating these parameters is difficult due to the lack of identifiability of the parameters, the uncertainty in the tumor segmentations, and the model approximation, which cannot perfectly capture the dynamics of the tumor. Therefore, we propose a method for conducting the Bayesian personalization of the tumor growth model parameters. Our approach estimates the posterior probability of the parameters, and allows the analysis of the parameters correlations and uncertainty. Moreover, this method provides a way to compute the evidence of a model, which is a mathematically sound way of assessing the validity of different model hypotheses. Our approach is based on a highly parallelized implementation of the reaction-diffusion equation, and the Gaussian Process Hamiltonian Monte Carlo (GPHMC), a high acceptance rate Monte Carlo technique. We demonstrate our method on synthetic data, and four glioblastoma patients. This promising approach shows that the infiltration is better captured by the model compared to the speed of growth.