What is the relationship between Maxwell–Boltzmann statistics and the grand canonical ensemble? In the grand canonical ensemble one derives the expectation value $\langle \hat n_r\rangle^{\pm}$ for fermions and bosons of sort $r$:
$$ \langle \hat n_r\rangle^{\pm} \ \propto \ \frac{1}{\mathrm{exp}[(\varepsilon_r-\mu)/k_B T] \mp 1} . $$
For $(\varepsilon_r-\mu) / k_B T\gg 0$, we find 
$$ \langle \hat n_r\rangle^{\pm} \ \approx \ \frac{1}{\mathrm{exp}[(\varepsilon_r-\mu)/k_B T]} \ \propto \mathrm{exp}[-(\varepsilon_r-\mu)/k_B T].$$
The same motivation seems to be found in this Wikipedia article. However, on the same page, right at the beginning, that intuitive statement is made: 

In statistical mechanics, Maxwell–Boltzmann statistics describes the average distribution of non-interacting material particles over various energy states in thermal equilibrium, and is applicable when the temperature is high enough or the particle density is low enough to render quantum effects negligible.

Now from my derivation above, it seems that "temperature is high enough" does the opposite of helping $(\varepsilon_r-\mu) / k_B T\gg 0$ to be fulfilled. What is going on?
 A: You are absolutely right that the limit in which this approximation holds is
$$\beta(\epsilon - \mu) \gg 1 \,,$$
which is not trivially the 'high-temperature limit', and indeed looks rather like the low temperature limit. However, it also looks like the limit of large negative $\mu$. If we want to know how temperature will affect the exponent, we need to know how temperature will effect the chemical potential. To proceed, suppose we're dealing with a gas of non-interacting particles. The grand potential is, in this limit,
$$ \Phi = -k_B T \ln \mathcal{Z} = -k_B T \int_0^\infty \ln \mathcal{Z}_\epsilon  \,g(\epsilon)\,\mathrm{d}\epsilon \simeq -k_B T \int_0^\infty \ln \bigg(1 + \exp(-\beta(\epsilon - \mu))\bigg)\,g(\epsilon) \,\mathrm{d}\epsilon \,,$$
where $\mathcal{Z}_\epsilon$ is the grand partition function associated with the energy level $\epsilon$ and $g(\epsilon)$ is the density of states. The integral is essentially just a sum of the partition functions due to each energy level. To get to the final expression we have assumed that we can approximate the grand partition function like so:
$$ \mathcal{Z}_\epsilon = \sum_{n} \bigg(\exp(-\beta(\epsilon - \mu))\bigg)^n \simeq 1 +  \exp(-\beta(\epsilon - \mu)) \,,$$
which corresponds to the limit stated at the top. As a brief detour, if we want to find the average occupancy of the energy level $\epsilon$, we can use
$$ \langle N_\epsilon \rangle = -\left(\frac{\partial \Phi_\epsilon}{\partial \mu}\right)_{T,V} \simeq \exp(-\beta(\epsilon - \mu))\qquad \mathrm{where} \qquad \Phi_\epsilon = -k_B T \ln \mathcal{Z}_\epsilon\,,$$
which is the Maxwell-Boltzmann distribution we were expecting (in the second equality we have Taylor expanded the logarithm in accordance with $\beta(\epsilon - \mu) \gg 1$). Now the density of states for a three-dimensional gas in a box can be obtained by standard means --- I won't bother going through it here, but the end result is:
$$ \Phi = -k_B TV\left(\frac{mk_B T}{2 \pi \hbar^2}\right)^{3/2} \exp(\beta \mu) \equiv  -\frac{k_B T V}{\lambda^3} \exp(\beta \mu) \,,$$
where the thermal wavelength $\lambda$ has been defined appropriately. From here we can write
$$N_\mathrm{tot} \equiv N = -\left(\frac{\partial \Phi}{\partial \mu}\right)_{T,V} =  \frac{V}{\lambda^3} \exp(\beta \mu) \,,$$
and hence
$$ \boxed{\mu = k_B T \ln \left(\frac{N \lambda^3}{V}\right) \,.}$$
Now to answer your question. The condition at the top can be considered the limit of $\beta \mu$ being large and negative. We see from the above that
$$ \beta \mu = \ln\left(\frac{N \lambda^3}{V}\right) \qquad \mathrm{where} \qquad \lambda = \left( \frac{2 \pi \hbar^2}{mk_B T}\right)^{1/2} \,.$$
This quantity will be large and negative when the argument of the logarithm is small. This will be the case for a) low densities $N/V$, b) high temperatures $T$ and/or c) high-mass particles.
You should think of the underlying situation in which the classical limit holds as when the number of thermally accessible states vastly exceeds the number of particles. This is because under such circumstances we can ignore multiple occupation of energy levels, which means we can ignore the fine details of particle indistinguishability. In the canonical distribution, when the number of states vastly exceeds the number of particles we can account for indistinguishability with a simple (but approximate) correction of division of the partition function by $N!$ --- we must do this even in the classical case, otherwise we run into all sorts of problems like the Gibbs paradox. However, when states start to become multiply occupied, this simple prescription fails, and we need to be more sophisticated in our consideration of particle indistinguishability.
If you imagine our gas particles as being wavepackets with a width of $\lambda$ as defined above, then you can think of each particle as occupying a volume $\lambda^3$. This has a nice interpretation --- the quantity $N \lambda^3/V$ that appears in the expression for the chemical potential can be thought of as the fraction of space occupied by the particles. The classical limit corresponds to this quantity being small, so that it's very unlikely for two particles to be in the same place --- i.e., be in the same state (here I'm essentially considering the states of our system to be position eigenstates rather than the usual energy eigenstates). If this quantity becomes larger, we start to get 'multiple-occupation', and so we imagine our classical approximation will break down. This is consistent: when $N \lambda^3 /V \sim 1$, the argument of the logarithm in the chemical potential is no longer large and negative, and so indeed the condition at the very top of this page breaks down.
Hope this helps!
A: The confusion arises because there are two kinds of classical limits, depending of the system under study.
Let's start with fermions, which distribution is $n_F(\epsilon)=\frac{1}{e^{(\epsilon-\mu)/T}+1}$. The first classical limit (corresponding to the case mentioned in the question) is $T\gg \epsilon-\mu$. This corresponds to the case where the temperature is large compared to all energy scales in the problem, and we get $n_F(\epsilon)=\frac{1}{2}$, which makes sense since one can only have one or zero fermion in any state, so at very high temperature, when all state are equiprobable, we get $1/2$. However, this case is not the one for which you would expect to recover the Boltzmann statistics, because for a gas of free fermion, the kinetic energy $\epsilon_k$ is not capped (and therefore $T\gg \epsilon$ is not possible).
The usual classical gas correspond to the limit $\mu<0$ and $|\mu|\gg T$, which gives back the expected distribution, and a very dilute gas. One shows that this is indeed this limit which gives back the usual classical results of Maxwell (for example, the virial expansion, etc.).
Now bosons, with distribution $n_B(\epsilon)=\frac{1}{e^{(\epsilon-\mu)/T}-1}$. The Maxwell classical gas is obtain in the same limit than for fermions, where the $\pm 1$ does not count.
Wikipedia quotes the limit $T$ large as classical for the following reason : for photons, $\mu=0$, and the distribution is thus $n_B(\epsilon)=\frac{1}{e^{\epsilon/T}-1}$. In the limit $\gg \epsilon$, one gets $n_B(\epsilon)=\frac{T}{\epsilon}$ which is the classical distribution in energy of light, and gives, for instance, the Rayleigh-Jeans law. This is, of course, not the classical distribution of point-like particle in a gas.
A: The introductory paragraph you quote with horror says temperature ''high enough'' to avoid quantum effects.  (It did not say anything like ''arbitrarily large''.)  If the temperature is too low, things like Bose--Einstein condensation can occur, which invalidate Maxwell--Boltzmann statistics.  The temperature should be high enough so that it is unlikely to have a quantum effect, but not so high that pair production occurs (yet another quantum effect).  These conditions have nothing to do with your analysis of the validity of dropping the plus or minus one in the denominator, which is yet another condition for the validity of the Maxwell--Boltzmann statistics.  Another condition is that the interaction between the particles should be weak: these are all independent conditions.
Boltzmann's constant is rather small by macroscopic standards: $1.3806488 \times 10^{-23} \mathrm{m^2} \mathrm{kg/s^2 K}$ so you can see that $T$ would have to be enormous before it would make the quantity you are worried about, $(\epsilon−\mu)/\mathrm{kT}$, much less than 22.  
A: Maxwell Boltmann statistics can be derived using canonical ensemble or microcanonical ensemble or grandcanonical ensemble. And also Fermi Dirac statistics as well as Bose Einstein statistics can also be derived from any one of the microcanonical, canonical or grandcanonical ensembles
