The Fourier coefficient path-integral representation of the quantum density matrix is used to carry out direct, accurate calculations of coupled vibrational partition functions. The present implementation of the Fourier path-integral method incorporates two noteworthy features. First, we use a Gaussian in Fourier space as a probability density function, which is sampled using the Box-Muller algorithm. Second, we introduce an adaptively optimized stratified sampling scheme in Cartesian coordinates to sample the nuclear configurations. We illustrate these strategies by applying them to a coupled stretch-bend model which resembles two of the vibrational modes of H 2O. We also apply a simple, yet accurate method for estimating the statistical error of a Metropolis integration, and we compare the Box-Muller and Metropolis sampling algorithms in detail. The numerical tests of the new method are very encouraging, and the approach is promising for accurate calculations of quantum free energies for polyatomic molecules.