Quantum Monte Carlo (QMC) methods such as Variational Monte Carlo, Diffusion Monte Carlo or Path Integral Monte Carlo are the most accurate and general methods for computing total electronic energies. Taking many-body hydrogen at high pressure as an example, each method has a limited range of applicability, particularly at finite temperature. We have introduced a method to perform a coupled QMC for the electrons and another MC simulation for the ions (CEIMC). Using quantum Monte Carlo, one estimates the Born-Oppenheimer energy which is then used in a Metropolis simulation of the ionic degrees of freedom. We have shown that one can modify the usual Metropolis acceptance probability to eliminate the bias caused by noise in this energy difference, thus allowing more noisy estimates of the energy difference and thereby drastically reducing the sampling time of the electronic degrees of freedom. We introduce an optimal importance sampling to compute the needed energy difference. We have performed simulations of liquid H2 and metallic H on a parallel computer. For the atomic case, we introduce a new trial function, to avoid the necessity of determining the LDA orbitals as the protons are moved. There are some advantages of the CEIMC method relative to other simulations concerning how the quantum effects of the ionic degrees of freedom can be included and how the boundary conditions on the phase of the wavefunction can be integrated over.