We present a Bayesian reconstruction of parton distribution functions (PDFs) from Lattice QCD data using Gaussian Processes (GPs). Within the pseudo-distribution framework, matrix elements of bilocal operators, computed on the lattice, are related to PDFs via a Fourier transform. This gives rise to an ill-posed inverse problem, as we aim to infer a smooth, continuous function from a finite set of noisy and correlated data points. Gaussian Processes provide a robust and flexible tool to model this uncertainty, allowing for the inclusion of prior physical knowledge such as positivity, asymptotic behavior, or sum rules. This approach opens a path toward interpretable and uncertainty-aware PDF extraction, grounded in modern probabilistic inference.