Ferrimagnets are described by molecular field models, which operate by summing sublattice magnetic moments with Brillouin functions. They are commonly initialized with the fully saturated ("spontaneous") moment M_s(T), obtained by extrapolating the high-field portion of M(H) hysteresis curves to zero-field. However, the bulk polycrystalline ferrimagnetic samples used in our experiments were measured in their remanent state, so the directly observed moment M_r(T) can differ significantly from M_s(T), leading to model and experiment mismatches of up to ~3x. We derive a practical connection between M_s(T) from the molecular field model and the experimental M_r(T). Our approach introduces a correction factor β(T) = M_r/M_s derived from hysteresis curves of our sample that can be folded back into the molecular field fit. Applying this correction to polycrystalline terbium iron garnet yields good agreement between the modeled and measured magnetization in the relevant temperature range, especially near the compensation temperature, which is where the richest physics resides. In addition, the corrected model reproduces the measured neutron spin rotation of our samples that depend on the internal magnetization of the sample. This explicit treatment of the remanent state reconciles the molecular field model with measurements on actual rare-earth iron garnet targets, providing a backbone to the results of the NSR-Ferrimagnets collaboration and other exotic force-searching experiments using ferrimagnets.