  We introduce a new type of Gutzwiller variational wavefunction for correlated
electrons coupled to phonons, able to treat on equal footing electronic and
lattice degrees of freedom. We benchmark the wavefunction in the infinite-$U$
Hubbard-Holstein model away from half-filling on a Bethe lattice, where we can
directly compare with exact results by Dynamical Mean-Field Theory. For this
model, we find that variational results agree perfectly well with the exact
ones. In particular the wavefunction correctly describes the crossover to a
heavy polaron gas upon increasing the electron-phonon coupling.
