Most previous analyses of the astrophysical distribution of merging binary black hole systems have focused on fitting parametric phenomenological models to the observed data (e.g., [1]} and references therein) or directly compared with detailed simulations (e.g., [2]}, [3]}).
However, non-parametric methods have been employed to infer the binary black hole mass [4]}, [5]}, [6]}, [7]} and spin distributions [5]}.
