Though half of cosmic starlight is absorbed by dust and reradiated at long wavelengths (3μ m – 3 mm), constraints on the infrared through millimeter galaxy luminosity function (the ‘IRLF’) are poor in comparison to the rest-frame ultraviolet and optical galaxy luminosity function, particularly at z ⩾ 2.5. Here we present a backward evolution model for interpreting number counts, redshift distributions, and cross-band flux density correlations in the infrared and submillimeter sky, from 70μm – 2 mm, using a model for the IRLF out to the epoch of reionization. Mock submillimeter maps are generated by injecting sources according to the prescribed IRLF and flux densities drawn from model spectral energy distributions that mirror the distribution of SEDs observed in 0 < z 0 < 5 dusty star-forming galaxies (DSFGs). We explore two extreme hypothetical case-studies: a dust-poor early Universe model, where DSFGs contribute negligibly (< 10%) to the integrated star-formation rate density at z > 4, and an alternate dust-rich early Universe model, where DSFGs dominate > 90% of z > 4 star-formation. We find that current submm/mm datasets do not clearly rule out either of these extreme models. We suggest that future surveys at 2 mm – both from ALMA and single-dish facilities – will be crucial to measuring the IRLF beyond z > 4.