We study how to numerically simulate quantum fermions out of thermal equilibrium, in the context of electroweak baryogenesis. We find that by combining the lattice implementation of Aarts and Smit [1] with the low cost fermions of Borsanyi and Hindmarsh [2], we are able to describe the dynamics of a classical bosonic system coupled to quantum fermions, that correctly reproduces anomalous baryon number violation. To demonstrate the method, we apply it to the 1+1 dimensional axial U(1) model, and perform simulations of a fast symmetry breaking transition. Compared to solving all the quantum mode equations as in [1], we find that this statistical approach may lead to a significant gain in computational time, when applied to 3+1 dimensional physics.