The recently introduced Gaussian Process State (GPS) provides a highly flexible, compact and physically insightful representation of quantum many-body states based on ideas from the zoo of machine learning approaches. In this work, we give a comprehensive description how such a state can be learned from given samples of a potentially unknown target state and show how regression approaches based on Bayesian inference can be used to compress a target state into a highly compact and accurate GPS representation. By application of a type II maximum likelihood method based on Relevance Vector Machines (RVM), we are able to extract many-body configurations from the underlying Hilbert space which are particularly relevant for the description of the target state, as support points to define the GPS. Together with an introduced optimization scheme for the hyperparameters of the model characterizing the weighting of modelled correlation features, this makes it possible to easily extract physical characteristics of the state such as the relative importance of particular correlation properties. We apply the Bayesian learning scheme to the problem of modelling ground states of small Fermi-Hubbard chains and show that the found solutions represent a systematically improvable trade-off between sparsity and accuracy of the model. Moreover, we show how the learned hyperparameters and the extracted relevant configurations, characterizing the correlation of the wavefunction, depend on the interaction strength of the Hubbard model as well as the target accuracy of the representation.