The finite-temperature properties of Bose quantum fluids axe studied by applying a variational density-matrix approach. Our ansatz for the density matrix is similar to that introduced by Campbell, Kurten, Ristig, and Senger [Phys. Rev. 30, 3728 (1984)] but has been generalized to include three-body correlations. An expression for the finite-temperature elementary-excitation energy, epsilon(k), is derived from an ansatz for the variational entropy. At zero temperature epsilon(k) reduces to the same expression obtained by Brillouin-Wigner perturbation theory. The hypernetted-chain (HNC) equations provide the relations between the n-particle distribution functions and the correlation functions. A substantial set of elementary diagrams have been included to minimize the loss of precision associated with the HNC approximation. We report values for the static structure function, radial distribution function, excitation energies, particle-hole interaction, Helmholtz free energy, entropy, and isothermal sound velocity. We present a detailed study of the spectra, three-body correlations, and elementary diagram dependence of the thermodynamic spinodal curve. For relevant densities, we find that three-body correlations tend to increase the liquid-gas critical temperature T(c) while adding elementary diagrams tends to reduce T(c). Including multiphonon terms in the spectrum produces several features in the variational entropy that were not observed at the Feynman level of excitations. Entropy isotherms show a rather abrupt transition from a low-density regime, where multiphonon contributions are insignificant, to a high-density regime where the isotherms are strongly influenced by multiphonon scattering processes.