I think some information should be added to the previous answers, which contain good general points but miss two more specific issues relevant for answering your question.
The question about how large has to be $N$ for using statistical mechanics with confidence certainly doesn't have a sharp number as an answer, and indeed, the use of statistical mechanics for the ordinary matter in a lab is justified by numbers of degrees of freedom of the order of Avogadro number.
However, experiments on small atomic clusters and more than half a century of computer simulation have shown that the real answer may be more subtle.
First of all, we have to distinguish between the values of $N$ providing a good approximation for the thermodynamic limit ($N\rightarrow \infty$) and the values of $N$ enabling us to use statistical methods safely.
Indeed, the behavior of systems as small as a few tens of atoms may provide a satisfactory qualitative and quantitative description of thermodynamic properties almost everywhere in the phase diagram, provided the effect of boundary surfaces is minimized (for example, through periodic boundary conditions). The key lesson is that an important parameter is the correlation length. The size of the system has to be larger than the correlation length $\xi$. Except for continuous phase transitions, the number of degrees of freedom inside a volume $\xi^D$, where $D$ is the dimensionality of the space, is not huge and much below the Avogadro number.
We have also to recall that a statistical description of the behavior of small systems is possible independently of the problem of approximating the behavior of bulk systems. In that case, we can continue using the usual ensembles, ut we need to choose the most adapted to the physical conditions.
A second point, specific to gravitational systems, is that application of statistical mechanics methods should take into account the peculiarities of long-range gravitational systems (non-extensiveness of the energy and entropy and negative specific heat). The presence of such oddities does not hamper the use of statistical mechanics methods but requires a more careful analysis than usual. Also in such a case, the usual equivalence of the ensembles at the thermodynamic limit is lost. Also, the usual way of analyzing the size dependence of fluctuations requires a deep revision. Still, a statistical mechanics description of stellar systems can be done safely even with numbers much below the Avogadro number.