I introduce, for one degree-of-freedom harmonic oscillator stationary states, a modification of the quasi-classical approximation for expectation values of observables which are polynomial functions of the position and momentum. After time averaging a suitable function determined by the observable, instead of evaluating the average at the energy eigenvalue, if one then averages over the energy with a certain gamma distribution, one gets more accurate expectation values for large quantum numbers. I consider Monte Carlo evaluations of the averages with the random phase points chosen according to the respective distributions. Although the Monte Carlo error is greater for the modification than for the quasi-classical approximation when the same number of phase points is used, it is less than that for a Monte Carlo evaluation of a certain exact expression in which the same time averaged function is averaged over the energy with Wigner's quasiprobability distribution.