We present a new approach for modelling galaxy/halo bias that utilizes the full non-linear information contained in the moments of the matter density field , which we derive using a set of numerical simulations . Although our method is general , we perform a case study based on the local Eulerian bias scheme truncated to second-order . Using 200 N -body simulations covering a total comoving volume of 675 h ^ { -3 } { Gpc } ^ { 3 } , we measure several 2 - and 3 -point statistics of the halo distribution to unprecedented accuracy . We use the bias model to fit the halo-halo power spectrum , the halo-matter cross spectrum and the corresponding three bispectra for wavenumbers in the range 0.04 \lesssim k \lesssim 0.12 h { Mpc } ^ { -1 } . We find the constraints on the bias parameters obtained using the full non-linear information differ significantly from those derived using standard perturbation theory at leading order . Hence , neglecting the full non-linear information leads to biased results for this particular scale range . We also test the validity of the second-order Eulerian local biasing scheme by comparing the parameter constraints derived from different statistics . Analysis of the halo-matter cross-correlation coefficients defined for the 2- and 3-point statistics reveals further inconsistencies contained in the second-order Eulerian bias scheme , suggesting it is too simple a model to describe halo bias with high accuracy .