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-dependentdensity matrix and for the mass operator. This allows us to develop a series of approximations for the operator in terms of exact propagators. The crudest approximation is a coherent potential approximation applicable at a general temperature