The XMM-LSS, XMM-COSMOS, and XMM-CDFS surveys are complementary in terms of sky coverage and depth. Together, they form a clean sample with the least possible variance in instrument effective areas and PSF. Therefore this is one of the best samples available to determine the 2-10 keV luminosity function of AGN and its evolution. The samples and the relevant corrections for incompleteness are described. A total of 2887 AGN is used to build the LF in the luminosity interval 10^42-10^46 erg/s, and in the redshift interval 0.001-4. A new method to correct for absorption by considering the probability distribution for the column density conditioned on the hardness ratio is presented. The binned luminosity function and its evolution is determined with a variant of the Page-Carrera method, improved to include corrections for absorption and to account for the full probability distribution of photometric redshifts. Parametric models, namely a double power-law with LADE or LDDE evolution, are explored using Bayesian inference. We introduce the Watanabe-Akaike information criterion (WAIC) to compare the models and estimate their predictive power. Our data are best described by the LADE model, as hinted by the WAIC indicator. We also explore the 15-parameter extended LDDE model recently proposed by Ueda et al., and find that this extension is not supported by our data. The strength of our method is that it provides: un-absorbed non-parametric estimates; credible intervals for luminosity function parameters; model choice according to which one has more predictive power for future data.