We present a method for calculating the infrared emission from a population of dust grains heated by starlight, including very small grains for which stochastic heating by starlight photons results in high transients. Because state-to-state transition rates are generally unavailable for complex molecules, we consider model PAH, graphitic, and silicate grains with realistic vibrational mode spectra and realistic radiative properties. The vibrational density of states is used in a statistical-mechanical description of the emission process. Unlike previous treatments, our approach fully incorporates multiphoton heating effects, important for large grains or strong radiation fields. We discuss how the temperature of the grain is related to its vibrational energy. By comparing with an exact statistical calculation of the emission process, we determine the conditions under which the thermal and the continuous cooling approximations can be used to calculate the emission spectrum. We present results for the infrared emission spectra of PAH grains of various sizes heated by starlight. We show how the relative strengths of the 6.2, 7.7, and 11.3um features depend on grain size, starlight spectrum and intensity, and grain charging conditions. We show results for grains in the cold neutral medium, warm ionized medium, and representative conditions in photodissociation regions. Our model results are compared to observed ratios of emission features for reflection nebulae and photodissociation regions, the Milky Way, normal spiral galaxies, and starburst galaxies.