New variational wave functions are introduced to describe the Hubbard model in strong correlation. Intersite correlation (or Jastrow) factors are taken into account beyond the Gutzwiller wave function (GF). A study is made on the one-dimensional lattice and the two-dimensional square lattice at half-filling as well as less-than-half-filling. Much improvement can be achieved thereby in the description of physical quantities such as the momentum distribution and spin and charge density correlation functions. In addition, the ground state energy is lowered appreciably compared with the GF. In this connection, the Mott-Hubbard-type metal-insulator transition is also re-examined.