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. For the first time, 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 10 8.5 M , 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 ∼10 per cent higher star formation and ∼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.With the recent successes in modeling dwarf galaxies, various groups are pushing the limits of resolution even further, into the ultra-faint dwarf (UFD) galaxy range (M star 10 5 M ; e.g.