  A recently developed method, introduced in Phys. Rev. Lett. 94 (2005) 180403,
Phys. Rev. B 72 (2005) 064302, Phys. Lett. A 344 (2005) 84, systematically
improved the convergence of generic path integrals for transition amplitudes.
This was achieved by analytically constructing a hierarchy of $N$-fold
discretized effective actions $S^{(p)}_N$ labeled by a whole number $p$ and
starting at $p=1$ from the naively discretized action in the mid-point
prescription. The derivation guaranteed that the level $p$ effective actions
lead to discretized transition amplitudes differing from the continuum limit by
a term of order $1/N^p$. Here we extend the applicability of the above method
to the calculation of energy expectation values. This is done by constructing
analytical expressions for energy estimators of a general theory for each level
$p$. As a result of this energy expectation values converge to the continuum as
$1/N^p$. Finally, we perform a series of Monte Carlo simulations of several
models, show explicitly the derived increase in convergence, and the ensuing
speedup in numerical calculation of energy expectation values of many orders of
magnitude.
