Do Bose-Einstein and Fermi-Dirac statistics hold for situations where the grand canonical density is inappropriate? The grand canonical density operator describes a system (that may or may not be macroscopic) which is in contact with a thermal reservoir and a particle reservoir (so that in particular $T$ and $\mu$ of the system is a well-defined constant).
The Bose-Einstein and Fermi-Dirac statistics, and their classical analogue Maxwell-Boltzmann statistics are derived from the grand canonical density.
However, I see many people apply these statistics to any system in thermal equilibrium regardless of the boundary conditions of the system.
This motivates my question: Do Bose-Einstein and Fermi-Dirac statistics hold for situations where the grand canonical density is inappropriate?
Here, I mean:
Bose-Einstein statistics : $\left<  \hat {n_i}  \right> =1/(e^{\beta ( \epsilon_k - \mu) } -1)$
Fermi-Dirac statistics = $\left<  \hat {n_i}  \right> =1/(e^{\beta ( \epsilon_k - \mu) } +1)$
Maxwell-Boltzmann statistics : $\left<  \hat {n_i}  \right> = e^{\beta (\mu - \epsilon_k)}$
 A: As usual in Statistical Mechanics, different ensembles provide different values for thermodynamic and statistical properties. It is only for a subset of the possible Hamiltonians and suitably large systems that these different values go to the same limit. The technical procedure is the so-called thermodynamical limit.
Once one knows that the thermodynamic limit exists for the Hamiltonian of interest, it is usually possible to show that different ensembles provide the same thermodynamic and statistical descriptions independently on the starting ensemble at the thermodynamic limit. This the case for the ideal gases (bosons or fermions).
We also note that even though the thermodynamic limit implies an infinite size, numerical deviations from the limit values can be very small, even for a few tens of particles if periodic boundary conditions are used.
The above considerations apply as well to the ideal gas both in the classical and quantum regimes. In the case of the ideal gas, it is possible to show explicitly how the difference of description in different ensembles vanishes at the thermodynamic limit.
Thus, in the usual introduction to the behavior of Fermi and Bose gases, the equivalence of the ensembles is taken for granted, and the most convenient ensemble for calculations is chosen. In the case of quantum statistics, the Grand-canonical ensemble gets rid of the constraint about the number of particles in the summations simplifying the mathematical manipulations. However, it is possible to do calculations in the canonical ensemble and you may find indications and links in the answers to this related but different question.
Moreover, a derivation of the Bose-Einstein distribution also in the microcanonical ensemble can be found here.
Summarizing, the answer to your question is that in the case of systems made by a few particles we do expect deviations from the grand-canonical expressions even for the average occupancy of a state. Those deviations can be expressed in terms of a finite-size correction to the chemical potential.
A: It is somewhat unclear what you are meaning to ask in your question, since you appear to have misunderstood the meaning of the terms "Bose-Einstein statistics" and "Fermi-Dirac statistics."  The "statistics" are fundamental properties of the particles; they are not limited to thermodynamics or statistical mechanics.  Bose-Einstein statistics says that identical particles with integer spin are indistinguishable, and any number of them can be in the same state ("bosons").  Fermi-Dirac statistics says that identical particles with half-integer spin are also indistinguishable, but no more than one of them can be in a given quantum state ("fermions").  That is all "Bose-Einstein statistics" and "Fermi-Dirac statistics" mean.
Using these statistical properties, it is possible to calculate boson and fermion distribution functions, using statistical mechanics.  The easiest way to do this, when the quantum statistics of the particles is important, is using the grand canonical ensemble.  The advantage of the grand canonical ensemble is that the calculation is almost completely automatic, with the thermodynamic properties following straight from the grand partition function.  However, you can get the same results in the canonical ensemble or the microcanonical ensemble, although it may take a little more work.  In the microcanonical case, you can find the same distribution functions by calculating the most likely distributions, then showing that for large systems, almost all possible distributions are effectively identical to the most likely ones.  So for a thermodynamically large system, the "Bose-Einstein" and "Fermi-Dirac distribution functions" are valid for any system of identical quanta in equilibrium, regardless of whether temperature $T$ is fixed or energy $E$ is fixed, and regardless of whether chemical potential $\mu$ is fixed or particle number $N$ is fixed.
A: In the vast majority of applications of thermodynamics, there is no near-infinite reservoir that sets the temperature or chemical potential; such a reservoir isn't needed, it's mostly a teaching tool.
For example, if I place a gas in a box so that the gas molecules collide elastically with the walls, then there can be no energy transfer to or from the walls. So the walls can't act as an energy reservoir, but that doesn't mean the gas's temperature isn't defined. Given enough time, it will approach the Maxwell-Boltzmann distribution for some temperature $T$. Intuitively, because the gas is macroscopic, each part of the gas effectively acts as a near-infinite reservoir for the rest of the gas. Energy transfers between different parts of the gas until the entropy is maximized, which corresponds to a uniform temperature. Similarly, particle flow between different parts of the box, yielding a uniform chemical potential.
Of course, this argument only works if the system you're considering is macroscopic. If you have a microscopic system, you absolutely do need to consider the boundary conditions (i.e. how it is coupled to its macroscopic surroundings, which act as a reservoir).
