  Cosmological simulations are reaching the resolution necessary to study
ultra-faint dwarf galaxies. Observations indicate that in small populations,
the stellar initial mass function (IMF) is not fully populated; rather, stars
are sampled in a way that can be approximated as coming from an underlying
probability density function. To ensure the accuracy of cosmological
simulations in the ultra-faint regime, we present an improved treatment of the
IMF. We implement a self-consistent, stochastically populated IMF in
cosmological hydrodynamic simulations. We test our method using high-resolution
simulations of a Milky Way halo, run to $z=6$, yielding a sample of nearly 100
galaxies. We also use an isolated dwarf galaxy to investigate the resulting
systematic differences in galaxy properties. We find that a stochastic IMF in
simulations makes feedback burstier, strengthening feedback, and quenching star
formation earlier in small dwarf galaxies. For galaxies in halos with mass
$\lesssim10^{8.5}$ M$_\odot$, a stochastic IMF typically leads to lower stellar
mass compared to a continuous IMF, sometimes by more than an order of
magnitude. We show that existing methods of ensuring discrete supernovae
incorrectly determine the mass of the star particle and its associated
feedback. This leads to overcooling of surrounding gas, with at least
${\sim}$10 per cent higher star formation and ${\sim}$30 per cent higher cold
gas content. Going forward, to accurately model dwarf galaxies and compare to
observations, it will be necessary to incorporate a stochastically populated
IMF that samples the full spectrum of stellar masses.
