A robust uncertainty estimate in global analyses of Parton Distribution Functions (PDFs) is essential at the Large Hadron Collider (LHC), especially in view of the high-precision data anticipated by experimentalists in the High-Luminosity phase of the LHC. A Bayesian framework to determine PDFs provides a rigorous treatment of uncertainties and full control on the prior, though its practical implementation can be computationally demanding. To address these challenges, we introduce a novel approach to PDF determination tailored for Bayesian inference, based on the use of linear models. Unlike traditional parametrisations, our method represents PDFs as vectors in a functional space spanned by specially chosen bases, derived from the dimensional reduction of a neural network functional space, providing a compact yet versatile representation of PDFs. The low-dimensionality of the preferred models allows for particularly fast inference. The size of the bases can be systematically adjusted, allowing for transparent control over underfitting and overfitting, and facilitating principled model selection through Bayesian workflows. In this work, the methodology is applied to a fit of Deep Inelastic Scattering synthetic data, and thoroughly tested via multi-closure tests, thus paving the way to its application to global PDF fits.