This paper presents a formal foundation for the recent extension [J. Chem. Phys. 110, 3647 (1999)] of the centroid molecular dynamics (CMD) method to systems obeying Bose-Einstein statistics. It is shown that the introduction of centroid phase space coordinates corresponding to individual physical particles allows one to obtain (exact) canonical averages within the framework of the bosonic CMD method. It is also shown that formally exact expressions for quantum mechanical Kubo transformed correlation functions can be written in terms of individual particle centroids and that a CMD approximation can be introduced. Calculations for a bosonic trimer are used as an illustration of the new concepts introduced in this work.