We present a novel Path Integral Monte Carlo scheme to solve the Fr\"ohlich polaron model. At intermediate and strong electron-phonon coupling, the polaron self-trapping is properly taken into account at the level of an effective action obtained by a preaveraging procedure with a retarded trial action. We compute the free energy at several couplings and temperatures in three and two dimensions. Our results show that the accuracy of the Feynman variational upper bound for the free energy is always better than 5% although the thermodynamics derived from it is not correct. Our estimates of the ground state energies demonstrate that the second cumulant correction to the variational upper bound predicts the self energy to better than 1% at intermediate and strong coupling.