An ambitious goal in cosmology is to forward model the observed distribution of galaxies in the nearby Universe today from the initial conditions of large-scale structures. For practical reasons, the spatial resolution at which this can be done is necessarily limited. Consequently, one needs a mapping between the density of dark matter averaged over similar to Mpc scales and the distribution of dark matter haloes (used as a proxy for galaxies) in the same region. Here, we demonstrate a method for determining the halo mass distribution function by learning the tracer bias between density fields and halo catalogues using a neural bias model. The method is based on the Bayesian analysis of simple, physically motivated, neural network-like architectures, which we denote as neural physical engines, and neural density estimation. As a result, we are able to sample the initial phases of the dark matter density field while inferring the parameters describing the halo mass distribution function, providing a fully Bayesian interpretation of both the initial dark matter density distribution and the neural bias model. We successfully run an upgraded BORG (Bayesian Origin Reconstruction from Galaxies) inference using our new likelihood and neural bias model with halo catalogues derived from full N-body simulations. In preliminary results, we notice there could potentially be orders of magnitude improvement in modelling compared to classical biasing techniques.