Ok, I dug out our old stat mech/thermo textbook. YES, Maxwell-Boltzmann statistics definitely apply to stars in a globular cluster or galaxy, but you have have to pare back the results to the absolute most general.
Sears and Salinger go through an excellent derivation of Maxwell-Boltzmann statistics as well as the Maxwell-Boltzmann distribution function. The most general results leave the distribution function as a function of completely unspecified energy levels, and everything I saw looks like it is absolutely applicable to big astrophysical things (e.g. stars) clustering together into even bigger astrophysical things (e.g. globular clusters, galaxies).
I had some qualms, though. The energy levels are quantized in their treatment. And the particles are non-interacting. However, it looks like M-B stats are still applicable. To put the icing on the cake, later in the text they go through a derivation of the usual ideal thermodynamic piston, except this one is IN A GRAVITATIONAL FIELD, BABY! They assume a uniform field, and use that to derive the Newtonian hydrostatic equation (my link is for the general relativistic generalization of the Newtonian; the Wikipedia link for Newton's version was unsatisfactory) from a pure thermodynamic point of view, which, as a professional physicist, made me nearly literally stand up and clap.
So- at the very least, I can see someone assuming that stars are floating around passively, i.e. not interacting gravitationally with each other, but subject to a magical, arbitrary gravitational field matching the real-life solution. They would then derive a purely thermodynamic equation showing the distribution of stars in that magical gravitational field. Then, they would calculate the gravitational field of stars distributed how they just calculated. Then, they would show that the generated gravitational field matches the original, magical, arbitrary field.
This is called generating a self-consistent solution. I'd have to check with a mathematician, but I believe the equations you would use have only a single solution, so even if you kinda-sorta cheated in solving them, your solution is still The solution.
If you were even smarter, you might be able to generalize Sears and Salinger's gravity piston to self-gravitating particles and derive the solution directly. Not sure that's possible, but maybe.
I don't think that an expanding Universe can be considered in thermal equilibrium, except on short time scales. I mean, the CMB started out in short-term equilibrium, then just look at what happened!