Variational Monte Carlo method for fermionic models combined with tensor networks and applications to the hole-doped two-dimensional Hubbard model