© 2020 authors. Published by the American Physical Society. Published by the American Physical Society under the terms of the "https://creativecommons.org/licenses/by/4.0/"Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. We present a novel, nonparametric form for compactly representing entangled many-body quantum states, which we call a "Gaussian process state."In contrast to other approaches, we define this state explicitly in terms of a configurational data set, with the probability amplitudes statistically inferred from this data according to Bayesian statistics. In this way, the nonlocal physical correlated features of the state can be analytically resummed, allowing for exponential complexity to underpin the ansatz, but efficiently represented in a small data set. The state is found to be highly compact, systematically improvable, and efficient to sample, representing a large number of known variational states within its span. It is also proven to be a "universal approximator"for quantum states, able to capture any entangled many-body state with increasing data-set size. We develop two numerical approaches which can learn this form directly - a fragmentation approach and direct variational optimization - and apply these schemes to the fermionic Hubbard model. We find competitive or superior descriptions of correlated quantum problems compared to existing state-of-the-art variational ansatzes, as well as other numerical methods.