The Maxwell-Boltzmann distribution is derived by maximizing the entropy of a ##N##-body state in the (grand-)canonical ensemble of non-interacting particles in a box of volume ##V##. A convenient single-particle basis is chosen in the previous posting, namely momentum space. Let's consider scalar particles for convenience. For the many-body system we use the corresponding Fock (occupation-number) basis of the common eigenvectors of the number occupying each single-particle state labeled by the complete set of discrete momenta ##\vec{p} \in 2 \pi/L \mathbb{Z}^3## (I make the simplification to set ##L_1=L_2=L_3=L##, i.e., use a cube as the "quantization volume").
Instead of evaluating the simple partition sum of the grand-canonical ensemble, we rather calculate the functional
$$Z[\beta(\vec{p})]=\mathrm{Tr} \exp\left \{-\sum_{\vec{p}} [\beta(\vec{p}) E(\vec{p}) -\alpha] \hat{N}(\vec{p}) \right \}.$$
At the end we get the usual partition sum by setting ##\beta=\frac{1}{k_{\text{B}} T}## and ##\alpha(\vec{\beta})=\beta \mu##, where ##T## is the absolute temperature, ##\mu## the chemical potential for the conserved particle number and ##E(\vec{p})=\vec{p}^2/(2m)##. Since the ##\hat{N}(\vec{p})## all commute with each other we have
$$Z[\alpha(\vec{p})]=\prod_{\vec{p}} \sum_{N(\vec{p})=0}^{\infty} \exp\{[-\beta E(\vec{p})+\alpha(\vec{p})]N(\vec{p}) \} = \prod_{\vec{p}} \frac{1}{1-\exp\{[-\beta E(\vec{p})+\alpha(\vec{p})]\}}.$$
We have made use of the fact that spin-0 particles are bosons, and thus each ##N(\vec{p}) \in \mathbb{N}_0##.
Now it's easier to use
$$W[\alpha(\vec{p})]=\ln Z[\alpha(\vec{p})]=-\sum_{\vec{p}} \ln \{1-\exp\{[-\beta E(\vec{p})+\alpha(\vec{p})]\}.$$
The mean number of particles in the state ##\vec{p}## is given by
$$n(\vec{p}) = \langle \hat{N}(\vec{p}) \rangle=\left [\frac{\partial W}{\partial \alpha(\vec{p})} \right]_{\alpha(\vec{p})=\beta \mu}=\frac{\exp [-\beta E(\vec{p})+\beta \mu]}{1-\exp [-\beta E(\vec{p})+\beta \mu]}=\frac{1}{\exp [\beta E(\vec{p})-\beta \mu]-1}.$$
Now making ##L## very large we can take ##\vec{p} \in \mathbb{R}^3## making use of the above derived single-particle density we find
$$\mathrm{d} N=f(\vec{p}) \mathrm{d}^3 \vec{p}= f_{\text{B}}[E(\vec{p})] V \frac{\mathrm{d}^3 \vec{p}}{(2 \pi \hbar)^3}.$$
For a very diluted gas, for which ##\beta \mu## is a very small negative number, we can neglect the 1 in the denominator of the Bose-Einstein distribution, and we get the classical limit,
$$f(\vec{p})=\frac{V}{(2 \pi \hbar)^3} exp[-\beta E(\vec{p})+\beta \mu].$$
This is the phase-space distribution function of the ideal gas.
Now to get the distribution function for ##v=|\vec{p}|/m=p/m## we note that ##\mathrm{d}^3 \vec{p}=\mathrm{d} p \mathrm{d}^2 \Omega p^2=\mathrm{d} v \mathrm{d}^2 \Omega m^3 v^2## and thus
$$\mathrm{d} N=f(\vec{p}) \mathrm{d}^3 \vec{p} =\frac{V}{(2 \pi \hbar)^3} \exp(-\beta m v^2/2+\alpha) m^3 v^2 \mathrm{d} p \mathrm{d}^2 \Omega$$
Integrating over the whole solid angle and using that the integrand doesn't depend on the angles we simply get
$$\frac{\mathrm{d} N}{\mathrm{d} v} = \frac{V}{(2 \pi \hbar)^3} 4 \pi m^3 v^2 \exp(-\beta m v^2/2+\alpha).$$
You have to be aware that if you change the variable of interest for a distribution, you always get the original distribution multiplied with the appropriate Jacobian changing the "volume element" from one set of variables to another.
You can also check that the energy-distribution function is
$$\frac{\mathrm{d} N}{\mathrm{d} E} \propto \sqrt{E} \exp(-\beta E+\alpha)$$
using
$$E=m v^2/2 \; \Rightarrow \; v=\sqrt{2 m E} \; \Rightarrow \; \mathrm{d} v=\mathrm{d} E \sqrt{2 m}\frac{1}{2 \sqrt{E}}=\mathrm{d} E \sqrt{\frac{m}{E}}.$$