On the statistical meaning of density of states (DOS) According to the so-called law of the unconscious statistician:

The expected value $\langle \cdot \rangle$ of a measurable function of ${\displaystyle X}$, ${\displaystyle g(X)}$, given that $X$ has a probability density function ${\displaystyle f(x)}$, is given by the inner product of $f$ and ${\displaystyle g}$
$${\displaystyle \langle g(X)\rangle =\int _{\mathbb {R} }g(x)f(x)\,dx}$$

Suppose now I have an ensemble of $N$ non-interacting particles at equilibrium temperature $T$. The average energy of the system can be calculated as follows:
$${\displaystyle U =\int E\,f(E)\,g(E)\,{\rm {d}}E}$$
while
$$N =\int f(E)\,g(E)\,{\rm {d}}E$$
here $f(E)$ represents the probability density function (PDF) (e.g. Fermi-Dirac PDF, Bose-Einstein PDF), whereas $g(E)$ is the density of states.
My question is: why doesn't $f(E)$ alone play the role of $f(x)$ in the statistical definition? Why is $g(E)$ necessary too?
 A: The reason is that $f(E)$ is not a probability density function (this is a misinterpretation of some sloppy sources), but it’s the mean occupation number of a state, that is, it’s the mean number of one-particle states that occupy a level with energy $E$.
Therefore, since there can be $g(E)$ degenerate levels with energy $E$, the average number of particles between $E$ and $E+\mathrm{d}E$ is $g(E)f(E)\mathrm{d}E$.
A: Average value of a variable $x_i$ having probability $P(x_i)$ for $N$ variables can be given as,$$\bar x=\sum_{i=1}^Nx_iP(x_i)=\sum_{i=1}^Nx_i\frac{n_i}{N}$$where $n_i$ is number of  variables having value $x_i$. The total value of variables $X$, is either given by summing of all values of all variables or multiply number of variables to average value,$$X=N\bar x=\sum_{i=1}^Nx_iNP(x_i)$$Now if $N$ is break into summation and write $p(x_i)=\frac{P(x_i)}{N}$. Generally summation of $N$ is given in terms of either $n_i$ or $x_i$. So $X$ is given in terms of,$$X=\sum_{i=1}^NN^2x_ip(x_i)$$where $N^2x_i=g(x)dx$ and $p(x_i)=f(x)$.
But Planck calculated the average energy of cavity, so naturally multiply with number of modes or oscilator to have total energy. But the value came so small, so multiply with weighted function of $g(E)$ which is number of modes per unit volume of cavity. It is clear that probability density is multiplied with $N$. It is not necessary that to have total value of variables, $g(x)$ must be employed.
In case of Planck's oscillator, probably his idea was multiply energy term with square of modes but that gave only average energy because that was divided again by number of modes. If simply add the number of oscillators multiplied with minimum energy of mode, that gives total energy of cavity, where $nh\nu$ is value of energy for $n^{th}$ level and $\frac{n}{N}=\exp{\frac{-nh\nu}{kT}}$ is probability of any mode to occupy that.$$\begin{aligned}E&=h\nu+2h\nu\ldots=\sum_{i=0}^{\infty}n_ih\nu\\&=Nh\nu\left(1+\exp{\frac{-h\nu}{kT}}+\exp{\frac{-2h\nu}{kT}}+\ldots\right)\\&=Nh\nu\frac{1}{\exp{\frac{h\nu}{kT}}-1}\Rightarrow \frac{E}{N}=\bar E=\frac{h\nu}{\exp{\frac{h\nu}{kT}}-1}\end{aligned}$$Above expression shows that this function is similar to final expression by Planck, but that one is upscaled by again square of number of modes and variable of integration. This weighted up the function. Planck thought that this gave average value, but it gives only average value of particular energy level.
So density of states or number of states multiplied because to equate or upscale theoretical value. Otherwise no need to calculate because in Planck's calculation average value multiplied with number of modes, gives total energy, but to reuse number of states, they made it average energy.
