A variational Monte Carlo method for bosonic lattice models is introduced. The method is based on the Baeriswyl projected wave function. The Baeriswyl wave function consists of a kinetic energy based projection applied to the wave function at infinite interaction, and is related to the shadow wave function already used in the study of continuous models of bosons. The wave function at infinite interaction, and the projector, are represented in coordinate space, leading to an expression for expectation values which can be evaluated via Monte Carlo sampling. We calculate the phase diagram and other properties of the bosonic Hubbard model. The calculated phase diagram is in excellent agreement with known quantum Monte Carlo results. We also analyze correlation functions.