Path integral hybrid Monte Carlo (PIHMC) algorithm for strongly correlated Bose fluids has been developed. This is an extended version of our previous method [S. Miura and S. Okazaki, Chem. Phys. Lett. 308, 115 (1999)] applied to a model system consisting of noninteracting bosons. Our PIHMC method for the correlated Bose fluids is constituted of two trial moves to sample path-variables describing system coordinates along imaginary time and a permutation of particle labels giving a boundary condition with respect to imaginary time. The path-variables for a given permutation are generated by a hybrid Monte Carlo method based on path integral molecular dynamics techniques. Equations of motion for the path-variables are formulated on the basis of a collective coordinate representation of the path, staging variables, to enhance the sampling efficiency. The permutation sampling to satisfy Bose-Einstein statistics is performed using the multilevel Metropolis method developed by Ceperley and Pollock [Phys. Rev. Lett. 56, 351 (1986)]. Our PIHMC method has successfully been applied to liquid helium-4 at a state point where the system is in a superfluid phase. Parameters determining the sampling efficiency are optimized in such a way that correlation among successive PIHMC steps is minimized. (C) 2004 American Institute of Physics.