  Gravitational waves and electromagnetic signals from merging neutron star
binaries provide valuable information about the the properties of dense matter,
the formation of heavy elements, and high-energy astrophysics. To fully
leverage observations of these systems, we need numerical simulations that
provide reliable predictions for the properties of the matter unbound in these
mergers. An important limitation of current simulations is the use of
approximate methods for neutrino transport that do not converge to a solution
of the transport equations as numerical resolution increases, and thus have
errors that are impossible to quantify. Here, we report on a first simulation
of a binary neutron star merger that uses Monte-Carlo techniques to directly
solve the transport equations in low-density regions. In high-density regions,
we use approximations inspired by implicit Monte-Carlo to greatly reduce the
cost of simulations, while only introducing errors quantifiable through more
expensive convergence studies. We simulate an unequal mass neutron star binary
merger up to $5\,{\rm ms}$ past merger, and report on the properties of the
matter and neutrino outflows. Finally, we compare our results to the output of
our best approximate `M1' transport scheme, demonstrating that an M1 scheme
that carefully approximates the neutrino energy spectrum only leads to $\sim
10\%$ uncertainty in the composition and velocity of the ejecta, and $\sim20\%$
uncertainty in the $\nu_e$ and $\bar\nu_e$ luminosities and energies. The most
significant disagreement found between M1 and Monte-Carlo results is a factor
of $\sim 2$ difference in the luminosity of heavy-lepton neutrinos.
