The canonical one-band Hubbard model is studied using a computational method that mixes the Monte Carlo procedure with the mean field approximation. This technique allows us to incorporate thermal fluctuations and the development of short-range magnetic order above ordering temperatures, contrary to the crude finite-temperature Hartree-Fock approximation, which incorrectly predicts a N\'eel temperature $T_N$ that grows linearly with the Hubbard $U/t$. The effective model studied here contains quantum and classical degrees of freedom. It thus belongs to the "spin fermion" model family widely employed in other contexts. Using exact diagonalization, supplemented by the traveling cluster approximation, for the fermionic sector, and classical Monte Carlo for the classical fields, the Hubbard $U/t$ vs. temperature $T/t$ phase diagram is studied employing large three and two dimensional clusters. We demonstrate that the method is capable of capturing the formation of local moments in the normal state without long-range order, the non-monotonicity of $T_N$ with increasing $U/t$, the development of gaps and pseudogaps in the density of states, and the two-peak structure in the specific heat. Extensive comparisons with determinant quantum Monte Carlo results suggest that the present approach is qualitatively, and often quantitatively, accurate, particularly at intermediate and high temperatures. Finally, we study the Hubbard model including plaquette diagonal hopping (i.e. the $t-t^\prime$ Hubbard model) in two dimensions and show that our approach allows us to study low temperature properties where determinant quantum Monte Carlo fails due to the fermion sign problem. Future applications of this method include multi-orbital Hubbard models such as those needed for iron-based superconductors.