We consider the density matrices that arise in the statistical mechanics of the electron-phonon systems. In the path integral representation the phonon coordinates can be eliminated. This leads to an action that depends on pairs of points on a path, that depends explicitly on time differences, and that contains the phonon occupation numbers. The integral is reduced to a standard form by scaling to the thermal length. We use the technique of integration by parts and add specially chosen generating functionals to the action. We set down functional derivative equations for the source-dependent density matrix and for the mass operator.