The accuracy of some density functional (DF) models widely used in material science depends on empirical or free parameters that are commonly tuned using reference physical properties. Grid-search methods are the standard numerical approximations used to find the optimal values of the free parameters, making the computational complexity to scale with the number of points in the grid. In this report, we illustrate that Bayesian optimization (BO), a sample-efficient machine learning algorithm, can calibrate different density functional models, e.g., hybrid exchange-correlation and range-separated density functionals. Using the atomization energies and bond lengths from the Gaussian-1 (G1) and Gaussian-2 (G2) databases, we show that BO optimizes the free parameters of the hybrid exchange-correlation functionals with approximately SS evaluations of the error function. We also show that selecting exchange-correlation functionals for different physical systems can be done with BO. We jointly optimize and select the form of the exchange-correlation functionals and the free parameters by minimizing the root-mean-square error functions for the G1 and G2 data set atomization energies. The calibrated DF model is more accurate on average than standard DF methods, e.g., PBE0 and B3LYP.