If it is still actual for you. The probability density
$$
e^{-\beta H\left( \mu\right) }
$$
is the so called canonical (Gibbs) distribution.There are plenty of methods how
to derive it. I can reproduce the simplest one.
Let's imagine that your system has the Hamiltonian $H\left( \mu\right) $ and you would like to study it for a certain temperature. In order to set the temperature you put your system inside a thermostat, so that your system exchanges only energy with the
thermostat but the volume and the number of particles are constant. Let's
suppose that the thermostat is a big tank filled by an ideal gas so that its
energy:
$$
h=\sum_{i=1}^{N}\frac{P_{i}^{2}}{2m}.
$$
The total system (your system+thermostat) is isolated thus the total energy is
fixed. Therefore, the distribution with respect to the total energy is a
delta-function:
$$
\rho\left( E\right) =\Lambda\delta\left( h+H-E\right) ,
$$
where $\Lambda$ is a some normalization factor such that
$$
\int \rho\left( E\right)\,d\Gamma =1,\qquad\left( 1\right)
$$
where $d\Gamma$ is an element of the full phase space:
$$
d\Gamma=d\mu\prod_{i=1}^{N}d^{3}P_{i}\,d^{3}Q_{i}.
$$
Let's integrate out all degrees of freedom of the thermostat:
$$
\rho\left( H\right) =\int\prod_{i=1}^{N}d^{3}P_{i}\,d^{3}Q_{i}
\,\Lambda\delta\left( H+\sum_{i=1}^{N}\frac{P_{i}^{2}}{2m}-E\right) =\Lambda
V^{N}\int\prod_{i=1}^{N}d^{3}P_{i}\,\,\delta\left( H+\sum_{i=1}^{N}
\frac{P_{i}^{2}}{2m}-E\right) ,
$$
where $V$ is the volume of thermostat. The integration measure can be
simplified as follows:
$$
\int\left[\prod_{i=1}^{N}d^{3}P_{i}\right] f(\epsilon)=\frac{2\pi^{3N/2}}{\Gamma\left( 3N/2\right)
}\int\left[\epsilon^{3N-1}d\epsilon\right] f(\epsilon),
$$
where
$$
\epsilon^{2}=\sum_{i=1}^{N}P_{i}^{2}.
$$
Hence the integration can be preformed as follows:
$$
\rho\left( H\right) =\Lambda V^{N}\frac{2\pi^{3N/2}}{\Gamma\left(
3N/2\right) }\int\,d\epsilon\,\epsilon^{3N-1}\,\delta\left( H+\frac
{\epsilon^{2}}{2m}-E\right) =\Lambda V^{N}\frac{2m\pi^{3N/2}}{\Gamma\left(
3N/2\right) }\,\left( E-H\right) ^{\frac{3N}{2}-1}.
$$
Let's now consider the $N\rightarrow\infty$ limit so that
$$
\frac{E}{N}\approx\frac{h}{N}=\frac{3T}{2}.
$$
The distribution takes the form:
$$
\rho\left( H\right) \sim\left( 1-\frac{H}{E}\right) ^{3N-2}\approx\left(
1-\frac{H}{\frac{3N}{2}T}\right) ^{\frac{3N}{2}-1}\approx\exp\left(
-\frac{H}{T}\right) .
$$
The normalization factor can be found from the normalization condition (1).
Finally the probability density for the energy of your system takes the form:
$$
\rho\left( H\right) =\frac{e^{-\beta H\left( \mu\right) }}{Z},\quad Z=\int
d\mu\,e^{-\beta H\left( \mu\right) }.
$$
In fact, the result is independent of the nature of the thermostat, see e.g. L.D. Landau, E.M. Lifshitz, Volume 5 of Course of Theoretical Physics, Statistical physics Part 1 Ch. III, The Gibbs distribution.
