A statistical sampling method is proposed for computing oscillatory integrals associated with the semiclassical initial value representation. The semiclassical expression is rewritten as an integral over a phase distribution P(s). The phase distribution is obtained from Metropolis sampling of trajectories according to a properly chosen weight function. The averaging of oscillatory integrals is converted into a Monte Carlo algorithm where one diffuses through trajectory space. A histogram of phases is collect from importance sampling. Techniques of Metropolis Monte Carlo such as umbrella (or biased) sampling are generalized to the present context. From example calculations, phase distributions are seen to be multi-peaked, thus clearly demonstrating the origin of quantum interference. Trajectories that are responsible for the interference patterns can be collected using this method.