A new method is proposed for the calculation of full density matrix and thermodynamic functions of a many-boson system. Explicit expressions are obtained in the pair correlations approximation for an arbitrary temperature. The theory is self-consistent in the sence that the calculated properties at low temperatures coincide with that of Bogoliubov theory and in the high-temperature limit lead to the results for classical non-ideal gas in the random phase approximation. The phase transition is also revealed as a concequence of Bose-Einstein condensation deformed by interatomic interactions. All the final formulae are written solely via the liquid structure factor taken as a source information instead of the interatomic potential and, therefore, interconnect only observable quantities. This gives also a possibility to study such a strongly non-ideal system as liquid 4 He.