Theory uncertainty from missing higher orders is usually estimated using scale variation. While scale variation is certainly a good tool to guess the size of the next perturbative order, it lacks of a statistical (probabilistic) interpretation and it often underestimates the actual uncertainty. Having a statistical definition of theory uncertainties is not only need for a fair comparison with data in precision physics: it is also a useful tool for extracting PDFs accounting also for theory uncertainties. A few years ago Cacciari and Houdeau proposed a Bayesian model to give a statistical meaning to theory uncertainties. The idea was excellent, but it has some limitations, which in turn make it not very well performing, especially for LHC physics. In this talk I will present two new Bayesian models: one is an improved version of the Cacciari-Houdeau model, and the other is inspired by scale variation. I will further show how scale dependence can be removed from a finite-order result within the context of these models. As a proof of concept, I will apply the methods to inclusive Higgs production in gluon fusion, for which four perturbative orders are known, and which is characterised by large perturbative corrections. The results are very interesting and promising.