The momentum distribution function of a parabolically confined gas of bosons with harmonic interparticle interactions is derived. In the Bose-Einstein condensation region, this momentum distribution substantially deviates from a Maxwell-Boltzmann distribution. It is argued that the determination of the temperature of the boson gas from the Bose-Einstein momentum distribution function is more appropriate than the currently used fitting to the high momentum tail of the Maxwell-Boltzmann distribution.