A scheme for calculating thermally averaged observables for quantum dissipative systems is presented. The method is based on a wavefunction with equal amplitude and random phase composed of a complete set of states, which is then propagated in imaginary time beta/2. Application to a Surrogate Hamiltonian simulation of a molecule subject to an ultrafast pulse coupled to a bath is studied. Compared to Boltzmann thermal averaging the method scales more favorably with an increase in the number of bath modes. A self-averaging phenomenon was identified which reduces the number of random sets required to converge the thermal average. (C) 2003 Elsevier B.V. All rights reserved.