The density matrix, rho, of a model polariton system is obtained numerically from a master equation which takes account of pumping and losses. In the stationary limit, the coherences between eigenstates of the Hamiltonian are three orders of magnitude smaller than the occupations, meaning that the stationary density matrix is approximately diagonal in the energy representation. A weakly distorted grand canonical Gibbs distribution fits well the occupations.