  The properties of black-hole and neutron-star binaries are extracted from
gravitational-wave signals using Bayesian inference. This involves evaluating a
multi-dimensional posterior probability function with stochastic sampling. The
marginal probability density distributions from which the samples are drawn are
usually interpolated with kernel density estimators. Since most post-processing
analysis within the field is based on these parameter estimation products,
interpolation accuracy of the marginals is essential. In this work, we propose
a new method combining histograms and Gaussian Processes as an alternative
technique to fit arbitrary combinations of samples from the source parameters.
This method comes with several advantages such as flexible interpolation of
non-Gaussian correlations, Bayesian estimate of uncertainty, and efficient
re-sampling with Hamiltonian Monte Carlo.
