  We present a method that generalises the standard mean field theory of
correlated lattice bosons to include amplitude and phase fluctuations of the
$U(1)$ field that induces onsite particle number mixing. This arises formally
from an auxiliary field decomposition of the kinetic term in a Bose Hubbard
model. We solve the resulting problem, initially, by using a classical
approximation for the particle number mixing field and a Monte Carlo treatment
of the resulting bosonic model. In two dimensions we obtain $T_c$ scales that
dramatically improve on mean field theory and are within about 20% of full
quantum Monte Carlo estimates. The `classical approximation' ground state,
however, is still mean field, with an overestimate of the critical interaction,
$U_c$, for the superfluid to Mott transition. By further including low order
quantum fluctuations in the free energy functional we improve significantly on
the $U_c$, and the overall thermal phase diagram. The classical approximation
based method has a computational cost linear in system size. The methods
readily generalise to multispecies bosons and the presence of traps.
