  A technique to simulate the grand canonical ensembles of interacting Bose
gases is presented. Results are generated for many temperatures by averaging
over energy-weighted stochastic paths, each corresponding to a solution of
coupled Gross-Pitaevskii equations with phase-noise. The stochastic gauge
method used relies on an off-diagonal coherent-state expansion, thus taking
into account all quantum correlations. As an example, the second order spatial
correlation function and momentum distribution for an interacting 1D Bose gas
are calculated.
